Abstract. The scattered seismic waves of fractured porous rock are strongly affected by the wave-induced fluid pressure diffusion effects between the compliant fractures and the stiffer embedding background. To include these poroelastic effects in seismic modeling, we develop a numerical scheme for discrete distributed large-scale fractures embedded in fluid-saturated porous rock. Using Coates and Schoenberg’s local effective medium theory and Barbosa’s dynamic linear slip model characterized by complex-valued and frequency-dependent fracture compliances, we derive the effective viscoelastic compliances in each spatial discretized cell by superimposing the compliances of the background and the fractures. The effective governing equations of the fractured porous rock are then characterized by the derived anisotropic, complex-valued, and frequency-dependent effective compliances. We numerically solved the effective governing equations by mixed-grid stencil frequency-domain finite-difference method. The good consistency between the scattered waves off a single horizontal fracture calculated using our proposed scheme and those calculated using the poroelastic linear slip model shows that our modeling scheme can properly include the FPD effects. We also find that for a P-point source, the amplitudes of the scattered waves from a single horizontal fracture are strongly affected by the fluid stiffening effects due to fluid pressure diffusion, while for an S-point source, the scattered waves are less sensitive to fluid pressure diffusion. In the case of the conjugate fracture system, the scattered waves from the bottom of the fractured reservoir and the reflected waves from the underlying formation are attenuated and dispersed by the FPD effects for both P- and S-point sources. The proposed numerical modeling scheme can also be used to improve migration quality and the estimation of fracture mechanical characteristics in inversion.