We present a novel numerical algorithm for the simulation of seismic
wave propagation in porous media, which is particularly suitable
for the accurate modelling of surface wave-type phenomena. The differential
equations of motion are based on Biot's theory of poro-elasticity
and solved with a pseudospectral approach using Fourier and Chebyshev
methods to compute the spatial derivatives along the horizontal and
vertical directions, respectively. The time solver is a splitting
algorithm that accounts for the stiffness of the differential equations.
Due to the Chebyshev operator the grid spacing in the vertical direction
is non-uniform and characterized by a denser spatial sampling in
the vicinity of interfaces, which allows for a numerically stable
and accurate evaluation of higher order surface wave modes. We stretch
the grid in the vertical direction to increase the minimum grid spacing
and reduce the computational cost. The free-surface boundary conditions
are implemented with a characteristics approach, where the characteristic
variables are evaluated at zero viscosity. The same procedure is
used to model seismic wave propagation at the interface between a
fluid and porous medium. In this case, each medium is represented
by a different grid and the two grids are combined through a domain-decomposition
method. This wavefield decomposition method accounts for the discontinuity
of variables and is crucial for an accurate interface treatment.
We simulate seismic wave propagation with open-pore and sealed-pore
boundary conditions and verify the validity and accuracy of the algorithm
by comparing the numerical simulations to analytical solutions based
on zero viscosity obtained with the Cagniard-de Hoop method. Finally,
we illustrate the suitability of our algorithm for more complex models
of porous media involving viscous pore fluids and strongly heterogeneous
distributions of the elastic and hydraulic material properties.