Although purely analytical models can provide rough qualitative predictions in the field of electrohydrodynamics (EHD), more sophisticated numerical approaches are necessary to quantitatively study the involved phenomena [1, 2]. Considering the computational cost and complexities associated with the mesh-free numerical methods, mesh-based methods are usually more efficient for fluid dynamics applications. Nevertheless, in the case of multiphase EHD flows, the difference in the material properties of the phases imposes discontinuities in the field variables (e.g. pressure, electric field). In this sense, the accuracy of the solution of EHD problems depends on the sharp representation of the strong (jump) and weak discontinuities in the field variables. So far, different numerical techniques have been proposed in the literature to represent such discontinuity, for example, Weighted Harmonic Averaging Method (WHAM), the Ghost Fluid Method (GFM), Immersed Interface Method (IIM), to name just a few. However, these schemes can accurately capture the electric field only in cases of small permittivity ratio, perfect dielectric fluids, or via a computationally expensive refinement process. On the other hand, the Enriched Finite Element Method (EFEM) can be acquired as a viable option for EHD problems. EFEM relies on the enrichment of the shape functions for the elements cut by the phase interface. In this work, such enrichment is proposed to accurately capture the weak discontinuity in the electric potential (or equivalently the jump in the electric field), adopting the ideas previously explored for representing pressure gradient discontinuity in two-phase flows [3]. The main advantage of this method is that the enrichment functions do not depend on the neighboring elements, and therefore, the associated additional degrees of freedom (DoF) can be condensed at the elemental level. This feature makes EFEM one of the most efficient techniques for multi-phase problems. Although this technique has been widely used for multi-phase CFD applications [4, 5], the employment of EFEM for EHD applications has scarcely been addressed in the literature. In this sense, the present work is among the very first applications of the EFEM method to EHD problems.