Viscoelastic representation of surface waves in patchy saturated poroelastic media

Abstract

Wave-induced flow is observed as the dominated factor for P wave propagation at seismic frequencies. This mechanism has a mesoscopic scale nature. The inhomogeneous unsaturated patches are regarded larger than the pore size, but smaller than the wavelength. Surface wave, e.g., Rayleigh wave, which propagates along the free surface, generated by the interfering of body waves is also affected by the mesoscopic loss mechanisms. Recent studies have reported that the effect of the wave-induced flow in wave propagation shows a relaxation behavior. Viscoelastic equivalent relaxation function associated with the wave mode can describe the kinetic nature of the attenuation. In this paper, the equivalent viscoelastic relaxation functions are extended to take into account the free surface for the Rayleigh surface wave propagation in patchy saturated poroelastic media. Numerical results for the frequency-dependent velocity and attenuation and the time-dependent dynamical responses for the equivalent Rayleigh surface wave propagation along an interface between vacuum and patchy saturated porous media are reported in the low-frequency range (0.1–1,000 Hz). The results show that the dispersion and attenuation and kinetic characteristics of the mesoscopic loss effect for the surface wave can be effectively represented in the equivalent viscoelastic media. The simulation of surface wave propagation within mesoscopic patches requires solving Biot’s differential equations in very small grid spaces, involving the conversion of the fast P wave energy diffusion into the Biot slow wave. This procedure requires a very large amount of computer consumption. An efficient equivalent approach for this patchy saturated poroelastic media shows a more convenient way to solve the single phase viscoelastic differential equations.

1 Introduction

For usual subsurface unconsolidated sedimentary media, the presence of pore fluid can dramatically influence the wave propagation in a porous medium. The bulk modulus of the medium becomes frequency-dependent and attenuation effects arise due to the distributions of the fluid in pore space. It is particularly interesting to consider the problem of a gas–liquid mixture saturated medium. In this case, even more dissipative mechanisms have to be considered beyond the global interaction between the pore fluid and the solid frame in case of pore space only saturated by liquid in terms of Biot’s theory (Biot 1956). Recent studies (Pride et al. 2004) have shown that the main cause of attenuation at low frequencies (seismic frequencies) is patchy saturation at the mesoscopic scale, larger than the pore size but smaller than the wavelength (typically tens of centimeters), in which the fluid distribution scale is large enough so that the wave-induced pore pressure gradient cannot equilibrate during a wave period. Patches of nonuniform saturation always occur at the gas–liquid contacts. Between full gas and full liquid saturation, typically a transition zone exists. During the underground liquid motion, the pressure variation might lead to the diffusive formation of free gas pockets (Johnson 2001).

This low-frequency loss effect of partially gas-saturation has been taken into account for many years. The mesoscopic effects are firstly modeled by isolated spherical gas patched in the liquid saturated background by White (1975). Since then, significant progress has been made by considering various patch distributions and flow regimes (e.g., White et al. 1975; Dutta and Odé 1979a, b; Dutta and Seriff 1979; Norris 1993; Gelinsky et al. 1998; Johnson 2001; Müller and Gurevich 2004, 2005; Müller et al. 2008; Gurevich et al. 2009; Vogelaar et al. 2010). More references for this topic are found in Toms et al. (2006). The gas pocket is subjected to the pore pressure variation of waves, whose wavelength is larger than the size of the inhomogeneity, and perturbed by the oscillations. The generated waves on this scale induce the fluid flow between the pocket and background matrix, which causes the intrinsic attenuation. The real-data, laboratory observations, and numerical simulations are consistent with a partially saturated reservoir model in which poroelastic effects caused by wave-induced fluid flow of different fluid phases can modify the omnipresent wave background behavior (e.g., Winkler and Nur 1979; Knight et al. 1998; Carcione et al. 2003; Wenzlau and Müller 2009; Saenger et al. 2009).

To investigate the mesoscale attenuation effect, we have to solve the differential equations in very small grids, which always faces large computational consumption and unstable problems (Picotti et al. 2007). An equivalent approach has to be suggested. For global homogeneous saturation of the Biot model, it was reported that the fast waves can be effectively approximated by a linear viscoelastic relaxation models (e.g., Geertsma and Smit 1961; Ben-Menahem and Singh 1981; Carcione and Quiroga-Goode 1996; Carcione 1998, 2007). This approximate method is also extended to describe the wave scattering problems (Morochnik and Bardet 1996). Recent studies also show that the mesoscopic saturation effect on body waves can be equivalently represented by viscoelastic models (Liu et al. 2009a, 2010; Picotti et al. 2010, 2012). Despite all the efforts and attention to find effective representation of mesoscopic saturation for wave propagation, to our best knowledge, there is no study concerning the equivalent approach to represent the surface wave propagation and attenuation along an interface. Although great progress has been made by considering different boundaries (e.g., Zhang et al. 2011, 2012; Feng and Johnson 1983a, b; Gubaidullin et al. 2004; Allard et al. 2004; Markov 2009a), and multi-fluid distributions (e.g., Dai et al. 2006; Chao et al. 2006; Lo 2008) on surface waves. Markov et al. (2010) also pointed out that at low frequencies, zero-order approximation of the surface wave propagation corresponds to an elastic medium, as the body waves in Gassmann equivalent medium. In order to further apply poroelasticity in the practical data processing in an economic way, it is necessary to find an effective way to represent the attenuation of the surface wave in the poroelastic media, especially by considering the dominant mesoscopic loss mechanism in the seismic frequency band. In this paper, we concentrate on the low-frequency effective effects of the patchy saturated gas pockets on the near surface wave propagations. Therefore, the vacuum/patchy saturated poroelastic medium interface is adopted. This paper is organized as follows. In Sect. 2, we review the models for the wave propagation in patchy saturated media of gas–liquid mixture saturation. In Sect. 3, we introduce a linear viscoelastic relaxation functions to equivalently represent the attenuation of the surface wave. Some numerical examples of the frequency and time-dependent characteristics of the equivalent models are present and analyzed in Sect. 4. Finally, the discussions and conclusions are summarized in Sect. 5.

2 Bulk modulus of a patchy saturated medium

A patchy saturation model of a nonrigid porous medium fully saturated by a fluid that contains gas pockets (radius a) larger than the typical pore size is considered. The interaction among the individual gas pocket is neglected by defining a liquid influence shell (radius b) surrounding each pocket (Fig. 1). The gas friction is $S_1=\frac{a^3}{b^3}$, and the liquid friction is S₂ = 1 – S₁. The radius b is chosen for the volume of each sphere $\frac{4}{3}\mathrm{\pi }{b}^3$ equals to the volume of the unit cell of the cubic lattice.

Fig. 1

Geometry of a cubic lattice of periodic spherical gas pocket with radius a, separated by distance 2b′ (from Vogelaar et al. 2010). Each gas pocket is surrounded by a liquid shell with radius b, so that the volume of the cube equals to the volume of the sphere $V_{{\text{b}}^{\prime }}=V_{\text{b}}$

The external pressure field is assumed to be spatially homogeneous at the scale of the mesoscopic inhomogeneity. We presume the wave frequency f is low enough that the Biot theory (Biot 1956) is around its low-frequency limit, i.e., f ≪f_c, where the Biot critical frequency is defined by $f_{\text{c}}=\frac{\mathit{\varphi }\mathit{\eta }}{2\mathit{\pi }C\mathit{\kappa }{\mathit{\rho }}_{\text{f}}}$. Therefore, the whole medium is on its quasi static state. The physical properties about the rock frame are porosity $\mathit{\varphi }$, tortuosity C, and the permeability κ; those about the pore fluid are viscosity η and density ρ_f.

The starting conditions are essentially the Biot theory at quasi static state by setting the all higher order inertial terms to zero and by taking the dynamical permeability equal to its steady one (Vogelaar et al. 2010). The quasi static Biot equations are degenerated as

$$\nabla ·\mathit{\tau }=0,$$

(1)

$$\frac{\mathit{\kappa }}{\mathit{\eta }}\nabla p=-i\mathit{\omega }\mathit{w},$$

(2)

where w is the relative displacement of the fluid, defined as $\mathit{w}=\mathit{\varphi }(\mathit{U}-\mathit{u}).\mathit{u}$ and U, u are the solid and fluid displacements, respectively. The solid stress τ and the pore fluid pressure p are written by

$$\mathit{\tau }=2\mathit{\mu }{e}_{i,j}+{\mathit{\delta }}_{i,j}({\mathit{\lambda }}_{\text{c}}-\mathit{\alpha }M\mathit{\varsigma }),$$

(3)

$$p=-\mathit{\alpha }Me+M\mathit{\varsigma },$$

(4)

where $e_{i,j}=\frac{1}{2}({\mathit{\partial }}_j{u}_i+{\mathit{\partial }}_i{u}_j),e=\nabla ·\mathit{u}$ and $\mathit{\varsigma }=-\nabla ·\mathit{w}.{\mathit{\lambda }}_{\text{c}}$ and μ are lame coefficients. α and M are the Biot coefficients. Explicit expressions are given in terms of the bulk moduli of the pore fluid, the solid, and the matrix K_f,s,m, respectively

$$\mathit{\alpha }=1-\frac{K_{\text{m}}}{K_{\text{s}}},$$

(5)

$$\frac{1}{M}=\frac{\mathit{\alpha }-\mathit{\varphi }}{K_{\text{s}}}+\frac{\mathit{\varphi }}{K_{\text{f}}},$$

(6)

The Gassmann bulk modulus can be written by

$$K_{\text{G}i}=K_{\text{m}}+{\mathit{\alpha }}^2{M}_i,i=1,2,$$

(7)

where i denotes the different saturated regions in the patch for the gas and the liquid. At low-frequency range, the fluid pressure has enough time to equilibrate at a uniform state (relaxed state). The effective fluid bulk modulus is defined as the Reuss average $K_{\text{W}}={\left(\frac{S_1}{K_{\text{f1}}}+\frac{S_2}{K_{\text{f2}}}\right)}^{-1}.$ The effective modulus of the whole composite at the low-frequency limit, referred as the static limit, is given by

$$K_{\text{GW}}=K_{\text{m}}+{\mathit{\alpha }}^{2}M(K_{\text{W}}),$$

(8)

In this low-frequency limit, the slow wave is diffusive. The diffusivity is given by

$$D(K_{\text{f}})=\frac{\mathit{\kappa }M\left(K_{\text{m}}+\frac{4}{3}\mathit{\mu }\right)}{\mathit{\eta }\left(K_{\text{G}}+\frac{4}{3}\mathit{\mu }\right)},$$

(9)

which is also a solution to Eqs. (1) and (2) when the fluid and solid move out-of-phase (Johnson 2001).

On the other hand, when the frequency is sufficiently high, the fluid pressure cannot equilibrate, and this unrelaxed state induces fluid pressure gradient. The fluid pressure is not uniform, but it can be assumed to be constant within each phase. The composite bulk modulus can be given by Hill theorem, ignoring the multiphase flow effect, referred as no-flow limit

$$K_{\text{GH}}={\left(\frac{S_1}{K_{\text{G1}}+\frac{4}{3}\mathit{\mu }}+\frac{S_2}{K_{\text{G2}}+\frac{4}{3}\mathit{\mu }}\right)}^{-1}-\frac{4}{3}\mathit{\mu },$$

(10)

The macroscopic bulk modulus between these two frequency limits can be obtained by considering the representative volume comprising a single gas pocket and a liquid shell surrounding the pocket on the volume variation as a response to an external oscillating field. Solving Eqs. (1) and (2) yields the solid displacement as a function of the outer applied pressure p_e and hence the effective bulk modulus of the representative volume. The effective bulk modulus can be written by

$$K(\mathit{\omega })=-\frac{b}{3u_2}{p}_{\text{e}},$$

(11)

where u₂ is the complex-valued radial solid displacement at the outer boundary of the unit cell.

Once the effective bulk modulus is obtained, velocity $c={\left[\Re \left(\frac{k(\mathit{\omega })}{\mathit{\omega }}\right)\right]}^{-1},$ and attenuation $\mathit{\theta }={tan}^{-1}\left|\frac{\Im (k_{}^2(\mathit{\omega }))}{\Re (k_{}^2(\mathit{\omega }))}\right|={tan}^{-1}\left(\frac{1}{Q}\right)$ of the wave propagation are computed using the effective complex wavenumber:

$$k_1(\mathit{\omega })=\mathit{\omega }\sqrt{\mathit{\rho }/\left(K(\mathit{\omega })+\frac{4}{3}\mathit{\mu }\right)},$$

(12)

Subscript 1 describes the fast wave on the macroscale. $\mathit{\rho }=(1-\mathit{\varphi }){\mathit{\rho }}_{\text{s}}+\mathit{\varphi }((1-S_1){\mathit{\rho }}_{\text{f2}}+S_1{\mathit{\rho }}_{\text{f1}}),$ and ρ_f1,2 denotes the gas (subscript 1) or the liquid (subscript 2) density. We consider the approaches by White (1975) as corrected by Dutta and Seriff (1979), and Johnson (2001) to represent the complex-valued frequency-dependent crossover phenomenon, which are outlined in Appendix 1 in detail. It was shown by Tserkovnyak and Johnson (2003) that in case of considering the fluid surface tension across the different fluid patch interfaces, the transition behavior between the two limits will be affected and rescaled. However, this mechanism does not affect the analytical structure of the frequency-dependent bulk modulus. In general, the surface tension affects the slow P wave significantly which only perturbs the fast P wave slightly (also see Markov and Levin 2007; Markov 2009b). How this phenomenon affects the surface wave needs to be investigated further. It also should be noted that in the homogeneous solid frame background, the two frequency limits the shear wave motion coincident. Therefore, the model presented here assumes that the shear wave is influenced by the presence of the gas phase only due to changes in density, which is also demonstrated by numerical experiments (Rubino et al. 2009; Liu et al. 2009b). The frequency-dependent mechanisms incorporated in this model will be experimentally corroborated by impulse-induced vacuum/patchy saturated solid interface surface wave experiments.

3 Viscoelastic representation of a patchy saturated medium

It is well known that the velocities and attenuation of the fast waves in a poroelastic medium can be well fitted by the standard linear viscoelastic elements (Zener elements). The frequency-dependent modulus of a viscoelastic element of the Zener model can be expressed as

$$M(\mathit{\omega })=M_0\left(\frac{1+\text{i}\mathit{\omega }{\mathit{\tau }}_{\mathit{\epsilon }}}{1+\text{i}\mathit{\omega }{\mathit{\tau }}_{\mathit{\sigma }}}\right),$$

(13)

where τ _ε and τ _σ are relaxation time parameters for the stress and strain, respectively. M₀ is the relaxed modulus at low-frequency limit. The quality factor associated with M(ω) equals to $\Re (M)/\Im (M)$. Its minimum value Q₀ that represents the maximum damping locates at the frequency of ω₀ = 2πf₀. The corresponding relaxation time parameters can be written by

$${\mathit{\tau }}_{\mathit{\epsilon }}=\frac{1}{{\mathit{\omega }}_0{Q}_0}\left(\sqrt{Q_0^2+1}+1\right),$$

(14)

$${\mathit{\tau }}_{\mathit{\sigma }}=\frac{1}{{\mathit{\omega }}_0{Q}_0}\left(\sqrt{Q_0^2+1}-1\right),$$

(15)

The complex velocities in the viscoelastic media are related to relaxed modulus for two mechanisms for the two fast waves:

$${\tilde{v}}_{\text{p}}=v_{\text{p0}}\sqrt{M_1/M_{01}},$$

(16)

$${\tilde{v}}_{\text{s}}=v_{\text{s0}}\sqrt{M_2/M_{02}},$$

(17)

where v_p0 and v_s0 are quasi static limit phase velocities of the fast P wave and the shear wave, respectively.

For the patchy saturated poroelastic media, the mesoscopic loss mechanism triggered around the quasi static limit of wave propagation. It also means that the slow P wave is mainly at its diffusion character, which controls the wave-induced flow between the different fluid saturated patches. It can be thought that around the quasi static limit, the combinations of the propagation of the fast P wave and the diffusion of the slow P wave make the frequency-dependent effective bulk modulus of the representative patchy saturated volume at low frequencies. Therefore, the patchy saturated media have the nature that one dominant effective P wave can be propagatory, similar to the viscoelastic media. Therefore, it is possible to represent the P wave propagations in a patchy saturated medium by fitting the different wave modes by Eq. (16). We choose

$$M_{01}=K_{\text{GW}}+\frac{4}{3}\mathit{\mu },$$

(18)

and definite Q₀ and ω₀ corresponding to the peak location of the P wave attenuation to represent the wave propagation. Picotti et al. (2010) have showed the viscoelastic representations of the fast P wave in patchy saturated media. It is noted that they chose the high frequency unrelaxed modulus M_u to write Eq. (13), by $M(\mathit{\omega })=M_{\text{u}}\frac{{\mathit{\tau }}_{\mathit{\sigma }}}{{\mathit{\tau }}_{\mathit{\epsilon }}}\left(\frac{1+\text{i}\mathit{\omega }{\mathit{\tau }}_{\mathit{\epsilon }}}{1+\text{i}\mathit{\omega }{\mathit{\tau }}_{\mathit{\sigma }}}\right).$ Therefore, they obtained the bad fit of the velocity at low frequencies (Fig. 4 in Picotti et al. 2010). Due to the patchy saturation being accounted for the attenuation at low frequencies, it is believed that the chosen relaxed modulus to describe the transition behavior can gain better velocity, representing at low frequencies that we are interested in, which will be also illustrated by our numerical experiments in the next section.

In this research, the shear wave is assumed not to be perturbed by the fluid distributions. The propagation of the shear wave is described by the Biot model. Therefore,

$$M_{02}=\mathit{\mu },$$

(19)

and definite Q₀ and ω₀ corresponding to the attenuation peak of the global fluid flow.

In order to obtain the interface effect for the surface wave propagation, the free surface boundary conditions are introduced. For the effective viscoelastic medium in Cartesian coordinate,

$${\mathit{\tau }}_{zz}{|}_{z=0}=F(t)\mathit{\delta }(x)\mathit{\delta }(y)$$

(20)

$${\mathit{\tau }}_{xz}{|}_{z=0}=0$$

(21)

$${\mathit{\tau }}_{yz}{|}_{z=0}=0,$$

(22)

where δ is Dirac delta function. By the constraint of the boundary conditions, the equivalent representation of the Rayleigh surface wave can be obtained by the elastic wave modes of the P and shear waves. The mathematical procedure involving the numerical solution of the boundary problem at the interface shows in Appendix 2.

4 Numerical examples

In this section, the numerical results of the viscoelastic representation of the surface wave modes that propagate along a vacuum/patchy saturated half space interface are discussed. We consider a typical clay rock, whose physical properties are obtained from the soil sample near a highway construction site. The properties of the gas and liquid correspond to the air and the water near surface. The material properties of the clay and the fluids are given in Table 1.

Table 1 Material properties of a patchy saturated medium

Figure 2 shows the phase velocities (a) and loss angles (b) of the P waves in the White and Johnson patchy saturation models and the corresponding viscoelastic equivalent representations, for the patch scale of 0.4 m when the gas friction is S₁ = 0.5. The velocity and attenuation dispersion curves of the two patchy saturation models cases are quite similar. The best fits of the Zener model are obtained with parameters f₀ = 37.7 Hz and Q₀ = 9.6 for the White model, and f₀ = 42 Hz and Q₀ = 10.1 for the Johnson model. It is noted that the Zener fits for the both models are better in the loss angle curves than in the velocity curves. Moreover, the curves are more reliable at the seismic frequencies (<100 Hz). It is verified that our choice of the relaxed modulus at quasi static limit to represent the viscoelastic elements is more appropriate than the choice of the unrelaxed modulus at no flow limit in the analysis of Picotti et al. (2010) for the equivalent represent in the seismic frequency band. Figure 3 shows the phase velocities (a) and loss angles (b) for the shear wave in case of the gas friction S₁ = 0.5. Because the reference critical frequency f_c is far greater than 1000 Hz, the shear wave propagates around its quasi static limit with nondistinctive dispersion and small attenuation. At low frequencies, the Zener fits in velocity and attenuation are in great agreement.

Fig. 2

Fit of both the phase velocity (a) and loss angle (b) curves of the P waves in a mixture of water and air patchy saturated porous solid by using viscoelastic representation. The patchy saturation models are obtained by White (1975) as corrected by Dutta and Seriff (1979), and Johnson (2001). The material properties are shown in Table 1. S₁ = 0.5 and b = 0.4 m. Fitting parameters of the Zener elements are f₀ = 37.7 Hz, Q₀ = 9.6 for the White model, and f₀ = 42.7 Hz, Q₀ = 10.1 for the Johnson model

Fig. 3

Fit of both the phase velocity (a) and loss angle (b) curves of the shear wave in a mixture of water and air patchy saturated porous solid by using viscoelastic representation. The material properties are shown in Table 1. S₁ = 0.5. The shear wave propagates around its quasi static limit

Figure 4 shows the phase velocities (a) and loss angles (b) in case of Fig. 1, except for the gas friction S₁ = 0.1. The transition behaviors of the mesoscopic loss at low frequencies are present stronger and triggered at lower frequency than those for S₁ = 0.5. The best fits of the Zener model are obtained with the parameters f₀ = 4.37 Hz and Q₀ = 4.40 for the White model, and f₀ = 4.50 Hz and Q₀ = 4.41 for the Johnson model. Distinctive differences between the patchy saturation models and the viscoelastic equivalent models are observed at the high frequency range (>100 Hz). Figure 5 shows the phase velocities (a) and loss angles (b) for the shear wave in case of the gas friction S₁ = 0.1. The shear wave also propagates around its quasi static limit with weaker dispersion and attenuation. The Zener fits in velocity and attenuation are in great agreement as well.

Fig. 4

Fit of both the phase velocity (a) and loss angle (b) curves of the P waves in a mixture of water and air patchy saturated porous solid by using viscoelastic representation. The patchy saturation models are obtained by White (1975) as corrected by Dutta and Seriff (1979), and Johnson (2001). S₁ = 0.1 and b = 0.4 m. Fitting parameters of the Zener elements are f₀ = 4.37 Hz, Q₀ = 4.40 for the White model, and f₀ = 4.50 Hz, Q₀ = 4.41 for the Johnson model

Fig. 5

Fit of both the phase velocity (a) and loss angle (b) curves of the shear wave in a mixture of water and air patchy saturated porous solid by using viscoelastic representation. S₁ = 0.1. The shear wave propagates around its quasi static limit

Figure 6 shows the phase velocities (a) and loss angles (b) for the Rayleigh waves in the White and Johnson patchy saturation models and the corresponding viscoelastic equivalent representations, for the patch scale of 0.4 m when the gas friction is S₁ = 0.5. It is noted that one attenuation coefficient peak appears at low frequencies corresponding to the mesoscopic loss mechanism. The transition behavior is weaker both present in velocity dispersion and attenuation. At the frequency range of 100–1,000 Hz, the attenuation coefficient curves turn to increase because of the global flow Biot mechanism. The viscoelastic equivalent models represent this changes on the whole. Figure 7 shows the phase velocities (a) and loss angles (b) for the Rayleigh waves in case of Fig. 6, except for the gas friction S₁ = 0.1. The low-frequency transition variations in velocity and attenuation are similar to those of the P wave, which illustrates that the low-frequency attenuation effect of the surface wave is dominated by the P wave. The viscoelastic equivalent models also present these changes in velocity and attenuation. It can be concluded that the Zener viscoleastic equivalent model can effectively represent the properties of the Rayleigh wave propagation considering the mesoscopic loss mechanism at low frequencies reliably.

Fig. 6

Fit of both the phase velocity (a) and loss angle (b) curves of the Rayleigh surface waves in a mixture of water and air patchy saturated porous solid by using viscoelastic representation. The patchy saturation models are obtained by White (1975) as corrected by Dutta and Seriff (1979), and Johnson (2001). The material properties are shown in Table 1. S₁ = 0.5 and b = 0.4 m

Fig. 7

Fit of both the phase velocity (a) and loss angle (b) curves of the Rayleigh surface waves in a mixture of water and air patchy saturated porous solid by using viscoelastic representation. The patchy saturation models are obtained by White (1975) as corrected by Dutta and Seriff (1979), and Johnson (2001). S₁ = 0.1 and b = 0.4 m

We extend the study to realistic time-dependent responses of the full wave modes in the White and Johnson patchy saturation models and the corresponding viscoelastic equivalent representations. In these cases, we calculate the vertical displacement responses corresponding to the boundary conditions of Eqs. (20) ~ (22). The point force impulse is selected as a Ricker wavelet written by

$$F(t)=F_{\text{n}}(1-2{\mathit{\pi }}^2{f}_{\text{s}}^2{t}^2)exp({\mathit{\pi }}^2{f}_{\text{s}}^2{t}^2),$$

(23)

where f_s = 35 Hz is the peak frequency. The magnitude F_n = 100 N and time shift 0.1 s are selected. Using a fast Fourier transform algorithm (Zhang et al. 2012), the dynamical response is obtained by multiplication of the spectra of the Green’s functions and the source function (Eq. (23)). Figure 8 shows the results at the time windows for the Rayleigh wave modes at z = 0 and offset r = 200 m. The shear waves are overlapped in the strong Rayleigh waves. The arrival time and the amplitude of the Rayleigh wave characterize the velocity and attenuation in patchy saturated medium, which are consistent with the previous frequency-dependent velocity dispersion and attenuation analyses. The distinguishable differences of the patchy saturated models and the viscoelastic equivalent models cannot be observed. The viscoelastic equivalent representations for the patchy saturated models provide the enough reliable approximations for the Rayleigh wave forms.

Fig. 8

Fit of the wave forms of the Rayleigh waves in a mixture of water and air patchy saturated porous solid by using viscoelastic representation at S₁ = 0.5 (a) and S₁ = 0.1 (b). The patchy saturation models are obtained by White (1975) as corrected by Dutta and Seriff (1979), and Johnson (2001). The material properties are shown in Table 1. The source impulse is 35 Hz Ricker wavelet, and offset is 200 m

5 Discussion and conclusions

In this work, we have studied the equivalent representations of the Rayleigh surface wave that propagates along a free surface between vacuum and a patchy saturated solid, which can effectively represent the attenuation of the surface waves avoiding large computer consumptions. The numerical results show several interesting features when a mixture of water and air inhomogeneously distributed in a mesoscopic scale, i.e., the patchy saturated gas pocket is larger compared to the dimensions of the wavelength of the fast P wave and shear wave and less compared to the dimensions of the each pore size. At low frequencies, the mesoscopic patchy saturation effect is a major loss mechanism for the P wave propagation. The viscoelastic representation can well fit the velocity and attenuation dispersion curves, especially at low frequencies when we choose the relaxed modulus to describe the viscoelastic elements. The Rayleigh wave shows observable transition behaviors in the seismic frequency band, which suggests to be dominated by the P wave affected by the patchy saturation. The velocity dispersion and attenuation become remarkable as the gas friction is decreased. And the surface wave mode shows a well-defined maximum in the attenuation at low frequencies and a turn to increase at high frequencies. These frequency-dependent properties of the surface waves are well represent by the viscoelastic element fitting. The efficiency of the viscoleastic representation is observed in the time-dependent dynamical responses of the surface wave, as well.

It can be concluded that the viscoelastic equivalent representations for the surface wave propagations in patchy saturated poroelastic media provide the well reliable approximations at low frequencies. The viscoelastic representation is more convenient, since the viscoelastic modeling reduces the resource consumption of poroelasticity and supplies an applicable approach to build a simple model for further applications.

Appendices

Appendix 1: analytical expressions of K(ω)

1.1 White’s model

White (1975) firstly developed the theory for a gas-filled sphere of radius a inside a water-filled sphere of outer radius b (a < b). Assuming that the dry-rock, grain modulus, and permeability of the different regions are the same, the dynamic bulk modulus as a function of frequency is given by (Dutta and Seriff 1979)

$$K(\mathit{\omega })=\frac{K_{\text{GH}}}{1-W{K}_{\text{GH}}},$$

(24)

where K_GH is given by Eq. (10). Moreover,

$$W=\frac{3ia\mathit{\kappa }(R_1-R_2)}{b^3\mathit{\omega }({\mathit{\eta }}_1{Z}_1-{\mathit{\eta }}_2{Z}_2)}\left(\frac{M_1}{K_{\text{G1}}}-\frac{M_2}{K_{\text{G2}}}\right),$$

(25)

$$R_1=\frac{(K_{\text{G1}}-K_{\text{m}})(3K_{\text{G2}}+4\mathit{\mu })}{K_{\text{G2}}(3K_{\text{G2}}+4\mathit{\mu })+4{\mathit{\mu }}_{\text{m}}(K_{\text{G1}}-K_{\text{G2}})S_1},$$

(26)

$$R_2=\frac{(K_{\text{G2}}-K_{\text{m}})(3K_{\text{G1}}+4\mathit{\mu })}{K_{\text{G2}}(3K_{\text{G2}}+4\mathit{\mu })+4{\mathit{\mu }}_{\text{m}}(K_{\text{G1}}-K_{\text{G2}})S_1},$$

(27)

$$Z_1=\frac{1-exp(-2{\mathit{\gamma }}_{1}a)}{({\mathit{\gamma }}_{1}a-1)+({\mathit{\gamma }}_{1}a+1)exp(-2{\mathit{\gamma }}_{1}a)},$$

(28)

$$Z_2=\frac{({\mathit{\gamma }}_{2}b+1)+({\mathit{\gamma }}_{2}b-1)exp[2{\mathit{\gamma }}_2(b-a)]}{({\mathit{\gamma }}_{2}b+1)({\mathit{\gamma }}_{2}a-1)-({\mathit{\gamma }}_{2}b-1)({\mathit{\gamma }}_{2}a+1)exp[2{\mathit{\gamma }}_2(b-a)]},$$

(29)

where K_Gi are the Gassmann modulus given by Eq. (7) and the parameter M _i is given by Eq. (6).

And

$${\mathit{\gamma }}_i={\left(\frac{i\mathit{\omega }{\mathit{\eta }}_i}{\mathit{\kappa }{K}_{\text{A}i}}\right)}^{1/2},i=1,2,$$

(30)

$$K_{\text{A}i}=\frac{K_{\text{m}}}{K_{\text{G}i}}{M}_i,i=1,2,$$

(31)

where η _i are the fluid viscosities.

1.2 Johnson’s model

Johnson developed a simple model for the dynamic bulk modulus K(ω), which describes the crossover between the two frequency limits. This theory has two parameters, S/V, the ratio of the surface area of a patch to the sample volume, and T, which is related to the geometry of the patches. These two parameters appear in the expressions for the high- and low-frequency limits. The low- and high-frequency limits can be given by

$$\underset{\mathit{\omega }\to 0}{lim}K(\mathit{\omega })=K_{\text{GW}}[1-i\mathit{\omega }T+\cdots ],$$

(32)

$$\underset{\mathit{\omega }\to \infty }{lim}K(\mathit{\omega })=K_{\text{GH}}[1-G{(i\mathit{\omega })}^{-1/2}+\cdots ]$$

(33)

where

$$G={\left|\frac{\mathit{\Delta }p}{p_{\text{e}}}\right|}^2\frac{S}{V}\sqrt{D_{\text{e}}},$$

(34)

The effective diffusivity D_e is

$$D_{\text{e}}={\left[\frac{\mathit{\kappa }{K}_{\text{GH}}}{{\mathit{\eta }}_1\sqrt{D_1}+{\mathit{\eta }}_2\sqrt{D_2}}\right]}^2,$$

(35)

where D _i is given by Eq. (9). The resulting discontinuity in pore pressure Δp, relative to the applied external stress p_e, is constant along the interface between the two regions and is given by

$$\frac{\mathit{\Delta }p}{p_{\text{e}}}=\frac{\mathit{\alpha }{M}_2\left(K_{\text{G1}}+\frac{4}{3}\mathit{\mu }\right)-\mathit{\alpha }{M}_1\left(K_{\text{G2}}+\frac{4}{3}\mathit{\mu }\right)}{S_1{K}_{\text{G1}}\left(K_{\text{G2}}+\frac{4}{3}\mathit{\mu }\right)+S_2{K}_{\text{G2}}\left(K_{\text{G1}}+\frac{4}{3}\mathit{\mu }\right)},$$

(36)

K(ω) has the expression

$$K(\mathit{\omega })=K_{\text{GH}}-\frac{K_{\text{GH}}-K_{\text{GW}}}{1-\mathit{\xi }+\mathit{\xi }\sqrt{1-i\mathit{\omega }\mathit{\sigma }/{\mathit{\xi }}^2}},$$

(37)

where

$$\mathit{\sigma }={\left[\frac{K_{\text{GH}}-K_{\text{GW}}}{K_{\text{GH}}G}\right]}^2,$$

(38)

and

$$\mathit{\xi }=\frac{K_{\text{GH}}-K_{\text{GW}}}{2K_{\text{GW}}}\left(\frac{\mathit{\sigma }}{T}\right),$$

(39)

The parameters ξ and σ have a precise physical significance, named as a shape parameter and frequency scale, respectively. The two Johnson’s parameters can be defined for this spherical patch as

$$S/V=3\frac{a^2}{b^3},$$

(40)

$$T=\frac{K_{\text{GW}}{\mathit{\varphi }}^2}{30\mathit{\kappa }{M}_2^3{\mathit{\varphi }}^3}\{[3{\mathit{\eta }}_2{g}_2^2+5({\mathit{\eta }}_1-{\mathit{\eta }}_2)g_1{g}_2-3{\mathit{\eta }}_1{g}_1^2]a^5-15{\mathit{\eta }}_2{g}_2(g_2-g_1)a^3{b}^3+5g_2[3{\mathit{\eta }}_2{g}_2-(2{\mathit{\eta }}_1-{\mathit{\eta }}_2)g_1]a^2{b}^3-3{\mathit{\eta }}_2{g}_2^2{b}^5\},$$

(41)

where

$$g_i=\frac{(1-K_{\text{m}}/K_{\text{s}})(1/K_{\text{W}}-1/K_{\text{f}i})}{1-K_{\text{m}}/K_{\text{s}}-\mathit{\varphi }{K}_{\text{m}}/K_{\text{s}}+\mathit{\varphi }{K}_{\text{m}}/K_{\text{W}}},i=1,2,$$

(42)

It is easily understood that the physical significance of the parameter T represents the diffusion time for equilibrating stress in the porous skeleton over the size of the fluid patch (Johnson 2001). For the corrected model considering fluid surface tension influence (Tserkovnyak and Johnson 2003), the static limit bulk modulus of K_GW is replaced by the revised modulus K₀, in which the effective fluid modulus (see Eq. (59) in Tserkovnyak and Johnson 2003) replaces the original one. The nondimensional parameter $\mathit{\Delta }=\frac{s\mathit{\varphi }{\mathit{\sigma }}_0}{\mathit{\kappa }{K}_{\text{f}}{k}_i}$ can identify the effects of this mechanism (Markov and Levin 2007; Markov 2009a). Where s is a shape factor which becomes much less than one as the pore structure becomes irregular (Nagy and Blaho 1994). σ₀ is the surface tension between the two fluids, and k _i is the wavenumber. These effects may appear in the patchy saturation for a definite rock sample, which can be extended straightforwardly and will be studied in the future.

Appendix 2: displacement potential decomposition for the surface wave

For Rayleigh type surface wave that propagates along a free surface between vacuum and a patchy saturated porous half space, we consider a cylindrical coordinate to describe the axisymmetric problem. We introduce the transformation $f^{(\text{m})}(k)={\int }_0^{+\infty }rf(r)J_{\text{m}}(kr)\text{d}r$ to obtain the Hankel spectrum of the function f(r) · J_m denotes the m-th order first kind Bessel function and k denotes the horizontal wave number.

In the frequency domain, by the potential decompositions for the displacement in the effective viscoelastic medium (z > 0) can be written in Hankel spectra by

$$u_z^{(0)}=-A_0{\mathit{\gamma }}_p{\text{e}}^{-{\mathit{\gamma }}_{p}z}+C_{0}k{\text{e}}^{-{\mathit{\gamma }}_{s}z},$$

(43)

$$u_{\text{r}}^{(1)}=-A_{0}k{\text{e}}^{-{\mathit{\gamma }}_{p}z}+C_0{\mathit{\gamma }}_s{\text{e}}^{-{\mathit{\gamma }}_{s}z},$$

(44)

where A₀ and C₀ are the undetermined coefficients. The wavenumbers in the z direction are related to the horizontal wavenumber k through the following relation

$${\mathit{\gamma }}_i=\sqrt{k^2-k_i^2},i=p,s,$$

(45)

where k _i is the complex wavenumber of the corresponding viscoelastic equivalent body modes, which can written by $k_i^{}=\sqrt{\mathit{\omega }/{\tilde{v}}_i}$.

The boundary conditions are rewritten in cylinder coordinate in the frequency and wavenumber domain by

$${\mathit{\tau }}_{zz}^{(0)}=-\frac{F(\mathit{\omega })}{2\mathit{\pi }},$$

(46)

$${\mathit{\tau }}_{rz}^{(1)}=0,$$

(47)

where F(ω) is the spectrum of F(t) in Eq. (20). By the constitution relation, the stress and pressure components can be written by

$${\mathit{\tau }}_{zz}^{(0)}=(-\mathit{\lambda }{k}_p^2+2\mathit{\mu }{\mathit{\gamma }}_p^2)A_0{\text{e}}^{-{\mathit{\gamma }}_{p}z}-2\mathit{\mu }k{\mathit{\gamma }}_s{C}_0{\text{e}}^{-{\mathit{\gamma }}_{s}z},$$

(48)

$${\mathit{\tau }}_{rz}^{(1)}=\mathit{\mu }[2k{\mathit{\gamma }}_p{A}_0{\text{e}}^{-{\mathit{\gamma }}_{p}z}+({\mathit{\gamma }}_s^2+k^2)C_0{\text{e}}^{-{\mathit{\gamma }}_{s}z}],$$

(49)

Substitution Eqs (46) and (47) into Eqs. (48) and (49), the coefficients can be obtained as

$$A_0=-\frac{(k^2+{\mathit{\gamma }}_s^2)}{\mathit{\mu }R}\frac{1}{2\mathit{\pi }}F(\mathit{\omega }),$$

(50)

$$C_0=-\frac{2k{\mathit{\gamma }}_p^{}}{\mathit{\mu }R}\frac{1}{2\mathit{\pi }}F(\mathit{\omega }),$$

(51)

where

$$R={(k^2+{\mathit{\gamma }}_s^2)}^2-4k^2{\mathit{\gamma }}_p{\mathit{\gamma }}_s,$$

(52)

where R is the Rayleigh wave dominator, in which the horizontal wavenumber k of the surface wave mode can be defined in R = 0 at different frequencies. The time domain Green’s functions of the full wave modes along the free surface are obtained by inverse transformation of the Hankel spectra (Eqs. (43), (44), (48) and (49)) to the space–time domain. The numerical implement of the integrations can be put forward in the fast coverage algorithm (e.g., van Dalen et al. 2010; Zhang et al. 2012).