Effective Fluid Bulk Modulus in the Partially Saturated Rock and the Amplitude Dispersion Effects

Frequency dispersion is a well‐known effect in geophysics, which means that waves of different wavelengths propagate at different velocities. Amplitude dispersion is a less‐known effect, which means that waves of different amplitudes propagate at different velocities. Herewith, we consider the alteration of the interfacial energy during wave‐induced two‐phase fluid flow in a partially saturated rock and demonstrate that this leads to a nonlinear amplitude dispersion effect. When the wave amplitude is small, seismic waves cause bending of the interface menisci between immiscible fluids at the pore scale. However, when the wave amplitude is sufficiently large, the interface menisci will slip at the pore scale, causing attenuation of the elastic energy by the contact line friction mechanism. At the zero frequency limit, all viscous dissipation models predict zero attenuation of the elastic wave energy, while this approach predicts a nonzero attenuation due to a static contact angle hysteresis effect. Herein, we extend the Gassmann's theory with three extra terms, which can be obtained from standard laboratory tests: pore‐size distribution and interfacial tension between immiscible fluids and rock wettability (advancing and receding contact angles). We derive closed‐form analytical expressions predicting the effective fluid modulus in partially saturated rock, which falls between Voigt and Reuss averages. Next, we demonstrate that the nonlinear amplitude dispersion effect leads to energy transfer between different frequencies. This may explain the low‐frequency microtremor anomalies, frequently observed above hydrocarbon reservoirs, when the low‐frequency energy of ocean waves (0.1–1 Hz) is converted to higher frequencies (2–6 Hz) by partially saturated reservoirs.


Introduction
To extract fluid types or saturations from seismic or borehole sonic data, we need a procedure to model fluid effects on rock velocity and attenuation. Numerous techniques have been developed. However, Gassmann's (1951) equations are by far the most widely used relations to calculate seismic velocity changes caused by different fluid saturation in reservoirs (Berryman, 1999;Mavko et al., 2009;Müller et al., 2010;Pride et al., 2004). The importance of accurate saturation prediction grows as seismic data are increasingly used during reservoir monitoring and during monitoring of carbon dioxide storage (e.g., Bergmann & Chadwick, 2015;Dupuy et al., 2017;Furre et al., 2017;Grude et al., 2014). Gassmann's relations were originally derived to describe the change in rock modulus from one pure saturation to another-from dry to fully brine saturated, from fully brine saturated to fully oil saturated, and soforth. Hydrocarbon reservoirs and CO2 storage formations are almost never fully saturated either with oil or gas phases. It always contains some partial fractions of oil (or/and gas) and water at the pore scale. Domenico (1976) suggested that mixed gas-oil-brine saturations can also be modeled with Gassmann's relations, if the mixture of phases is replaced by a single fluid phase with effective (equivalent) bulk modulus. The effective bulk modulus of immiscible pore fluids is not well defined and remains an active topic of the research (e.g., Amalokwu et al., 2017;Behzadi et al., 2011;Brie et al., 1995;Cadoret, 1993;Domenico, 1976;Dvorkin & Wollner, 2017;Monachesi et al., 2020;Müller et al., 2008;Papageorgiou et al., 2016;Papageorgiou et al., 2018;Rozhko & Bauer, 2019;Rozhko, 2019;Wollner & Dvorkin, 2016;Wollner & Dvorkin, 2018a, 2018b. According to Mavko and Mukerji (1998), it is possible to define the upper and lower bounds for the effective bulk modulus of the pore fluid mixture. The upper Voigt bound is calculated by arithmetic (isostrain) volume average of fluid properties: where S we is the wetting phase saturation and K we and K nw are bulk moduli for the wetting and nonwetting fluid phases, respectively.
The lower Reuss (also known as Wood) bound is calculated by harmonic (isostress) volume average of fluid properties: These theoretical upper and lower bounds give a significant uncertainty in saturation prediction. Brie et al. (1995) suggested the empirical mixing law, based on experimental testing of sandstone, partially saturated with liquid and gas: where e is an empirical constant with a typical value equal to 3. Brie's empirical correlation recovers the Voigt upper bound when e = 1, but it is not consistent with Reuss lower bound for sufficiently large values of e and K nw . Mavko et al. (2009) suggested a practical solution to assume that the effective fluid modulus falls roughly between Reuss bound and Brie average with e = 3. Wollner and Dvorkin (2016, 2018a, 2018b argued that effective fluid bulk modulus can be approximated by a linear combination of Voigt and Reuss bounds, depending on the elastic moduli of the rock and the pore fluids. They suggested that the relative contrast between the gas and water is larger in soft rock as compared to rock with a stiffer frame modulus. Another important conclusion was recently made by Papageorgiou et al. (2016), who demonstrated that the empirical Brie parameter is related to the ratio of fluid pressure increments in the wetting and nonwetting fluid phases and may not be necessarily related to the patch size and frequency, as it is usually considered (e.g., Mavko et al., 2009;Mavko & Mukerji, 1998). This ratio is not a material parameter, which can be easily measured; thus, it is challenging to apply Papageorgiou et al. (2016) results for saturation prediction. Rozhko (2019) and Rozhko and Bauer (2019) predicted analytically the bulk modulus of the rock with an isolated partially saturated microcrack of given length and aspect ratio. Practically speaking, it is difficult to apply the model of Rozhko (2019) and Rozhko and Bauer (2019) because it depends on unknown microcrack geometry, which is extremely difficult to determine from laboratory data. Contrary, all input parameters in Gassmann's theory can be derived from the lab, because this theory does not depend on the pore geometry. In this paper we extend both approaches of Gassmann (1951) and Rozhko (2019Rozhko ( , 2020 by coupling of the interfacial energy to Gassmann's theory. Interfacial energy effects on seismic wave velocity dispersion have been reported in many publications (Knight et al., 2010;Knight & Nolen-Hoeksema, 1990; Moerig et al., 1996;Murphy, 1984;Murphy et al., 1986;Papageorgiou et al., 2016;Rozhko & Bauer, 2019;Vo-Thanh, 1995;Waite et al., 1997;Wang et al., 2015;Wyllie et al., 1958); see Rozhko (2019) for detailed literature review. The coupling is done by three extra parameters, which can be obtained from standard laboratory tests: pore-size distribution and interfacial tension between immiscible fluids and rock wettability (advancing and receding contact angles). Therefore, our new model shares the same input parameters as Gassmann's theory plus these three additional parameters. The boundary condition at the interface between two immiscible fluid phases is described by a static contact angle hysteresis phenomenon. We derive the closed-form analytical expressions predicting effective fluid bulk modulus in the partially saturated rock, which falls between Voigt and Reuss averages. It follows from the model that the effective bulk modulus of pore fluid depends on the wave amplitude. It decreases with increase of the wave amplitude and approaching the Reuss lower bound. When the wave amplitude is small, seismic waves cause bending of the interface menisci between immiscible fluids at the pore scale. When the wave amplitude is sufficiently large, the interface menisci will slip at the pore scale, causing attenuation of the elastic energy by the contact line friction mechanism (Rozhko, 2019(Rozhko, , 2020Rozhko & Bauer, 2019). In other words, when the wave amplitude is small, the interface menisci restrict the relative motion of fluids during wave-induced two-phase fluid flow. If the wave amplitude is sufficiently large, the relative motion of fluids would be enabled. At the zero frequency limit, all viscous dissipation models predict zero attenuation of the elastic wave energy (Mavko et al., 2009;Müller et al., 2010;Pride et al., 2004), while this approach predicts a nonzero attenuation due to a static contact angle hysteresis effect (Rozhko, 2019(Rozhko, , 2020Rozhko & Bauer, 2019). The nonlinear amplitude dispersion effect, predicted by our model, is well known in theoretical physics (Karpman, 1975;Landau & Lifshitz, 1987;Whitham, 2011). It brings a whole range of phenomena, which are not found in linear poroelastic materials; particularly, it may cause the energy transfer between different frequencies (Karpman, 1975;Landau & Lifshitz, 1987;Whitham, 2011). This may explain the low-frequency microtremor anomalies, frequently observed above hydrocarbon reservoirs (e.g.: Dangel et al., 2003;Goertz et al., 2012;Lambert et al., 2009;Riahi et al., 2012;Steiner et al., 2008;Saenger et al., 2009;Witten & Artman, 2011). Kazantsev (2018) and previous authors concluded that the physical mechanisms behind the low-frequency microtremor anomalies are unclear in the literature. Several authors argued that these anomalies could be related to the anthropogenic noise, geological structure, surface waves, human errors during interpretation, and so forth (e.g., Ali et al., 2010a;Hanssen & Bussat, 2008;Martini et al., 2013;Woith et al., 2014). In this study we demonstrate that these anomalies could be explained by the transfer of the low-frequency energy of ocean waves (0.1-1 Hz) to higher frequencies (2-6 Hz) by a partially saturated reservoir due to a nonlinear amplitude dispersion effect caused by a static contact angle hysteresis.
The article is constructed as follows. In section two we derive Gassmann's equation for a single-phase saturation using a Betti-Rayleigh reciprocity theorem. Next, in section three we explain the contact angle hysteresis effect and provide a short introduction to the capillarity of the rock. Next, in section 4 we couple the interfacial effects, outlined in section 3 to the Gassmann's theory, outlined in section 2. In section 5 we provide numerical results for the effective fluid bulk modulus, Skempton's B coefficients (for the wetting and nonwetting fluid phases) and attenuation factor. Next, in section 6 we demonstrate that a static contact angle hysteresis effect leads to the nonlinear energy transfer between different frequencies by calculating the spectrum of the bulk-volume strain rate. In section 7 we conclude the paper, while section 8 contains the nomenclature list. Auxiliary materials contain a Matlab script with implemented equations derived in this paper and a Maple script, which shows derivation of equations.

Gassmann's Theory
In this section we derive Gassmann's equations for fully saturated rock with one fluid phase. Walsh (1965) showed that an expression for the effective compressibility of a porous elastic material can be derived in terms of the compressibility of the solid material (C 0 ) and the rate of change of pore-volume δϑ P with external stress increment δσ. Consider a porous material with bulk-volume of ϑ B containing pores of pore-volume ϑ P . The expression for the effective compressibility can be written in the following form (Walsh, 1965): where equation (4) follows from the Betti-Rayleigh reciprocity theorem (e.g., Schmeling, 1985). The term on the left-hand side of equation (4) represents the effective compressibility of the rock, that is, while the term is the bulk-volume strain. Equation (4) is the general equation, which is independent on the boundary conditions for the fluid phase (drained or dry vs. undrained) and independent on the saturations (full or partial). This equation relates the change of bulk volume to the change of pore volume, while the change of pore volume will of course depend on the boundary conditions for fluid phases and on the saturation.
Small perturbations of the pore volume and bulk volume are calculated using equations and Here p fl is the fluid pressure, while δp fl is the wave-induced perturbation (increment) of the fluid pressure.
To derive Gassmann's equation, let us consider first the dry rock or, equivalently, drained boundary conditions. In this case, the wave-induced perturbation of the fluid pressure will be equal to zero, that is, δp fl = 0; consequently, equations (4) and (5) can be written in the form and The term, which appears on the right-hand side of equation (10), is related to the pore-volume compressibility of dry rock (or drained boundary conditions), which is defined as follows: Using the definition of porosity we can find the relationship between the bulk compressibility and pore-volume compressibility of dry rock: Equations (9), (10), (11), and (13) are well known and can be found in Zimmerman (1990).
Next, we consider undrained boundary conditions, when the pore fluid cannot flow in or out from the rock sample during wave-induced stress perturbation δσ. Consequently, the perturbation of the pore fluid pressure is not equal to zero. Next, we consider that the deformation of the pore volume is controlled by the effective stress (Zimmerman, 1990): where α is the Biot's poroelastic coefficient, where the typical values for α reported for different sandstones are very close to one (Zimmerman, 1990). Note here that we use a sign convention when compressive stresses and compressive strains are negative, while compressive pore pressure is positive.
Using the effective stress principle, the increment of the pore volume (see equation (8)) can be calculated using the chain rule as follows: Using eq. (11) we can re-write eq. (15) in the form: Next, the relation between the pore-volume increment and the fluid pressure increment during undrained boundary conditions can be calculated as where K fl is the bulk modulus of the pore fluid. Next, by combining equations (4), (8), (16), and (17), we derive the following expression for the bulk compressibility of fluid-saturated rock under undrained boundary conditions: Equation (18) is equivalent to the Gassmann equation, provided that bulk moduli are reciprocal to compressibilities (e.g., Avseth et al., 2010;Mavko et al., 2009).
Equation (18) explains how the undrained bulk modulus of the rock will change if the oil (or gas) saturated rock is replaced by water-saturated or dry rock and vice versa. However, this theory does not explain how to calculate the effective bulk modulus of partially saturated rock, which is the most common case.
To conclude this section, we point out that the classical Gassmann's theory suggests that the shear modulus of porous rock does not depend on the fluid saturation, while published laboratory data suggest that this conclusion is not always correct at low frequencies and for isotropic rock (e.g., Adam et al., 2006;Adam et al., 2009;Adam & Otheim, 2013;Bauer et al., 2016;Fabricius et al., 2010;Khazanehdari & Sothcott, 2003;Mikhaltsevitch et al., 2016). Consideration of the fluid effect on the shear modulus is outside of the scope of this paper. It has been addressed recently by Rozhko (2020).

Rock Capillarity and the Contact Angle Hysteresis Effect
In this section we will introduce some basics of reservoir engineering, which will be used further in our model (section 4). When the reservoir rock is partially saturated with two immiscible fluid phases, the fluid pressure in each fluid phase is not identical. The fluid pressure in the wetting phase is defined as p we , while fluid pressure in the nonwetting phase is defined as p nw . The capillary pressure is defined as (e.g., Barenblatt et al., 1990) At the equilibrium condition the fluid pressure within the pore space is not uniform. Small pore throats with the characteristic size smaller than a given capillary cutoff radius (r c ) are fully saturated with the wetting phase, while large pore throats with a size larger than a given cutoff radius are fully saturated with the nonwetting phase. The capillary cutoff radius is defined as (e.g., Washburn, 1921) where γ is the interfacial tension between the wetting and nonwetting fluid phases, while θ is the contact angle made by the wetting fluid phase with the surface of solid. Here we do not consider the adsorption layer of water molecules, which may cover the whole surface of the pore space in the rock. Starov and Velarde (2009) showed that the presence of thin water layer coating the surface of solid may affect contact angle measurements; therefore, the adsorption layer of water molecules can be considered as a part of a solid and not as a part of continuous water phase. Following Brown and Korringa (1975), we also treat bound (immobile) fluids, not connected with the pore as a part of a solid. We do not consider explicitly the surface chemistry effects (and chemical contamination) characterizing the physicochemical interactions between immiscible pore fluids and the pore surface. All physicochemical interactions are implicitly described by the surface tension and advancing and receding contact angles, which are sensitive to chemical composition of pore fluids and wettability of the rock.
The capillary pressure in the reservoir rock is commonly described by the Leverett J-function approach (Leverett, 1941;Murphy, 2013): where κ is the permeability of the rock, while ϕ is the porosity of the rock and J(S we ) is the dimensionless Leverett J-function, which depends on the wetting phase saturation S we . Brooks and Corey (1966) studied systematically the Leverett J-function for different types of oil and gas reservoir rocks and concluded that the Leverett J-function can be well described by the following power law correlation: where j o is the dimensionless constant, related to the capillary entry pressure, and λ is the pore size distribution index. For narrow distributions, λ is greater than 2; for wide distributions, λ is less than 2. The typical values for λ reported for different types of oil and gas reservoir rocks are in the range 1 < λ < 4.2 (Brooks & Corey, 1966).
Equation (21) for the capillary pressure can be simplified further, if we consider a Kozeny-Carman relation between permeability of the rock and porosity of the rock (e.g., Mavko et al., 2009): where κ 0 is the proportionality and unity factor [m 2 ]. In this case, the expression for the capillary pressure can be written as Next, to understand better the condition when the change of the wave-induced fluid-saturation is taking place at the pore scale, let us consider the contact angle hysteresis effect. Figures 1a shows schematically a liquid drop surrounded by air on the surface of solid. If we start to inject slowly more liquid (Figure 1a), the contact angle θ 1 will increase (δθ 1 > 0), while the contact line will remain pinned; that is, the radius of wetted area r 1 will remain the same (δr 1 = 0). When the contact angle reaches its critical value (θ 1 = θ a ), the contact line will slip into advancing direction. If we continue to inject very slowly (in a quasi-static manner), the contact line will move into advancing direction (δr 1 > 0) while the contact angle will not change (θ 1 = θ a ), that is, δθ 1 = 0. If we increase the injection rate, the contact line advancing velocity will also increase, in this case the advancing contact angle will increase δθ 1 > 0, because it depends on the contact line motion velocity.
If we start to withdraw liquid from the liquid drop (Figure 1b), the contact angle θ 2 will decrease (δθ 2 < 0), while the contact line will remain pinned; that is, the radius of wetted area r 2 will remain the same (δr 2 = 0). When the contact angle reaches its critical value (θ 2 = θ r ), the contact line will slip into receding direction. If we continue to withdraw very slowly the liquid, the contact line will move into receding direction (δr 2 < 0) while the contact angle will not change θ 2 = θ r , that is, δθ 2 = 0. The change of contact angles within a certain range (θ r ≤ θ ≤ θ a ) is called a contact angle hysteresis, which depends on the contact line velocity (Bormashenko, 2013;de Gennes et al., 2013). At zero velocity the difference between advancing and receding angles is not negligible and can be as large as tens of degrees (Ethington, 1990), and this phenomenon is called a static contact angle hysteresis. If these two angles are equal (θ r = θ a ), it means that the contact line pinning force is equal to zero. Indeed, the injection or withdrawal of fluids ( Figure 1) will lead to the contact line motion without any friction at the contact line location, while the capillary pressure is not equal to zero when θ r = θ a . When the contact line velocity is sufficiently large, the contact angles (θ a and θ r ) depend on the contact line velocity, and this phenomenon is called a dynamic contact angle hysteresis. At seismic frequencies (<200 Hz) when the rate of deformation is very slow, the dynamic effects of the contact angle hysteresis can be neglected, because dynamic corrections are much smaller than a static contact angle hysteresis. At acoustic frequencies (~10 kHz) the dynamic contact angle hysteresis may not be neglected. The Gassmann's theory is a quasi-static theory where the rate of deformation is assumed to be sufficiently slow so that there is enough time for the fluid pressure (in each fluid phase) to equilibrate within the pore space. Consequently, in the low-frequency limit we will consider only a static contact angle hysteresis. The contact angle hysteresis depends also on many other parameters: roughness and mineralogical composition of the substrate surface (i.e., pore surface or grain/crack surface); chemical treatment history of the substrate surface (i.e., wettability alteration history); temperature and fluid pressure; and external electromagnetic fields. These parameters may not be homogeneous within the pore space. We neglect those variations, similarly to Gassmann's theory, assuming a monomineralogical composition of the solid matrix. Furthermore, we neglect the external electromagnetic fields and electrical currents that affect the interfacial tension and contact angles, which are input parameters to our model.
Next, let us discuss the choice of the initial contact angle. In our calculations we assume that the initial contact angle is equal to the most stable angle (θ st ), which corresponds to the minimum of Gibb's free energy (e.g., Drelich, 2019). It is important to note here that any contact angles within the range between θ r and θ a are equilibrium contact angles; however, there is a most stable (i.e., energetically favorable) configuration of the equilibrium contact angle, which can be estimated as (Bormashenko, 2013;Drelich, 2019) Seismic waves may induce residual changes to contact angles, as it will be demonstrated below. Those changes will be metastable, and the contact angle may relax to the most stable configuration over period from hours to days (e.g., Drelich, 2019). Unfortunately, the relaxation time is relatively rare measured in the technical literature (e.g., Drelich, 2019).

Coupling of Interfacial Energy to Gassmann's Theory
Let us consider a periodic (sinusoidal) stress perturbation of amplitude Δσ applied to the external boundary of partially saturated rock (Figure 2a). Typical values of strain amplitudes induced by seismic waves, which can be recorded by seismometers, are very small, typically around 10 −8 to 10 −6 (dimensionless units) (e.g., Aki & Richards, 2002). For the rock with Young's modulus of E~10 GPa, this strain induces stress perturbations around Δσ~10 −2 to 10 −4 Pa. In calculations, we need to consider a much smaller increment of the wave amplitude (δσ ≪ Δσ) of the order δσ~1 Pa (Figure 2b) in order to keep changes of contact angles small.
Next, we are going to derive equations for a partially saturated rock with two immiscible fluids. Equations can be extended toward three-phase saturation; however, it is outside of the scope of this paper. Santos et al. (1990) argued that the deformation of pore volume is controlled by the effective stress (equation (14)), where the fluid pressure is equal to the volumetric averaged effective pressure (p fl ¼ p¯), that is, Thus, similarly to equation (16) the increment of pore volume is calculated as

10.1029/2019JB018693
Journal of Geophysical Research: Solid Earth On the other way, it is common to assume in the effective medium theories that there are no interactions of the stress fields surrounding various pores although those pores can be connected hydraulically (e.g., Mavko et al., 2009;Zimmerman, 1990). It implies that each pore can be mathematically described by an isolated inclusion, where the deformation of inclusion is controlled by the local fluid pressure and far-field stresses (e.g., Mavko et al., 2009;Zimmerman, 1990). Thus, pores saturated by the wetting and nonwetting phases will have different local fluid pressure: p we ¼ p¯− 1− S we ð Þ p cap and p nw ¼ p¯þ S we p cap , respectively. Thus, the deformation of the total pore volume (ϑ P ) and pore volumes occupied by the wetting (ϑ we ) and nonwetting (ϑ nw ) phases will be controlled by different effective stresses. Next, we are interested in the deformation of the total pore volume; therefore, we need to relate deformations of ϑ we and ϑ nw to the deformation of ϑ P by considering linear terms from Taylor series: and Note here the sum of ϑ we and ϑ nw is equal to V P as expected. The stress increment δσ induces perturbation to the equilibrium volumes of the wetting and nonwetting fluid phases, which can be calculated using the product rule as follows: Here the sum of δϑ we and δϑ nw is equal to δV P as expected; therefore, this approach is consistent both with the concept of volumetric-averaged fluid pressure (equation (26)) and effective medium theories with no interactions of the stress fields surrounding various pores (e.g., Mavko et al., 2009;Zimmerman, 1990). It must be noted here that the wetting phase saturation, which appears in the Leverett J-function, does not take into account the deformation of the rock; it is only related to the number of pores with a size smaller than a given cutoff radius r c , saturated by the wetting phase. Therefore, pore volumes calculated with equations (28) and (29) would be different from pore volumes calculated using equations, which does not consider the deformation of pore volumes: and For example, ϑ we calculated with equation (32) is larger than ϑ we calculated with equation (28), as expected, while those differences are negligibly small. Indeed, typical values for the capillary pressure in sandstone reservoirs are below~1 MPa (Brooks & Corey, 1966), while typical values for pore volume compressibilities of reservoir sandstone are around~10 −9 GPa −1 (e.g., Zimmerman, 1990); thus, the term p cap C ϕ~1 0 −3 represents a very minor correction in equations (28) and (29), while this correction is not small in equations (30) and (31) due to equation (27). Therefore, we will neglect those corrections and use equations (32)) and (33) for pore volumes, while keeping those terms during calculations of pore volume increments in equations (30) and (31).
Next, similarly to equation (17), we find relationships between (δϑ we , δϑ nw ) and (δp we , δp nw = δp we +δp cap ) as follows: and where K we and K nw are bulk moduli of the wetting and nonwetting fluid phases, respectively. Equations (34) and (35) correspond to undrained boundary conditions for both fluid phases. Next, let us consider the increment of porosity δϕ, calculated using the product rule and equations (12) and (4) as follows: Finally, we calculate the increment of capillary pressure using equations (24) and (36) as follows: At this point we have derived six linear equations (27), (30), (31), (34), (35), and (37) for seven unknown parameters: δϑ P , δϑ we , δϑ nw , δp we , δp cap , δθ, and δS we , caused by the wave-induced stress increment δσ. In the following two sections, we discuss the effect of boundary condition at the interface between immiscible fluids. Due to this boundary condition, described by a static contact angle hysteresis, the number of unknown parameters is reduced to six, which can be calculated by solving six linear equations.

Contact Line Is Pinned
When the contact line is pinned, we have six linear equations and six unknown parameters: δϑ P , δϑ we , δϑ nw , δp we , δp cap , and δθ, because δS we = 0. Solution to the system of equations (27), (30), (31), (34), (35), and (37) is given by and Using equations (4) and (5), we calculated the effective compressibility of partially saturated rock in the limit when the contact line is pinned: Finally, by equating equations (18) and (44), we derive the expression for the effective bulk modulus of the fluid mixture (fluid bulk modulus) in the limit when the contact line is pinned: It is interesting to note here that in the limit when pores are incompressible C ϕ → 0, equation (45) is equivalent to the Voigt upper bound, calculated by arithmetic volume average in equation (1). In the limit when pores are very soft C ϕ → ∞, equation (45) is equivalent to the Reuss lower bound, calculated by harmonic volume average in equation (2). Hill (1963) argued that both Voigt and Reuss bounds are rather poor when the bulk moduli differ by more than a factor two or so, as it is the case for gas-water mixtures. According to equation (45), those bounds are recovered only in extreme cases, which may explain why those bounds are poor according to Hill. Hill, however, did not consider the impact of interfacial energy in his energy approach during derivation of upper and lower bounds of elastic moduli. As it will be demonstrated in section 5, when the wave amplitude is sufficiently large, and the contact line is slipping, the effective fluid bulk modulus will approach the lower Reuss bound. Consequently, Hill's theoretical bounds are applicable only for elastic systems, because the interfacial energy may cause a deviation from these bounds (see Figure  12a in Rozhko, 2019).
Pore-volume compressibility of the rock, which appears in equation (45), can be obtained from standard laboratory tests (Zimmerman, 1990). It depends on the effective confining stress; see Table 1. Hence, the fluid mixing law also depends on the effective confining stress. Effective fluid bulk modulus increases with increase of the effective confining stress. Thus, equation (45) explains diversity of experimental data, commonly fitted by empirical Brie correlation, which is shown by black curve in Figure 3a for e = 3. Theoretical curves of Figure 3 "look like" a linear combination of Voigt and Reuss bounds, as suggested by Wollner and Dvorkin (2018a), who also argued that the relative contrast between the gas and water is larger in soft rocks as compared to rocks with a stiffer frame. Dvorkin (2016, 2018b) argued that the deviation from the harmonic average is related to the absence of hydraulic communication between domains in rock saturated with different fluids, often called as patchy saturation. In this paper we show that the deviation from the harmonic average is related to the absence of hydraulic communication at the pore scale, due to interface menisci, which restrict the relative motion of fluids during wave-induced two-phase fluid flow. Equation (45) is applicable for the limit, when the contact line is pinned, that is, when the wave amplitude is sufficiently small, that is, smaller than a certain critical value Δσ < Δσ c . Using equation (43), we can calculate the maximum value of the wave amplitude, when the contact line remains pinned, as follows: (46) Figure 3b shows Δσ c , for corresponding curves, shown in Figure 3a. Additional input parameters used in Figure 3 are given in Table 2. In section 5 we will come back to these equations, while in the next section we will consider the case when Δσ > Δσ c , that is, when the contact line is slipping. Note here that using equations (44) and (46), we can also calculate a critical strain-wave amplitude when the contact line remains pinned during seismic wave propagation.

Contact Line Is Moving
When the contact line is slipping, we have six linear equations and six unknown parameters: parameters δϑ P , δϑ we , δϑ nw , δp we , δp cap , and δS we , because δθ = 0. In this case the analytical solutions for the system of equations (27), (30), (31), (34), (35), and (37) are rather cumbersome. Therefore, we will not present these solutions here, which can be found in a Maple script in auxiliary materials. Instead, we will present only the solution for the effective fluid bulk modulus for the stress increment when the contact line is slipping: where and It is interesting to note here that the fluid bulk modulus, calculated with equation (47), is slightly lower than the Reuss lower bound but the difference is negligibly small (<0.2%) for any values of pore volume compressibilities and pores size distributions of sandstone (as one can double-check). Thus, the Reuss lower bound can be used as a good approximation and simplification to equation (47). In our previous publication (Rozhko, 2019(Rozhko, , 2020, we demonstrated that the difference with the Reuss bound may not be small under certain conditions in a material possessing cracks. However, the consideration of the dual porosity is outside of the scope of this paper. In section 4 we derived analytical solutions for small increments of δϑ P , δϑ we , δϑ nw , δp we , δp cap , δθ, and δS we , caused by the stress increment δσ shown in Figure 2b. After each time step of calculation, parameters ϑ P , ϑ we , ϑ nw , p we , p cap , θ, and S we need to be updated by calculated increments. All other parameters, such as porosity, pore volume, and bulk volume should be calculated after each time step using corresponding equations. In the auxiliary materials of this paper, we provided the Matlab script, which does the numerical integration, while the Maple script shows derivations of analytical solutions. Using a Matlab script, it is possible to calculate the effect of different wavelet (sinusoidal, Ricker, etc.) on seismic response.
During derivation of equations, it was assumed that all pores of different sizes have the same pore-volume compressibility C ϕ equal to the averaged pore-matrix compressibility. The real rock may possess a dual porosity with two distinct pore systems such as matrix and cracks, having different pore-volume compressibilities. Herein, we neglected the presence of microcracks. However, Rozhko (2020) demonstrated that the presence of microcracks is responsible for the shear-softening effect (mentioned in section 2). The shear softening means that the shear modulus of the fluid-saturated rock is smaller than the shear modulus of dry rock (Rozhko, 2020). When the presence of microcracks is neglected, we can consider that fluids do not affect the shear modulus of the rock.
Here we do not include the isotropic system of microcracks, considered by Rozhko (2020) because the number of equations and input parameters will be increased. In this paper we aim to keep the minimum number of input parameters to predict the effective fluid bulk modulus in the partially saturated rock and demonstrate the energy transfer between different frequencies due to the amplitude dispersion effect. Similarly, the aim of Rozhko (2020) paper was to demonstrate the shear-softening effect keeping the minimum number of input parameters; that is, only isotropic microcrack porosity was assumed.
The rock physics model developed in this paper is limited to isotropic system, but it can be generalized with little difficulty to anisotropic systems (Berryman, 1999;Brown & Korringa, 1975;Gassmann, 1951). The model can also be generalized to high frequencies, for example, by small modification of Santos et al. (1990) equations. They generalized Biot's theory by considering capillary pressure in the partially saturated rock. Their analysis excluded the contact line pinning force. They related perturbation of the capillary pressure to perturbation of the saturation (only), while according to our model, perturbation of the capillary pressure should be related (1) to perturbation of the saturation, (2) to perturbation of the contact angle, and (3) to perturbation of pore volumes, saturated by the wetting and nonwetting phases. Hence, little modifications to Santos et al. (1990) equations are required to include the contact line pinning force into the generalized poroelasticity theory. Note. Data are taken from Zimmerman (1990, p. 29, Sample 2). Here, the data at σ = − 20 MPa are interpolated.

Numerical Results
Let us consider numerically a periodic (sinusoidal) stress perturbation, of given amplitude Δσ = 4 kPa, applied to the external boundary of the partially gas-water-saturated rock with a given water saturation of S we = 0.8. The frequency of the applied sinusoidal stress perturbation is f = 0.5 Hz.
All other input parameters of the model are given in Table 2. Figure 4a shows a stress perturbation (Δσ) versus time applied to the external boundary of partially saturated rock under undrained boundary conditions for both fluid phases, while Figure 4b shows-Δϵ B -bulk-volume strain versus time. One may notice that the shape of the bulk-volume strain is not sinusoidal. Figure 4b shows-Δp we and Δp nw -changes of fluid pressure in the wetting and nonwetting fluid phases (respectively) versus time. Figure 4d shows-ΔS we and Δ ϑwe ϑP -changes of the wetting phase saturation versus time. It must be noted here that the wetting phase saturation, which appears in the Leverett J-function J(S we ), does not take into account the deformation of the rock (e.g., Coussy, 2004); it is only related to the number of pores with a size, smaller than a given cutoff radius r c , saturated by the wetting phase. So δS we does not represent the actual change of the wetting phase saturation due to deformation of pores; it represents only the change of pore numbers saturated by the wetting phase. Therefore, when the contact line is pinned, the number of pores, saturated by the wetting phase, does not change (δS we = 0), while the actual change of the wetting phase saturation is not equal to zero δ ϑwe ϑP ≠ 0, due to different compressibilities of the wetting and nonwetting fluid phases. The actual change of the wetting phase saturation can be calculated using the chain rule as follows: Hence, using this approach of equation ((48)), we can still calculate the change of the wetting phase saturation due to deformation of the rock. Therefore, Figure 4d shows two curves ΔS we and Δ ϑwe ϑP , representing different properties. The curve ΔS we is related to the change of the relative number of pores saturated by the wetting phase, while the curve Δ ϑwe ϑP is related to the actual change of the wetting phase saturation due to both deformation of pores and the change of the relative number of pores, saturated by the wetting phase. Figure 4e shows-Δθ-changes of the contact angle versus time. Figures 4d and 4e show that S we is increasing when θ = θ a , decreasing when θ = θ r , and does not change when θ r < θ < θ a , as expected, while the actual saturation Δ ϑwe ϑP is changing for any value of the contact angle. Calculations of Figure 4 show that the changes of pore pressure in each fluid phase are very small, that is, within the pressure range that is typical for linear seismic waves (Pride et al., 2004); however, the change of the contact angle is not small and can be sufficient to induce the contact line motion. Pride et al. (2004) argued, however, that as a wave passes, the menisci will bulge and change shape but will not slip away. Our calculations (Figure 4) show that the menisci may slip away, depending on the critical wave amplitude predicted by equation ((46)). Figure 5a shows the bulk-volume strain versus stress curve, that is, Δϵ B versus Δσ. The curve shows a bilinear behavior, which is related to the contact angle hysteresis effect. The local slope of the stress-strain curve is different when the contact line is pinned and when the contact line is slipping. Figure 5a shows the curve starts at initial conditions when Δϵ B = 0 and Δσ = 0 but does not return to its initial condition after one or several oscillation periods (see also Figure 4b). It means that passing seismic wave induces some residual strain. The model also predicts residual changes to all other  Tables 1 and 2. 10.1029/2019JB018693

Journal of Geophysical Research: Solid Earth
ROZHKO parameters, that is, fluid pressure, contact angles, and saturation as it can be seen from Figure 4. The residual changes of the fluid pressure, saturation, and porosity, caused by passage of seismic waves, were observed both in the lab and in the field; see Manga et al. (2012) for literature review. These residual changes typically recover to its initial values over certain time period, which could be explained by recovery of the contact angle to the most stable and energetically favorable configuration at the initial conditions (e.g., Drelich, 2019).
Next we can calculate the effective compressibility of the partially saturated rock by considering the slope of the largest diagonal of the parallelogram, shown in Figure 5a. Note here that, due to nonlinear effects, the definition of compressibility is not unique. Zimmerman (1990) discussed that the compressibility can be calculated by two different methods. The first method conserves the elastic strain, while the second method conserves the elastic energy. The difference between these two methods is small when the hysteresis (nonlinearity) is small. In this paper and in our previous papers (Rozhko, 2019;Rozhko & Bauer, 2019), we used the first method in calculation of bulk compressibility and the second method in calculation of the wave energy dissipation. The elastic energy, dissipated per cycle, is proportional to the area enclosed by the parallelogram, shown in Figure 5a, while the total elastic wave energy is proportional to the area below the largest diagonal of the parallelogram, shown in Figure 5a. Thus, from the geometry of this parallelogram we calculate both the effective compressibility of partially saturated rock (C p. sat = 0.27 GPa −1 ) and the attenuation factor (Q −1 = 0.26). By equating C p. sat , calculated above with equation (18), we calculate the effective bulk modulus of the gas-water mixture (K fl = 0.42 GPa). Using data, shown in Figure 5b, we can calculate Skempton's coefficients for the wetting (B we ) and nonwetting fluid phases (B nw ) (Müller et al., 2010;Pride et al., 2004;Skempton, 1954): Similarly to the compressibility, Skempton's coefficients are not uniquely defined due to nonlinearity effects. In this paper we defined Skempton's coefficients as the ratio at extreme points, that is, by choice of the slope of largest diagonal of the parallelogram, shown in Figure 5b, which gives values of B we = 0.516 and B nw = 0.087. Papageorgiou et al. (2016) argued that the effective compressibility of the fluid mixture should be related to the ratio of fluid pressure perturbation in each fluid phase (Δp we /Δ p nw ), while it was not explained how to predict this ratio, which follows from our model, that is, Δp we Next, we are going to investigate further how K fl , Q −1 , B we , and B nw depend on S we ∈ (0,1) and Δσ ∈ (10 2 ,10 5 ) Pa. Figure 6 shows the dependence of the effective fluid bulk modulus (ordinate) on the water saturation (abscissa) and on the dimensionless wave amplitude Δσ Δσc shown on color scale, where Δσ c versus saturation is shown by cyan curve (with σ = − 20 MPa) in Figure 4b and calculated using equation (46). The cyan curve in Figure 6 is the same as cyan curve in Figure 3a. Voigt, Reuss, and Brie (e = 3) curves are shown by magenta, red, and black colors (respectively) in Figure 6. Figure 6 shows that all data with Δσ ≤ Δσ c collapse on the cyan curve calculated using equation (45), Table 2 Input Parameters of the Model 0.026 · 10 while all data with amplitude Δσ ≥ 12Δσ c collapse on the Reuss lower bound. The data collapse for Δσ ≤ Δσ c is exact, while for Δσ ≥ 12Δσ c is asymptotic. The softening of the effective fluid bulk modulus is predicted for wave amplitudes within the range Δσ c ≤ Δσ ≤ 12Δσ c . Most changes of K fl are taking place within the range Δσ c < Δσ ≲ 5Δσ c . Calculations results (Figure 6) support the suggestion of Mavko et al. (2009) to assume that the effective fluid modulus should fall roughly between the Reuss and Brie (e = 3) averages. Figure 6 shows that this effect is controlled by the wave amplitude. It is interesting to note here that the effect of wave amplitude is opposite to the effect of the confining stress: When the wave amplitude is increasing, the fluid bulk modulus is decreasing, while when the confining stress is increasing (by absolute value), the value of K fl is also increasing, as shown in Figures 6 and 3a, respectively. Figure 7 shows the seismic attenuation factor Q −1 (color scale) as a function of water saturation (abscissa) and the dimensionless wave amplitude Δσ Δσc (ordinate). Figure 7 shows that when the contact line is pinned (Δσ < Δσ c ), seismic attenuation due to the contact line friction is equal to zero, while when the contact line is slipping (Δσ > Δσ c ), seismic attenuation is not equal to zero. The maximum value of the attenuation factor (Q −1 = 0.29) is achieved when Δσ ≈ 1.65Δσ c and S we ≈ 0.91. Figure 7 shows that the low-frequency attenuation can get very large when the contact line is slipping. This may explain the low-frequency attenuation "anomalies" observed on seismic data (Ahmad et al., 2018;Castagna et al., 2003;Ebrom, 2004;Korneev et al., 2004;Sell et al., 2018). Those features, observed on seismic data, are called in literature by "anomalies," because it is challenging to explain it using viscous dissipation models, which predict seismic attenuation to be proportional to the frequency (Q −1~f ) in the low-frequency limit, which gets negligibly small at the zero limit. This model and our previous publications (Rozhko, 2019(Rozhko, , 2020Rozhko & Bauer, 2019) predicted nonzero attenuation at the zero-frequency limit due to a static contact angle hysteresis effect, which was not previously considered in geophysical literature; however, this effect is well known in the physics literature (e.g., Bormashenko, 2013;de Gennes et al., 2013;Drelich, 2019). Figure 8 shows the dependence of Skempton's B we and B nw coefficients (ordinate) on the water saturation S we (abscissa) and the dimensionless wave amplitude Δσ Δσc (color scale). When the wave amplitude is small (Δσ < Δσ c ), Skempton's coefficients do not depend on the water saturation, while when Δσ > Δσ c , the difference between B we and B nw is decreasing with increase of the wave amplitude.
To conclude this section, let us consider two special cases: (1) γ = 0 and (2) θ a = θ r , but γ ≠ 0. In the first case, fluid phases are miscible; hence, the capillary pressure is equal to zero. In the second case fluid phases are immiscible, and the capillary pressure is not equal to zero, but the contact line pinning force is equal to zero; therefore, the contact line can slip without any resistance. In those two cases our model predicts zero attenuation due to contact line friction, and the effective modulus of pore fluid will follow the Reuss lower bound for any wave amplitude. Hence, the deviation of the effective fluid modulus from the Reuss lower bound is related to the contact line pinning force, which is sensitive both on the interfacial tension and upper and lower contact angles.
And finally, all input parameters considered by our model can be obtained from standard laboratory tests: rock compressibilities, porosity, and Biot's poroelastic coefficient from the Rock Mechanical lab tests (e.g., Zimmerman, 1990); pore-size distribution, interfacial tension, and contact angles from the SCAL (Special Core Analysis) lab tests (e.g., McPhee et al., 2015); and fluid compressibilities from the Pressure Volume Temperature lab tests (e.g., Danesh, 1998).

Discussions
Nonlinear amplitude dispersion effects have been studied in theoretical physics for many decades (e.g., Karpman, 1975;Whitham, 2011). The nonlinear amplitude dispersion effect brings a whole range of phenomena that are not present in linear systems (e.g., Whitham, 2011). At the same time, it brings additional computational challenges, as the principle of superposition is not applicable for nonlinear waves (e.g., Karpman, 1975;Landau & Lifshitz, 1987;Whitham, 2011). In this section we are not covering the whole range of nonlinear effects described in literature but focus only on a nonlinear energy interaction effect, which is responsible for the energy exchange between different frequencies. When the nonlinearity is quadratic, to the first approximation, the nonlinear wave propagating at frequency f = f 0 will generate high-order harmonics at frequencies f = f 0 ± k · f 0 , where k = 1,2,3,… is an integer number. If, for example, we are considering an attenuation factor for the frequency f = 2f 0 , it may appear that the attenuation factor for this frequency component is negative (Q −1 < 0), because the energy is transferred from the base frequency f = f 0 to the higher-order harmonic f = 2f 0 ; therefore, the amplitude of frequency component f = 2f 0 is increased. The analytical solution of Apostol (2003) may explain more details about the growth of wave amplitudes at frequencies f = 2f 0 and f = 3f 0 . If two nonlinear waves are propagating at frequencies f 1 and f 2 , the nonlinear interaction will lead the following frequencies in the wavelet spectrum f i ± kf j , where i & j = 1 or 2, while k is an integer number. When the nonlinearity is cubic, the nonlinear interaction will lead the following frequencies in the wavelet spectrum f i ± 2kf j . Several scholars have shown that additional frequencies (harmonics) are generated by propagation of seismic waves through a hydrocarbon-saturated reservoir (e.g., Johnson et al., 1996Johnson et al., , 2004Khan & Khan McGuire, 2005;Zhukov et al., 2007). Moreover, a negative attenuation is frequently reported in literature (e.g., Jannsen et al., 1985;Mateeva, 2002;Matsushima, 2006), and it is referred to as "nonphysical" effect, because it is challenging to explain it using a linear theory of poroelasticity. Mateeva (2002) proposed different reasons for the "nonphysical" negative attenuations, such as scattering effects, ambient noise, spectral distortion by widowing source conditions, geophone coupling, and the choice of receiver separation. However, Matsushima (2006) has demonstrated that even after subtracting the scattering attenuation from the total attenuation, the "physically unrealizable phenomenon" of negative intrinsic attenuation could not be removed from the vertical seismic profiles data recorded in methane hydrate-bearing sediments. The redistribution of elastic energy between different frequencies is frequently reported in literature during processing of time-lapse vertical seismic profile data from partially water and gas-saturated reservoirs (e.g., Castagna et al., 2003;Goloshubin et al., 1996;Goloshubin & Bakulin, 1998;Korneev et al., 2004). Furthermore, Goloshubin and Bakulin (1998) have demonstrated that it is not possible to explain the energy transfer using Biot's theory.
Next, we are going to demonstrate that our model predicts the nonlinear energy transfer between different frequencies. We will calculate the spectrum of the bulk-volume strain rate and compare it with published field data, showing low-frequency microtremor anomalies above hydrocarbon reservoirs (Dangel et al., 2003;Goertz et al., 2012;Riahi et al., 2012;Saenger et al., 2009;Steiner et al., 2008;Witten & Artman, 2011). Herein, we calculate the spectrum of the strain rate because seismic geophones are sensitive to the strain rate of the wave-induced deformation (e.g., Aki & Richards, 2002); therefore, the spectrum of the strain rate is more suitable (than  the spectrum of the strain) to compare with the field measurements. Figure 9a shows both the spectrum of applied stress perturbation Δσ at frequency f = 0.5 Hz, shown by red dashed curve and the spectrum of the volume strain _ Δϵ B , shown by black continuous curve, calculated after 50 oscillation periods. All input parameters are the same as for Figure 4. The spectrum is normalized to its peak value at f = 0.5 Hz. The spectrum of _ Δϵ B contains odd harmonics (0.5,1.5,2.5,3.5,…) Hz, typical for nonlinear waves with cubic nonlinearity. Published laboratory data show that odd harmonics tend to dominate in amplitude in rocks (e.g., Johnson et al., 1996). The effect of nonlinear energy transfer from lower to higher frequencies is frequently reported in literature as the low-frequency microtremor anomalies above hydrocarbon reservoir (Dangel et al., 2003;Goertz et al., 2012;Lambert et al., 2009;Riahi et al., 2012;Saenger et al., 2009;Steiner et al., 2008;Witten & Artman, 2011), which was first reported in the Russian literature in 1980s (see review in Nikolaevskiy, 2010;Kazantsev, 2018). Figure 9b shows an example of spectrum of the passive seismic wavefield vertical surface velocities. Red curve shows measurements over a known gas field; blue curve shows measurements over an area with no hydrocarbon potential (modified after Saenger et al., 2009). The figure shows a microtremor amplitude anomaly above gas field at the frequency range of 2-6 Hz. This amplitude anomaly is possible by the energy transfer from lower frequencies (0.01-1 Hz), where the low-frequency energy with a peak around 0.2 Hz is generated by ocean waves with much higher amplitude of energy spectra (Holzner et al., 2009). Ocean waves creating pressure forces on the sea floor (and shore), which will generate seismic waves that would travel across the Earth at frequencies (0.01-1 Hz). The energy of the low-frequency seismic wave, propagation through partially saturated reservoir, can be converted to higher frequencies due to the nonlinear amplitude dispersion effect caused by a static contact angle hysteresis effect. It must be noted here that the nonlinear energy transfer between different frequencies cannot be explained by visco-poro-elastic models (see review : Kazantsev, 2018). Broadhead (2010) writes that he was not able to reproduce the amplification effect by considering the capillary pressure in the oscillating oil bubble inside the pore. This author did not consider the contact line pinning and the contact angle hysteresis effects in his model. In our model, if we set θ a = θ r , the contact line pinning force will disappear, and the energy transfer between different frequencies will also disappear, while the capillary pressure will not disappear. Thus, this effect is related to both to the capillary pressure and to the contact angle hysteresis effects. Furthermore, Kazantsev (2018) argues that none of the available rock physics models can explain the nonlinear energy transfer from low to high frequencies; therefore, this effect is classified in his dissertation as a hypothetical effect. At the same time, Kazantsev argues that by using this "hypothetical" nonlinear energy transfer between different frequencies, it is possible to explain the low-frequency microtremor anomalies above hydrocarbon reservoirs. According to our model, this effect is highly sensitive to the wettability of the rock. Effects of wettability on seismic wave velocities were observed more than 60 years ago, after Wyllie et al. (1958), who were the first (to our knowledge) to report this effect (see also Waite et al., 1997;Wang et al., 2015;Rozhko, 2019Rozhko, , 2020; however, the wettability effects are not mentioned in the literature review by Kazantsev (2018). If the  rock is mixed wet (when contact angles are close to 90°), then equation (46) would predict very high value of the critical wave amplitude (Δσ c ), required for the contact line slippage. If the critical wave amplitude is higher than the amplitude of the natural noise, the low-frequency microtremor anomaly will not be observed above hydrocarbon reservoir. Thus, the low-frequency microtremor anomalies above hydrocarbon reservoirs should depend on the amplitude of natural and anthropogenic noise, which may vary with time due to weather (storm) and different anthropogenic activity during day and night. These temporal variations in amplitude of microtremor anomaly are frequently observed in the field (e.g., Ali et al., 2010b;Hanssen & Bussat, 2008).

Conclusions
In this paper we extended the Gassmann's theory by considering the interfacial phenomena effects during wave-induced two-phase fluid-flow in partially saturated rocks. The Gassmann's theory is extended by three additional terms: (1) interfacial tension between immiscible fluids, (2) wettability (advancing and receding contact angles), and (3) pore size distribution or capillary pressure, described by Leverett J-function/ Brooks-Corey correlation. Here the capillary pressure can be different during drainage and imbibition; therefore, our model also predicts the effect of different fluid distributions (during drainage and imbibition) on seismic velocity.
Also, the fluid saturation that appears in the Leverett J-function is only related to the number of pores, saturated by the wetting phase, and does not include the deformation of rock. Our model allows to calculate wave-induced changes of fluid saturation due to both rock deformation and alteration of pore numbers, saturated by the wetting phase.
Our model predicts a nonzero attenuation of elastic energy at the zero-frequency limit due to a static contact angle hysteresis effect. Contrary all conventional models, based on viscosity dissipation, would predict no fluid effect on the seismic wave attenuation at the zero-frequency limit.
The nonlinear or amplitude dispersion effect, predicted by our model, depends nonmonotonically on the wave amplitude. When a wave amplitude is smaller than a critical amplitude (predicted by the model), then there is no nonlinearity in the system. Contact lines remain pinned while interface menisci are bending during wave-induced two-phase fluid flow. The nonlinearity reaches its maximum value when the wave amplitude is about two times larger than a critical amplitude. In this case the nonlinearity is controlled by bending and slipping of interface menisci at the pore scale. The nonlinearity decreases when the wave amplitude is much larger than a critical amplitude, because the contact line slippage dominates over bending of interface menisci.
Due to nonlinear amplitude dispersion effect, our model predicts the nonlinear energy exchange between different frequencies, frequently observed during processing of seismic data from partially saturated rocks.