Numerical modeling of seismic waves in heterogeneous, porous reservoir rocks is an important tool for interpreting seismic surveys in reservoir engineering. Various theoretical studies derive effective elastic moduli and seismic attributes from complex rock properties, involving patchy saturation and fractured media. To confirm and further develop rock-physics theories for reservoir rocks, accurate numerical modeling tools are required. Our 2D velocity-stress, finite-difference scheme simulates waves within poroelastic media as described by Biot's theory. The scheme is second order in time, contains high-order spatial derivative operators, and is parallelized using the domain-decomposition technique. A series of numerical experiments that are compared to exact analytical solutions allow us to assess the stability conditions and dispersion relations of the explicit poroelastic finite-differ-ence method. The focus of the experiments is to model wave-induced flow accurately in the vicinity of mesoscopic heterogeneities such as cracks and gas inclusions in partially saturated rocks. For that purpose, a suitable numerical setup is applied to extract seismic attenuation and dispersion from quasi-static experiments. Our results confirm that finite-difference modeling is a valuable tool to simulate wave propa-gation in heterogeneous poroelastic media, provided the temporal and spatial scales of the propagating waves and of the induced fluid-diffusion processes are resolved properly.