Electrical resistivity tomography, as a non-invasive geophysical technique, has the potential to provide vast information about the heterogeneity, saturation distribution and contamination in the porous media. However, the sole use of this technique has its own limitations including the ambiguity regarding the interpretation of the results. One solution to this problem would be to pair field measured resistivity profiles with numerical simulations. Numerical simulations can aslo act as tools for calibration and to facilitate parameter estimation. Therefore, there is a need for efficient numerical models that can handle coupled fluid flow and electrical current in highly heterogeneous porous media, notably fractured porous media. We have developed an advanced scheme for the coupled modeling of electrical current distribution and variably saturated fluid flow in a fractured porous media with a hybrid dimension approach in a planar domain (i.e. 1D fractures embedded in 2D porous matrix). Mixed hybrid finite element method, known to be highly efficient in handling local heterogeneity and producing accurate velocity fields, has been used as the numerical method for discretization of both fluid flow (i.e. Richard's equation) and electrical current (i.e. ensemble of Ohm's law and Helmholtz equation). The results of the numerical solution is then validated by comparing against the results obtained by an FEM-based package (i.e. COMSOL Multiphysics®) for three synthetic configurations of which one is non-fractured, one has a single inclined fracture and the third having a network of orthogonal partially connected fractures. The developed tool is then used to investigate the effect of some geometric characteristics (e.g. angle and length) of the fracture on the normalized apparent resistivity in the top midpoint of the domain for the single fracture configurations. Finally, the effect of electrical dipole location and inclination on the simulated normalized apparent resistivity is investigated numerically for two configurations assimilating the cross-hole dipole-dipole configurations. The results show that the orientation and the length of the fracture and the geometrical features of the electrical dipole play important roles in the accuracy and resolution of the simulated electrical response attributed to the fracture. The present work can be used as a new tool for the coupled simulation of electrical current and variably saturated flow in highly heterogeneous porous media and the present scheme can be extended with facility to handle more realistic geometries.