We address a 2D finite difference (FD) frequency-domain modeling algorithm based on the theory of fractional diffusion of electromagnetic (EM) fields, which is generated by an infinite line source lying above a fractured geological medium. The presence of fractures in the subsurface, usually containing highly conductive pore fluids, gives rise to spatially hierarchical flow paths of induced EM eddy currents. The diffusion of EM eddy currents in such formations is anomalous, generalizing the classical Gaussian process described by the conventional Maxwell equations. Based on the continuous time random walk (CTRW) theory, the diffusion of EM eddy currents in a rough medium is governed by the fractional Maxwell equations. Here, we model the EM response of a 2D subsurface containing fractured zones, based on the fractional Maxwell equations. The governing equation in the frequency domain is discretized using the FD approach. The resulting equation system is solved by the multifrontal massively parallel solver (MUMPS). We find excellent agreement between the FD and analytic solutions for a rough half-space model. Then, FD solutions are calculated for a 2D fault zone model with variable conductivity and roughness. We illustrate a case in which a rough fault zone would not be resolved by classical diffusion modeling, even if its conductivity contrasts with the background.