Velocity dispersion and attenuation in residual saturated porous media strongly depend on the properties of the discontinuous fluid phase. We present a multiscale model based on an oscillator rheology for the fluid blobs/clusters, taking into account their distinct dynamic properties on the pore scale.


Propagation of seismic waves in partially saturated porous media depends on various material properties, including saturation, porosity, elastic properties of the skeleton, viscous properties of the pore fluids, and, additionally, capillary pressure and effective permeability. If the wetting fluid is in a discontinuous state (i.e., residual saturated configuration), phase velocities and frequency-dependent attenuation additionally depend on microscopical (pore-scale) properties such as droplet and/or ganglia size. To model wave propagation in residual saturated porous media, we developed a three-phase model based on an enriched continuum mixture theory capturing the strong coupling between the micro- and the macroscale. The three-phase model considers a continuous and a discontinuous part. The continuous part exhibits similar behavior as the poroelastic model introduced by Biot. The discontinuous part describes the movement of blobs/clusters of the wetting fluid and is based on an oscillator rheology. In comparison with other three-phase models, the presented one accounts for the heterogeneity of the discontinuous fluid clusters by use of their dynamic properties, i.e., their statistically distributed inertia, eigenfrequency, and damping effects. This heterogeneous and discontinuous distribution of the wetting fluid in the form of single blobs or fluid clusters is represented by a model-embedded distribution function of the cluster sizes. We define a dimensionless parameter that determines if the overall motion of the residual fluid is dominated by oscillations (underdamped, resonance) or not (overdamped). Our results show that the residual fluid has a significant impact on the velocity dispersion and attenuation no matter if it oscillates or not. For long wavelengths our model coincides with the Biot–Gassmann equations. We show under which conditions and how the classical biphasic models can be used to approximate the dynamic behavior of residual saturated porous media.

To understand and characterize the dynamical behavior of partially saturated porous media, such as soils, rocks, or organic materials, the partial properties of the solid skeleton, of the inherent pore fluids and, in addition, the main physical coupling effects between these interacting constituents have to be taken into account. Since the seminal work of Biot (1956a,b) and Frenkel (1944) in the middle of the last century, many theoretical and numerical studies about wave propagation phenomena in fully saturated porous media have been published (e.g., Bourbié et al., 1987; Stoll, 1989; Carcione, 2007, and references therein). To describe the macroscopic behavior of seismic waves propagating through partially saturated media (i.e., porous media saturated with a wetting and a nonwetting pore fluid), the biphasic Biot-type approach was extended to take into account quasistatic (Santos et al., 1990; Tuncay and Corapcioglu, 1996; Wei and Muraleetharan, 2002; Carcione et al., 2004; Lo and Sposito, 2005; Albers, 2009) and dynamic capillary effects (Lu and Hanyga, 2005). In contrast to Biot's poroelastic model, three longitudinal waves and one shear wave are predicted in such three-phase approaches. These three-phase models are able to describe phase velocities and attenuation as a function of saturation of the wetting fluid, sw. In addition to the frequency-dependence of attenuation and phase velocities, the results of these models also depend on capillary pressure, pc. Available experimental results of lab-scale wetting–dewetting experiments (soil-water characteristic curves) can be adopted straightforwardly. Standard pc(sw) functions of van Genuchten (1980) or Brooks and Corey (1964) are embedded in such models. In low-frequency (<5 kHz) laboratory experiments, Murphy (1982) systematically studied the dependence of the fast P-wave and the S-wave on saturation of the wetting fluid, sw, for water-saturated Massilon sandstones.

Common to all these continuum models, describing propagating waves on a macro scale, is the fact that they do not cover the case of residual saturation of the porous medium. In this case blobs or clusters of wetting fluid are disconnected or are connected via very thin water films. As such water films can be considered negligible in terms of their influence on wave propagation, the water distribution for the purpose of this study can be approximated as disconnected water blobs and clusters of different sizes and shapes, which occupy only a single or several pores. The size of the individual fluid clusters can be widely distributed, for example in a polydisperse granular medium. The properties (effective phase velocities and attenuation) strongly depend on this distribution of fluid cluster size. Propagating waves at low seismic frequencies (f < 100 Hz) are not affected by the discontinuous wetting fluid because the wavelength is much larger than the characteristic clusters’ sizes. On the other hand, a wave propagating through the continuous solid skeleton or through the continuous nonwetting fluid at higher frequencies is able to excite the blobs of wetting fluid through an exchange of momentum.

This phenomenon could be adopted in future applications. The nondestructive characterization of partially saturated soils in the vadose zone is still a challenge. The present mathematical model could be the basis of an inverse parameter identification tool, characterizing the soil with respect to its amount of saturation and determining the effective sizes of fluid clusters. Especially the information about the characteristic length scale of the wetting fluid in the residual state (i.e., specific surface areas) is one major attribute of the proposed model. Exciting and mobilizing discontinuous wetting fluids through an imposed flow of the nonwetting fluid is not yet used in technical applications. Nevertheless, first studies and experiments have been performed, for example as enhanced oil recovery (Beresnev and Johnson, 1994) or enhanced groundwater remediation (Reddi and Challa, 1994; Reddi and Wu, 1996). The acoustically enhanced ganglia mobilization and dissolution of dense non-aqueous phase liquids was successfully studied experimentally in the low-frequency range (10–225 Hz) in a monolayer of glass beads (Chrysikopoulos and Vogler, 2006) and in a sand core (Roberts et al., 2001). Whether such methods will turn out to be successful or not, the prospective knowledge about wave propagation in partially saturated media will help to characterize such materials. Theoretical investigations of the excitation and movement of blobs and clusters of wetting fluid in porous media were performed for idealized pore geometries on the pore scale (Hilpert et al., 2000; Beresnev, 2006; Hilpert, 2007; Holzner et al., 2009). Depending on the geometrical and surface properties of the pore space, Hilpert et al. (2000) distinguished between pinned and sliding motion of the contact lines between a single fluid cluster and the solid skeleton. When externally excited, a frequency-dependent resonance behavior of the blobs with a maximum response at the resonance frequency of the blob was observed. This resonance behavior will be one important element of the presented model.

On the macroscale, Frehner et al. (2009, 2010) investigated isotropic elastic solids, representing a residual saturated porous medium as a single-phase material. They proposed a spectral modification of the elastic wave around the resonance frequency of the trapped residual wetting fluid. The derived rheology consists of an elastic part for the solid and an oscillator part describing the discontinuous blobs or clusters of wetting fluid. Recently, Steeb et al. (2010) extended this approach to multiphase materials. For this, the continuous elastic part of the model was replaced by a poroelastic medium. Consequently, and well-known from biphasic poroelastic materials (Biot, 1956a), the continuous part of the proposed model is able to predict two longitudinal wave modes and one shear wave mode. Furthermore, the low- and high-frequency limits in this model coincide with Biot's theory. A distinct deviation from Biot's theory in phase velocities and attenuation is predicted around the resonance frequency of the blobs and clusters of wetting fluid.

The presented theory is based on the basic ideas of Frehner et al. (2009, 2010) and Steeb et al. (2010). It develops the physical framework from basic principles including the mathematical homogenization process of the new physical properties. The distinct fluid patches are modeled as harmonic oscillators and can account for three different types of attenuation: oscillations at their eigenfrequency, viscous effects inside the fluid body, and sliding motion relative to the solid. There are two key features of the model that support further investigation of residual saturated porous media. First, the distinction between the inherent attenuation mechanisms (energy loss due to oscillations or viscous attenuation) helps characterizing the dynamic behavior of the system. Second, the physical phenomena on the microscale are linked to the macroscale parameters by a concise upscaling technique. This allows for a predication of the dynamic values on the macroscale by knowledge about the microscopic morphology.

The presented model enhances existing macroscale wave propagation models fundamentally. Local or squirt flow models (e.g., Dvorkin and Nur, 1993) deal with an additional length on the microscale and introduce attenuation via an additional mode of motion. The physical relation to the fluid patches, like their geometry, surface tension, or viscosity, is not taken explicitly into account. The presented model does not enhance the intrinsic motion of all phases; it introduces new dynamical mechanisms for the wetting phase based on the distribution of different disconnected fluid clusters. Local flow is included within the possible damped oscillation modes of the discontinuous phases and contributes to the attenuation. In contrast, so-called “patchy-saturation” models consider a partially saturated porous medium with two fluid phases, which can both flow (e.g., White et al., 1975; Dutta and Odé, 1979a,b; Norris, 1993; Toms et al., 2006; Quintal et al., 2011). These patchy-saturation models predict attenuation and dispersion of seismic waves in the low-frequency range (<100 Hz), which is caused by wave-induced fluid flow (Pride et al., 2004). This wave-induced fluid flow is the result of pore pressure differences in the two fluid phases, which have different compressibilities (e.g., gas and water). Extended continuum models of higher order, like in Metrikine and Askes (2002) enhance the kinematical behavior of a single phase. For example, such classes of models take into account the motion of the microstructure of the material independently from the motion on the macroscale. Nevertheless, the aim of extended continua is different from the current work, because the latter includes new physical mechanisms. Again, the special case of distributed and disconnected fluid clusters with their own physical properties is not included in kinematically extended continua.

The goal of the present study is twofold. On one hand, we rigorously derive a model for residual saturated porous media based on first physical principles. Thus, we sketch the main physical assumptions of our approach in detail and we outline the capabilities and limits of the model. On the other hand, we show how the model can be linked to pore-scale properties of the medium. As the model is two-scale by nature, we outline the upscaling procedure of a heterogeneous distribution of blobs/clusters of wetting fluid (different in size/shape) on the level of a representative elementary volume (REV). Furthermore, we demonstrate how geometrical and material properties on the pore-scale are linked to the proposed macroscopical approach.

The model's development starts at the microstructure and focuses on the behavior of oscillating fluid clusters. An upscaling process conserves important information for wave propagation on a larger scale and connects the disconnected oscillator to the continuous phases. Results of exemplary residual saturated rocks show the capabilities and differences to classical models, followed by a discussion.

Microscale Modeling

According to their dynamical behavior, the model embeds discontinuous fluid clusters into the continuous phases as distinct oscillators. Many times, and especially in mixture theory, a link between two scales is applied for modeling purposes. In such a way Biot related the characteristic frequency to the microscopic flow through a duct (Biot, 1956b), and the partial and the effective density are linked via upscaling to the porosity.

Here, as well as in the examples before, an upscaling process connects the physics on the microscale with the dynamic parameters of the system on the macroscale. In this context, the microscale is defined by the length scale of the pore space, where fluid-flow is described by the Navier–Stokes equations. Additionally, the mechanical behavior of the porous solid skeleton is determined by common elasticity models. The principle of scale separation is applied as the macroscale is much bigger than the grain/pore size and acts on the assumption of a spatially homogenized system.

Foremost, a clear knowledge of the physical mechanisms on the microscale is essential for the development of a phenomenological continuum model. Only then, the derived macroscale model is able to capture the effective behavior of a unit cell (REV) including the microscale dynamics of clusters of wetting fluid. Surface tension, density, viscosity, and geometry of the wetting fluid are the most important factors on the microscale, resulting in (damped) oscillations of the fluid clusters. The following section gives a focused overview of the important physical equations describing an oscillating, microscopic fluid patch. The link to the macroscale is given by the upscaling process and an illustration of the new parameters.

The Physical Microscopic System

The microscopic system of interest consists of a discontinuous, incompressible wetting fluid in the form of bridges, blobs, ganglia, or clusters embedded in a solid matrix and surrounded by a compressible nonwetting fluid. In the following, we use the terms blob or cluster in a general meaning to account for the occurrence of possible geometries of the discontinuous fluid. We will discuss in detail how the distinct dynamic properties of the discontinuous fluid will be taken into account. In the present contribution, we assume isothermal conditions and describe the movement of the blob relative to the solid only.

This section describes the physical phenomena to motivate the dynamic behavior of an oscillating fluid cluster, to prove assumptions made during the reduction to an oscillator model, and to support physical understanding.

Surface Tension and Contact Angle

The dominating effect on the kinematics of a single blob of wetting fluid is surface tension. It results from the fact that the attraction forces of molecules are uniform in all directions inside the fluid, but different at the interface with a second material. A certain amount of energy is therefore needed to move a molecule from the interior to the boundary of the fluid cluster and release the connection to the former neighborhood. The required energy for surface increase is termed surface energy or, more commonly, surface tension, σ, and is introduced as a tangential force per unit length. Naturally, it depends on the properties of both materials on either side of the interface.

Furthermore, σ depends only on temperature as long as the characteristic length of the fluid cluster is greater than the Tolman length, δTolman (Tolman, 1949). Because the following investigation only focuses on fluid blobs with diameters d >> δTolman and because the system is treated isothermally, the surface tension is a constant material parameter; i.e., σ = constant.

The surface tension has two important consequences on the physics (de Gennes et al., 2004). First, the Young–Laplace equation relates the surface tension to a pressure difference, Δp, across the interface between two materials. Because the surface tries to minimize its energy, work is done against the arising pressure difference and yields the local balance equation, 
where H is the mean curvature and can be expressed in terms of the divergence of the normal vector, n, on the surface. Second, Young's equation determines the contact angle, Θ, at equilibrium at the contact line of three materials (i.e., two fluids and one solid) when all three surface tensions balance by 
In addition, the contact angle can be used instead of those surface tensions that are affected by the solid. This simplifies matters because measurements and experimental data are mostly provided in the form of surface tension between liquids and liquid/gas and contact angles between fluids and solids. Moreover, the movement of the solid is insignificantly influenced by surface tension. Therefore, the solid can be reduced to its influence on Θ.

Balance Equations and Boundary Conditions

Inside the wetting fluid, the continuity equation holds as well as the Navier–Stokes equations. The interface between the wetting and the nonwetting fluid is influenced by the pressure of the nonwetting fluid, pn, and additionally by the pressure difference due to surface tension, σ. The interface between the solid and the wetting fluid can exhibit a no-slip condition for pinned oscillations, or it can be a wetted wall leading to a shear stress depending on the relative velocity between the two media. For the interface between the nonwetting fluid and the solid, a no-slip condition is usually assumed while the influence of the solid on the blob of wetting fluid is reduced to the contact angle, Θeq, at equilibrium (Eq. [2]).

The Microscopic Oscillator

According to the physical mechanisms summarized above, the restoring force in an oscillating fluid blob is caused by surface tension. A movement of the blob away from the equilibrium position of minimal energy results in a change of the blob geometry. The increased surface or the changed curvature and the resulting change in pressure difference, Δp, results in a restoring force. Furthermore, the movement of the blob is attenuated by internal viscous flow (viscosity μwR) and, possibly, by a sliding motion relative to the slipping walls.

Considering small displacements around equilibrium and including inertia terms yields a single damped oscillator equation for the barycentric movement of a wetting blob i. If the displacement of the center of gravity relative to the solid's wall, uwi, is small, the momentum exchange with the solid balances the inertia term as follows: 
In Eq. [3], miωi2 accounts for the elastic restoring force, which is expressed in terms of the mass, mi, and the eigenfrequency, ωi, of the fluid blob and the parameter ci represents viscous damping effects. The latter can also be replaced by a mass specific damping parameter di with ci = midi. The physical behavior of the oscillator and its influence on the physical system depends on the dimensionless damping ratio ζi = ci/(2miωi). Overdamped oscillators (ζi > 1) are dominated by viscosity and show an exponential decay without oscillations. The fastest return to equilibrium occurs for critically damped oscillators (ζi = 1). The underdamped case (ζi < 1) still shows the influence of the stiffness including damped oscillations around equilibrium. This distinction and its influence on propagating waves in a residual saturated porous medium will be further discussed in section 4.

Upscaling Process

Because we are mainly interested in larger-scale geophysical applications, we consider scales much larger than the pore scale. Therefore, the physical properties of disconnected blobs of wetting fluid in an appropriate REV have to be upscaled and integrated into a macroscale continuum model.

If we sum up the microscale physics of the last sections and compare it to a classical harmonic oscillator, we observe that the dynamic behavior can be divided into three main parts:

  • pure pinned oscillation around equilibrium

  • viscous attenuation inside of the cluster

  • (viscous) attenuated sliding motion with respect to the solid wall

As the fluid patches transfer a wave only by connection to the continuous phases, they can store and dissipate energy by the three mentioned processes and influence dispersive wave propagation for those reasons. All blobs are characterized with respect to this behavior and determined by their inertia, damping (subsuming point two and three) and eigenfrequency. The main idea of this homogenization is now to account for the spatial distribution and the corresponding variety of blob types. Therefore, we conserve information about the microscopic oscillations of the blobs and the heterogeneous blob distribution in the considered REV (Steeb, 2010; Steeb et al., 2010). The entirety of oscillators is subjected to an upscaling process, but with the eigenfrequency as a differentiating factor. The eigenfrequency is used because it will turn out to be the important parameter for low saturated situations, when the system behaves differently from continuous, that is, viscous damping dominated ones. This means that a distinct averaged oscillator model is created for every possible eigenfrequency ωk (Fig. 1). The mathematical homogenization procedure is performed for the momentum balance equation. Based on the microscopical results above, three major assumptions are used in this process:

  • 1. The fluid blobs exchange momentum only with the solid skeleton (momentum exchangeforumla). This assumption results from the fact that the vadose zone is composed of the porous solid skeleton saturated with a gas (air) and a liquid. Due to the much lower density and much lower bulk modulus of the continuous gas phase, the exchange of momentum between the liquid and the gas phase can be neglected.

  • 2. The effective density of the different fluid clusters is assumed to be constant. Because the elastic behavior of the fluid clusters is mainly determined by the restoring force of the surface tension, the compressibility of the wetting phase is negligible and included in the general oscillation modes. Furthermore, changes of temperature and pressure are small with respect to their influence on the density.

  • 3. Scale separation is assumed because the wavelength λ of the incident wave is much larger than the characteristic size L of the discontinuous blobs; i.e., λ >> L. As an example, we refer to a typical reservoir rock (Berea sandstone). The lower limit of a fast P-wave is around cP1 ≈ 2600 m s−1. Thus, for a maximum ultrasound frequency of 100 kHz, the wavelength is λ ≈ 0.164 m. This frequency used, for example, in laboratory investigations, is still much larger than the characteristic length scale of the rock microstructure, i.e., the pore size, and, therefore larger than the characteristic length scale L of the blobs in the pores. Thus, a continuum approach is valid. On the other hand, the frequency of 100 kHz includes the expected resonance effects of the discontinuous phase (cf. Hilpert et al., 2000). For materials with significantly different properties, scale separation has to be proved separately.

It has to be noted that the homogenization procedure can also be executed without the first two assumptions, leading to more general results. This work focuses on the mentioned case, because it appears in many systems of the vadose zone and still allows a detailed characterization without unnecessary complication. A volume average is calculated for fluid blobs with the specific eigenfrequency ωk in the REV. The detailed homogenization process is outlined in Appendix A1 and leads to the nonequilibrium momentum exchange 
which is balanced by the inertia forces.

Note that uwkus is the averaged macroscale displacement of fluid blobs with eigenfrequency ωk relative to the solid skeleton. Besides the standard averaged quantities of displacement, effective density, ρwR, and volume fraction, nw, three new macroscopic values are introduced: ak is the volumetric portion of the discontinuous wetting fluid with a certain eigenfrequency ωk and ck is the corresponding averaged damping parameter due to inner viscous flow and sliding motion. (Appendix A1). They allow considering the required microstructural information of the discontinuous wetting fluid. This classification and differentiation of the discontinuous fluid patches becomes important for the characterization of systems with very low saturation.

As we will observe later in detail, a distribution of fluid blobs of different sizes can be described by a set of Eq. [8], which extends the classical poroelastic approach.

Macroscale Modeling

After upscaling the mechanical properties of the disconnected fluid clusters, the continuous phases will be introduced and subsequently linked to the discontinuous fluid phase. As a result, we obtain the final model.

General Definitions and Assumptions

On the macroscale, a REV of residual saturated porous material is composed of three constituents denoted by the indices α ∈ {s, n, w}: a continuous solid skeleton, s, a continuous nonwetting fluid phase, n, and a discontinuous wetting fluid phase, w (Fig. 1). Note that the discontinuous fluid phase can also be nonwetting in the present model. However, we focus here on a porous skeleton saturated with a discontinuous liquid and a continuous gas phase where the liquid constituent is usually the wetting phase.

The macroscale modeling approach is based on the thermodynamically consistent theory of porous media (Bowen, 1980; Ehlers and Bluhm, 2002), which extends the classical mixture theory (Truesdell, 1957) by the concept of volume fractions. The volume fractions, nα, of the constituents φα are the quotients of the volume occupied by the constituent in the REV and the volume of the REV, dv; i.e., nα = dvα/dv. According to the scale separation principle, the REV is much bigger than the grain size, such that the microscale can be assumed to be homogenized with respect to it. Thus, we introduce the partial density, ρα = dmα/dv, and the effective or true density of the constituents, ραR = dmα/dvα Based on the volume fractions, we introduce the porosity, ϕ = 1 − ns, and the amount of fluid saturation, sβ = nβ/ϕ, with β ∈ {n, w}. We assume that the wetting fluid is discontinuously distributed in the REV. Therefore, considering a small perturbation around the equilibrium state, the wetting fluid cannot flow through the pore space, but the single fluid blobs can deform and are therefore able to exchange momentum with the solid skeleton. In the following, we develop a phenomenological model, capturing the main effects of propagating seismic waves through such a REV, based on the microscopical considerations of the section Microscale Modeling.

Continuous Phases

The mathematical framework of the continuous phases is given by the balances of momentum; see Steeb et al. (2010) for a detailed description. We consider only small perturbations around an equilibrium state and restrict ourselves to the linearized form of the balance relations. Thus, we do not have to distinguish between material time derivatives and local, i.e., partial time derivatives, forumla. Neglecting body forces, the local form of the momentum balance equations is given for all constituents, φα as, 
Besides the partial stress tensor, Tα, and the inertia forces, forumla, the momentum of a constituent is balanced by the momentum exchange term,forumla. Note that forumla denotes here the exchange of momentum and not the fluid pressure, which is contained in the total partial stress tensor Tα. Next, we discuss this interaction term,forumla, between the constituents.

Field Equations

The macroscale model consists of the homogenized continuous phases (solid and nonwetting fluid) and the upscaled oscillators, describing the discontinuous blobs of wetting fluid. The total dynamic behavior of the residual saturated mixture in the REV is affected by the momentum balance of both continuous phases, φs and φn, respectively, and the oscillator-type rheology of the discontinuous blobs.

The field equations of the continuous phases and the different oscillators can be combined by the momentum exchange given in Eq. [4] (Appendix A2). Fluid–fluid interaction is neglected due to the big difference in density while a Darcy-type viscous coupling is assumed between the solid and the continuous fluid phase expressed by the viscous permeability function b0 = μnRϕ/ks. Here, we introduced the intrinsic permeability ks and the effective viscosity of the nonwetting phase μnR. A transformation to wave number-frequency space with a standard harmonic ansatz yields the corresponding eigenvalue problems for the longitudinal mode (index P; P-wave) and the transversal mode (index S; S-wave) of the displacement vector u 
where kj are the wave numbers and (see Appendix A3) 
The eigenvalue problem consists of the matrix A including complex densities forumla (see Eq. [8,9]). The defined complex densities account for inertia (∝ ω2) and viscous effects (∝ ω), respectively. Matrix Bj includes the material parameters (A, N, S, R) that are the same as in Biot's theory (Biot, 1956a; Bourbié et al., 1987; Stoll, 1989) or the theory of porous media (Steeb, 2010) and can be tied back to the intrinsic elastic properties of the continuous phases (Appendix A4).

The momentum balances of the oscillators relate their own displacement quantities to the displacements of the solid phase. As a result, the final system of Eq. [6] yields only two P-wave wave modes (one S-wave mode), which are related to the solid and the continuous nonwetting phase and is very similar to poroelastic models.

The oscillating fluid clusters contribute to the system by interaction with the solid phase. Mathematically they show an impact on the complex densities forumla, introduced as 
This determines the final model and Eq. [8] shows how the oscillator rheology is able to enhance the system by a variety of oscillating fluid blobs with different mass (αkρw), eigenfrequency (ωk) and damping behavior (ck) with only one addition inside of the governing equations.
In addition, approximating the discrete distribution of oscillating fluid blobs with a probability density function (PDF), α(ω), the sums in Eq. [8] approach integrals. forumla has to be defined between two positive boundary frequencies, 0 < dstart < dend, and has the dimension of [s]. The eigenvalue problem of the continuous case using a PDF is identical to the discrete one, except that the complex density, 
consists of integrals instead of sums.

It should be emphasized that the governing eigenvalue problem (Eq. [6] with Eq. [8] and [9]) is mathematically similar to the classical poroelastic problem where two longitudinal waves and one transversal wave are predicted (Biot, 1956a). The discontinuous wetting fluid influences the phase velocity and attenuation of the P1- and the S-wave, but it does not predict any further wave modes. The mathematical reduced eigenvalue problem of the three-phase system results from the constitutive coupling of the displacement of the solid and the blobs. Physically, no wave can propagate through the disconnected fluid phase. The blobs may affect the attenuation or wave speed but need to transmit the wave via their connection to the continuous phase(s).

Solution of the Dispersion Relations

The dispersion relations can be determined from Eq. [6] for a system with a discrete distribution of blobs of wetting fluid or its approximation by a PDF. Solving the eigenvalue problem with respect to the squared wave numbers kj2(ω) that are denoted as kj2(ω)=: ξj(ω), j ∈ {P, S}, leads to the three mentioned wave modes (P1, P2, and S-waves).

Longitudinal Mode



Transversal Mode



The quality factors of the wave modes, QP1,P2, (P-waves) and QS (S-wave), expressing the amount of intrinsic attenuation and the phase velocities, cP1,P2 and cS, are defined as




Damping and Characteristic Frequencies

To study the influence of damping of the model, we focus on a system with one single oscillator with damping c1 and eigenfrequency ω1, because this case is the most understandable. Figures 2 and 3 show respectively phase velocity dispersion and frequency-dependent attenuation of the P1-wave for this model. Only one dimensionless parameter, D, affects the characteristics of the P1- and the S-wave:



The characteristic damping parameter D represents the relationship between the two damping mechanisms of the discontinuous fluid due to viscous damping with nwc1 and to the oscillation at ω1 with the partial density ρw, respectively. For large values of D (log10(D) > 0.5), the phase velocity of the fast longitudinal wave, cP1, is qualitatively comparable to that of a classical poroelastic model (Fig. 2). However, the presented model exhibits a different damping behavior (Fig. 3) and the characteristic frequency is shifted to higher values compared to a classical poroelastic model (Fig. 2 and 3). If the oscillation dominates over damping (log10(D) < −0.1), the dispersion relation contains a distinct peak at resonance frequency. In general, the characteristic frequency of the oscillating wetting fluid can be described as



bosc can be interpreted as a more general parameter of momentum interaction compared to poroelastic models, where it only accounts for viscous attenuation and does not account for oscillations. Over all, the qualitative behavior of the oscillator in the system is different if it is underdamped (D << 1) or overdamped (D >> 1). Still, a fundamental difference to a discrete oscillator appears (Eq. [3]): the damping of the oscillating fluid blobs results in a higher critical frequency compared to a poroelastic model. In contrast, for a classical, damped harmonic oscillator (stiffness c, mass m, and undamped eigenfrequency ω) the characteristic frequency decreases with increasing damping; i.e.,forumla. This fundamental difference is due to the viscous coupling of the oscillating fluid blobs to the elastic porous solid. Although a real physical system consists of a manifold of different blobs, which can be described by a PDF, this simple example illustrates the influential physical phenomena and their consequences.

We have not taken into account the effect of tortuosity. Furthermore, the critical frequency ωc,osc (Eq. [14]) depends on the density ratio ρws, which can reach higher values for highly porous materials with a large amount of residual saturation.

Comparison with Biot's Theory

The presented model describes the influence of the continuous nonwetting fluid and of the discontinuous oscillating blobs/clusters of wetting fluid. Therefore, two different physical phenomena appear in the dispersion relation. On one hand, there is friction due to the movement of the viscous nonwetting fluid described by the permeability function, b0. On the other hand, the fluid blobs gather energy when they oscillate at their eigenfrequency and dissipate energy due to their own intrinsic viscous attenuation. In general, these second effects are not captured in Biot's two-phase poroelastic model. However, under certain conditions, a poroelastic model can be used to approximate the presented residual saturation model. These conditions are satisfied if the energy, stored by the blobs’ oscillations, becomes small compared to the one of general viscous damping (D >> 1). Such an approximation helps to better understand and simplify the more complex model. Therefore, we compare a one-oscillator model (ω1, c1) with two alternative poroelastic models I and II. Their mathematical description can be found in Appendix A3.

Model I treats the continuous fluid as the second phase of a two-phase poroelastic system. If one is interested in attenuation effects and phase velocities in a frequency range where this continuous, nonwetting fluid (gas) is dominating, it could be chosen as an appropriate approximation (Fig. 4 and 5). Note that this model is exact for the P2-mode, as this mode is not influenced by the wetting blobs.

In contrast to the first model, model II accounts for the influence of the wetting fluid (water) in case of D >> 1 (with boscnwc1), meaning that the oscillations are damped by the viscosity of the discontinuous wetting fluid. The second phase of this poroelastic system consists of the wetting fluid, composing the fluid blobs. If the nonwetting fluid is a gas, its influence on the first P- and the S-wave is small, and model II approximates well the residual saturation model for large values of D (Fig. 4 and 5).

We have introduced both models to separate the physical attenuation and dispersion effects of the continuous nonwetting fluid from the wetting blobs. Furthermore, simpler and well-understood poroelastic models can highlight the characteristic properties of the complex three-phase model.

Illustrative Example: Influence of a Heterogeneous Distribution of Blobs

So far, we only considered oscillating fluid blobs with a single eigenfrequency. However, the presented residual saturation model is capable of handling a distribution of eigenfrequencies, represented by a PDF (Eq. [9]). A PDF is a better representation of the natural distribution of blob/cluster sizes in a REV. The following example investigates a possible water cluster distribution in Berea sandstone in combination with air. We assume an artificial lognormal distribution of eigenfrequencies around a central frequency ω0 = 10e5·2π [1/s] :



Equation [15] is depicted in Fig. 6a for different values of s. Because the damping of the oscillating blobs and clusters of wetting fluid is not only a function of the fluid viscosity, but also of the friction at the pore walls, we assume that the damping parameter is a function of the size of the fluid blobs/clusters and, therefore, of the eigenfrequency. Here, we assume a linear relationship,



Equation [16] is depicted in Fig. 6a. The phase velocity, cP1, and the intrinsic attenuation, 1/QP1, of the P1-wave are depicted in Fig. 6b and 6c, respectively, for different values of s. For a very narrow distribution of eigenfrequencies, the phase velocity dispersion and the attenuation strongly resemble the single-oscillator case with a small value of parameter D (Fig. 2 and 3). For a wider distribution of eigenfrequencies, the curves resemble the single-oscillator case with a larger value of D. With increasing eigenfrequency distribution width, the attenuation distribution widens and the attenuation peak decreases and shifts to higher frequencies. This can be explained by the fact that with a wider distribution of eigenfrequencies, more high-frequency oscillators with a larger attenuation contribute to the total attenuation.


We presented a mathematical model that determines the interaction of the motion of a discontinuous wetting fluid (blob) on the microscale (i.e., pore scale) with the propagation of waves in a biphasic porous medium on the macroscale (i.e., scale of the wavelength). On the microscale, we described the motion of fluid blobs by the equation of a damped oscillator. There are two end-member scenarios for the motion of fluid blobs: first, if the damping ratio is smaller than one (underdamped) the fluid blobs oscillate, and large oscillation amplitudes can occur around the resonance frequency. Second, if the damping ratio of the harmonic oscillator is larger than one (overdamped), no oscillations and no resonance effects occur. In this case, the residual fluid blobs act as a viscous damper that absorbs energy from the passing wave. Our results show that for both, the overdamped and the underdamped situation, the presence of fluid blobs can have a significant impact on the wave velocities. The magnitude of the resonance frequency of fluid blobs has been investigated with analytical and numerical methods (Hilpert et al., 2000; Beresnev, 2006; Hilpert, 2007), and it was shown that the resonance frequency strongly depends, among other parameters, on the assumed size of the pore space in which the fluid blob oscillates (e.g., the pore itself, a connected pore channel or a microfracture). Hilpert et al. (2000) presented resonance frequencies for typical parameters of oil and pore size between 50 Hz and a few 1000 Hz. Furthermore, Hilpert (2007) showed that results from Lattice–Boltzmann simulations agree well with the analytically predicted resonance effects and frequency.

In the final analysis, the model predicts wave velocities and attenuation for different types of distributed fluid patches. First of all, this relationship allows characterizing the system. Conclusions about the distribution of the discontinuous phase are possible, especially if most parameters of the porous material are known. Often, for example, in the vadose zone, this is the case, and only the varying amount of saturation and especially the effective characteristic length of fluid clusters are unknown. To be more precise, for cases of residual saturations, the current model can be used as the theoretical basis for the inverse analysis of the determination of average blob sizes. Therefore, the identification process makes use of the previously discussed peaks of the attenuation spectrum. This may support the prediction of immobilized fluid by information about characteristic damping frequencies of remaining liquid blobs besides common form-locking domains (Fagerlund et al., 2007). The model can also be conducive to laboratory experiments, if properties like interfacial areas should be linked to measurable quantities of propagating waves. The presented theory allows one to combine and interpret measurements on both scales and to develop adequate constitutive equations.

There is still an open question of whether resonance effects can be actively used in, for example, enhanced oil recovery (e.g., Beresnev and Johnson, 1994) or groundwater remediation (e.g., Reddi and Challa, 1994; Reddi and Wu, 1996). For both of the above-mentioned applications it is necessary to mobilize residual organic liquids (usually oil) that are trapped in porous media (e.g., Li et al., 2005). For example, Pride et al. (2008) argued that resonance frequencies of residual oil blobs are higher than the seismic frequency band of interest for enhanced oil recovery, and they showed results of Lattice–Boltzmann simulations suggesting that oscillatory wave stimulation of reservoirs can improve oil recovery without resonance effects. It should be noted that in all cases, a foresighted treatment of the stimulated area has to be assured with respect to possible damages of the solid framework due to the stimulating waves. The presented model may help to characterize or describe the behavior of such areas with low saturation on the basis experiments, for example, ultrasound tests on the lab scale.

In general, for natural situations the pore space and the distribution of fluid blobs are heterogeneous. In our model we assume a certain resonance frequency of the residual fluid blobs, based on data like Hilpert et al. (2000). In the future, we will apply numerical pore-scale simulations based on the finite element method to calculate the resonance frequency of heterogeneous porous media in a similar way as has been done using Lattice–Boltzmann simulations (Hilpert, 2007). This technique can be applied if the pore space and the distribution of fluid blobs in a natural porous rock can be measured with three-dimensional imaging tools, such as, for example, synchrotron-based X-ray tomographic microscopy (Li et al., 2001). These results for the microscale resonance will then be combined with an accompanying eigenvalue analysis and will serve as input data for the macroscale model of wave propagation. Our theory predicts that the wave velocity in a poroelastic layer with a discontinuous wetting phase (e.g., residual liquid blobs) differs from the velocity for a continuous wetting phase. This prediction could be tested in either laboratory experiments or using so-called time-lapse seismics, where the wave velocity of reservoirs is recorded over a certain time interval during production. For example, our theory predicts that a “full” reservoir exhibits a different wave velocity than a reservoir containing residual discontinuous liquid blobs. Because the wave velocity of the reservoir changes, the reflection coefficient between the reservoir and the surrounding rock changes as well. The change in reflection coefficient can be predicted by our model and could be measured for a natural reservoir. Such changes in reflection coefficient due to change in fluid content have been observed for hydrocarbon reservoirs (e.g., Korneev et al., 2004) and also in reservoirs used for CO2 storage in the framework of CO2 sequestration (e.g., Rubino et al., 2011). The changes in reflection coefficient are usually explained with wave-induced fluid flow effects (e.g., Rubino et al., 2011). However, the theory we presented here might better describe the mechanism causing dispersion, attenuation and the related change in reflection coefficient for reservoirs with a residual and discontinuous wetting fluid phase.


We presented a mathematical model that couples the microscale (i.e., pore scale) motion of a discontinuous wetting fluid (blob) with the macroscale (i.e., scale of the wavelength) propagation of seismic waves in a porous medium. The macroscale balance equations are based on the theory of porous media, and the microscale motion is upscaled using a homogenization procedure and coupled to the macroscale equations using consistent momentum exchange terms. We formulated the resulting macroscale equations in the same mathematical structure as the classical Biot equations for wave propagation in a poroelastic medium, and we included all terms describing the microscale effects in the complex densities. Three probability density functions can describe and account for the variety of different fluid clusters with respect their distributed eigenfrequencies, mass, and attenuation within one additional term. This particular form of the equations provides a transparent way of describing the wave propagation phenomena investigated.

The presence of a discontinuous wetting fluid on the pore scale has a distinct influence on the phase velocities of the P1- and the S-wave and also on the corresponding attenuation in the porous medium. However, the discontinuous wetting fluid does not cause any additional wave modes as in continuous three-phase models, where waves are also transmitted by a second connected fluid phase with its own compressibility and inertia.

The motion of the discontinuous wetting fluid is described by the equation of a damped oscillator. The fluid can oscillate (underdamped) or not (overdamped) depending on the damping ratio of the oscillator. We derived a dimensionless parameter, D, similar to the damping ratio which (i) determines if wave propagation is affected by microscale fluid oscillations or viscous damping and (ii) controls the characteristic frequency at which dispersion occurs, due to microscale processes. For fluid oscillations, the dispersion curve exhibits strong variations around the resonance frequency whereas for fluid viscous damping, the dispersion curve has a similar shape as the dispersion in the classical Biot poroelastic theory. Our results indicate that seismic wave propagation in a porous medium can be affected by the presence of a discontinuous residual wetting fluid, whether the residual fluid oscillates or not.

Appendix A1

Homogenization of the Properties on the Microscale

The starting point is the exchange of momentum of one single (harmonic) microscopic oscillator



with its mass, mi, its eigenfrequency, ωi, the mass specific damping, di, and the displacement of its center of gravity relative to the motion of the surrounding solid material in the REV, forumla Moreover, a standard harmonic ansatz has been used. Summation over all oscillators with the same eigenfrequency ωk in one REV and dividing by the reference volume of the REV, dv, leads to the macroscopic momentum exchange for all non-wetting fluid clusters with the eigenfrequency ωk:



We have assumed that the blobs of wetting fluid have identical effective density, i.e., ρkwR = ρwR. This allows interchanging derivatives and a density weighted averaging of kinematical quantities. Furthermore, we assume that each point of the solid walls moves like the solid constituent in average



Additionally, we introduce the averaged displacement for all oscillators with the eigenfrequency ωk as



Since the influence of the pore-scale is determined by the different eigenfrequencies of the blobs of wetting fluid, the averaged mass specific damping parameter, dk, consists of the mass weighted average over all blobs of wetting fluid of the same eigenfrequency



The less these microscopic damping factors, di, diverge from each other, the less they deviate from the average. The used damping parameter is not mass specific; i.e., ck := ρkwRdk .

Appendix A2

Solution and Simplification of the Corresponding Eigenvalue Problem

The field equations of the continuous phases are extended by z distinct oscillators



Each oscillator equation is related to its characteristic viscous damping behavior and to its eigenfrequency. Besides a Darcy-type viscous momentum exchange, forumla, between the nonwetting fluid and the solid skeleton, expressed by the viscous permeability function b0 = μnRϕ0/ks (Biot, 1956a), the additional amount of momentum exchange,forumla, appears between the solid skeleton and the wetting fluid overall as well as the specific oscillators



These relationships for the interaction of the phases can be combined with the field equations to remove the unknowns forumla. Via usual splitting techniques and a standard harmonic ansatz (Steeb, 2010), the eigenvalue problem for z distinct oscillators can be derived from this balance of momentum as



where the symmetric matrices and the displacements are given by






Since the matrices forumla have only zero coefficients for the third and all further rows, one can write



by executing the matrix-vector multiplication for every oscillator with eigenfrequency ωk. Equation [A27] represents the equilibrium of momentum, which includes inertia terms, as well as friction and elastic forces. The latter depend on the microscopic eigenfrequency, ωk. If ak is non-zero, then



Replacing the unknown displacements of the oscillators by Eq. [A28] yields the reduced eigenvalue problem, Eq. [6] and [7], with only two degrees of freedom.

Appendix A3

Poroelastic Models for Comparison

Depending on whether the critical frequency of the fluid blobs, ωc,osc, is smaller or larger than the critical transition frequency of the continuous phase, ωn, we have to distinguish the mentioned cases, because the inherent fluid phases are attenuated in a different order.

Model I





Model II





Appendix A4

Biot Parameters

In the case of a homogeneous porous medium with connected pores, the parameters N, A, S, R, and P are generalized elastic material parameters. These quantities can be related to the intrinsic properties of the nonwetting fluid and the solid constituent (Biot and Willis, 1957). Therefore, we introduce the (dry) shear modulus of the skeleton, G, the (dry) bulk modulus of the solid skeleton, K, the bulk modulus of the grains composing the solid skeleton, Ks, and the bulk modulus of the non-wetting fluid, Kn. Furthermore, we define the effective porosity ϕR := ϕ + Kn/Ks(1 − ϕ − K/Ks) and the elastic parameters



Patrick S. Kurzeja gratefully acknowledges financial support in the form of a PhD scholarship from the Studienstiftung des deutschen Volkes (German National Academic Foundation). Stefan M. Schmalholz was supported by the University of Lausanne.

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Open access article.