To approximate seismic wave propagation in double-porosity media, one can use the effective Biot theory, which can explain the high level of attenuation observed at seismic frequencies, but which is unaccounted for with classic Biot theory. The governing equations of the effective Biot theory for homogeneous media with complex frequency dispersion characteristics need to be upscaled and extended to S-waves using a poroviscoelastic approach when simulating heterogeneous porous media. We have applied a frequency-space domain mixed grid finite-difference modeling method for this purpose and computed the wavefields for solid particle velocity, fluid flux, and pore pressure. A homogeneous full space model with a single Cole-Cole relaxation function (fractional Zener model) representing the attenuation mechanism is used to indicate that the numerical solutions match fairly well the analytical waveform solutions of the effective Biot theory for source center frequencies ranging from 25 to 10,000 Hz. The computed pore pressure wavefield for a two-layer example porous model is used to further support our numerical method. This model and a laterally heterogeneous three-layer porous medium model further illustrate the applicability of the procedure and even indicate the presence of the slow Biot wave near the interfaces as a result of mode conversion on reflection and transmission. The nearly perfectly matched layer (NPML) absorbing boundary condition is chosen to truncate the computational grid and avoid reflections from the model edges. We found and proved that the NPML performs identically to the PML for either elastic wave modeling or Biot porous medium wave modeling for 3D and 2D scenarios. Two methods to adjust the maximum damping factor (one for the source center frequency and the other for each frequency component) are suggested to adjust the damping factor function in the NPML condition to ensure its effectiveness for various frequency ranges.