Effects of pore fluids on quasi‐static shear modulus caused by pore‐scale interfacial phenomena

It is evident from the laboratory experiments that shear moduli of different porous isotropic rocks may show softening behaviour upon saturation. The shear softening means that the shear modulus of dry samples is higher than of saturated samples. Shear softening was observed both at low (seismic) frequencies and high (ultrasonic) frequencies. Shear softening is stronger at seismic frequencies than at ultrasonic frequencies, where the softening is compensated by hardening due to unrelaxed squirt flow. It contradicts to Gassmann's theory suggesting that the relaxed shear modulus of isotropic rock should not depend upon fluid saturation, provided that no chemical reaction between the solid frame and the pore fluid. Several researchers demonstrated that the shear softening effect is reversible during re‐saturation of rock samples, suggesting no permanent chemical reaction between the solid frame and the pore fluid. Therefore, it is extremely difficult to explain this fluid–rock interaction mechanism theoretically, because it does not contradict to the assumptions of Gassmann's theory, but contradicts to its conclusions. We argue that the observed shear softening of partially saturated rocks by different pore fluids is related to pore‐scale interfacial phenomena effects, typically neglected by the rock physics models. These interface phenomena effects are dependent on surface tension between immiscible fluids, rock wettability, aperture distribution of microcracks, compressibility of microcracks, porosity of microcracks, elastic properties of rock mineral, fluid saturation, effective stress and wave amplitude. Derived equations allow to estimate effects of pore fluids and saturation on the shear modulus and mechanical strength of rocks.

2. unrelaxed shear modulus of saturated rock is higher than shear modulus of dry rock μ sat > μ dry in the high frequency limit (e.g. Berryman 2005;Mavko et al. 2009).
There are several experimental studies which demonstrate a violation of our theoretical expectations due to socalled shear softening effect, when μ sat < μ dry , which was observed both at seismic and ultrasonic frequencies (e.g., Khazanehdari and Sothcott 2003;Adam, Batzle and Brevik 2006;Adam et al. 2009;Fabricius, Bächle and Eberli 2010;Adam and Otheim 2013;Bauer et al. 2016;Mikhaltsevitch, Lebedev and Gurevich 2016). Shear softening is stronger at seismic frequencies than at ultrasonic frequencies, where the softening is compensated by hardening due to unrelaxed squirt flow. Adam et al. (2006Adam et al. ( , 2009 reported shear modulus, measured on isotropic carbonate rocks at seismic and ultrasonic frequencies. Shear softening was observed for all samples at seismic frequencies and low confining pressure, whereas at high confining pressure shear softening is negligibly small for most of the samples. However, samples revealed either softening or stiffening at ultrasonic frequencies. Mikhaltsevitch et al. (2016) report both stiffening and softening of shear modulus both at seismic and ultrasonic frequencies for carbonate rocks. Furthermore, they demonstrated that stiffening or softening behaviour of shear modulus may depend on the pore fluid, which cannot be explained by Gassmann's theories. Vo-Thanh (1995) came to the same conclusion regarding the effect of pore fluid chemistry and saturation on acoustic velocities ultrasonic frequencies. Yin et al. (2019) report the shear softening of clay bearing sandstone at seismic and ultrasonic frequencies. Bauer et al. (2016) report experimental data at seismic and ultrasonic frequencies for sand, sandstone and shales suggesting that the Gassmann equations underestimate the observed velocity dispersion. Diethart-Jauk and Gegenhuber (2018) reported the shear softening behaviour at ultrasonic frequencies for different lithologies: limestone, dolomite, quartzite and basalt. Adam and Otheim (2013) reported shear softening effect observed at seismic frequencies on basalt rock and negligibly small stiffening effect at ultrasonic frequency. Lebedev, Wilson and Mikhaltsevitch (2014) reported softening of shear modulus of limestone sample observed at seismic frequencies.
Furthermore, several researchers reported the softening of elastic moduli observed during triaxial testing at large strains (e.g. Tutuncu, Podio and Sharma 1998;Risnes and Flaageng 1999;Baud, Zhu and Wong 2000;Risnes et al. 2003Risnes et al. , 2005David et al. 2015). Shear softening at ultrasonic frequencies for carbonate rocks were reported by Assefa, McCann and Sothcott (2003), Japsen et al. (2002), Røgen et al. (2005) and Sharma et al. (2013). Both shear softening and shear stiffening at ultrasonic frequencies for carbonate rocks were reported by Baechle et al. (2005Baechle et al. ( , 2009, Verwer, Braaksma and Kenter (2008), Fabricius et al. (2010), Vega, Prajapat and Al Mazrooei (2010), Verwer et al. (2010), Regnet et al. (2015) and Gegenhuber (2015), and for sandstone samples by Khazanehdari and Sothcott (2003), Bhuiyan and Holt (2016), Li et al. (2017Li et al. ( , 2018. To investigate if the reduction of frame modulus in carbonate rock is driven by chemical interaction, Baechle et al. (2005) conducted series of experiments in which saturation of carbonate samples were changed several times. They demonstrated that experimental results are reversible during resaturation of rock samples. These results imply that there is no permanent frame moduli reduction due to chemical interaction of pore fluid with the rock frame. Risnes and Flaageng (1999), Risnes et al. (2003Risnes et al. ( , 2005 obtained the same conclusion suggesting that the elastic moduli and strength are reversible during re-saturation of chalk with different pore liquids. These experiments suggest that dissolution and precipitation of calcite is not significant on the time scale when these experiments are conducted (Risnes and Flaageng 1999;Risnes et al. 2003Risnes et al. , 2005Baechle et al. 2005). Several researchers pointed out that the shear softening can be explained by chemical reaction of water with clay minerals (Yin et al. 2019). The chemical interaction of water with clay is well-documented in literature (e.g., Fjaer et al. 2008); however, in our opinion it will lead to irreversible effect, and thus not relevant to explain the reversible effect of pore fluid, reported by Risnes and Flaageng (1999), Risnes et al. (2003Risnes et al. ( , 2005 and Baechle et al. (2005). Murphy (1984) and Murphy, Winkler and Kleinberg (1986) conducted series of experiments to demonstrate that the water softening effect is taking place even in unconsolidated sandstones, which excludes watersoftening mechanisms, such as swelling of clay cement, osmotic suction and dissolution and precipitation of calcite. The literature outlined above shows that physicochemical reasons for shear softening are not clear as it was stated in several publications (e.g. Japsen et al. 2002;Adam et al. 2006;Fabricius et al. 2010), especially when there is no evidence to chemical interaction of pore fluid with the rock frame.
In this paper, we demonstrate that the shear softening effect is related to the effect of hysteresis of liquid bridges in rocks containing partially saturated cracks. Here we are focusing on shear moduli of partially saturated rock, whereas bulk moduli were investigated by our previous publication (Rozhko 2019). The hysteresis of liquid bridges (described in the next section) is related to interface phenomena effects between immiscible fluid phases (such as water, gas or oil) and the solid inside compliant pores (cracks). The article is organized as follows. First, we introduce the hysteresis of liquid bridges effect in isotropic rock containing partially saturated cracks and argue that due to this effect the relaxed shear modulus is sensitive to the fluid saturation. Second, we derive mathematical equations to describe this effect numerically. Third, we introduce Betti-Rayleigh-reciprocity theorem (e.g., Schmeling 1985) for calculation of the effective shear modulus and attenuation of partially saturated rock. Then, we present numerical investigation of different input parameters of the model on the shear softening effect and shear attenuation. Next, we made the dimensional analysis of derived non-linear equations to demonstrate that the shear softening effect can be described by relatively simple equation. And finally, we discuss applications of the theory to reservoir rocks with log-normal distribution of micro-crack aperture. Mavko et al. (2014, p. 138) suggested the following explanations for why Gassmann's relations only work at low frequencies and for isotropic rock. Imagine an isotropic sample of rock with cracks at all orientations. Under 'pure shear' loading, there is no volume change of the rock sample or the pore space, because some cracks open, whereas others close. According to Fig. 1(a), crack #1 decreases in volume, its pore pressure locally increases if the fluid cannot flow out of the crack, whereas crack #2 increases in volume, its pore pressure locally decreases if the fluid cannot flow into the crack. If the frequency is too high, there is a tendency for local pore pressures to increase in some pores and decrease in others: hence, the shear modulus depends on the fluid bulk modulus. Greater bulk modulus of fluid makes it more difficult to deform undrained pores, thus the shear modulus increases with increase of fluid bulk modulus at high frequencies. However, if the frequency is low enough, the fluid has time to flow and adjust: there is no net pore volume change and thus no change of pore pressure, therefore the shear modulus is independent of the fluids. If the rock is anisotropic, under pure shear conditions, some cracks will be more open in one direction, whereas others less closed in perpendicular direction: hence, the shear modulus depends on the fluid bulk modulus.

F L U I D E F F E C T O N R E L A X E D S H E A R M O D U L U S D U E T O H Y S T E R E S I S O F L I Q U I D B R I D G E S
Next, we apply the same approach as described above by Mavko et al. (2014) to demonstrate the opposite, that is why Gassmann's relations will not work at low frequencies for isotropic, but partially saturated rock (see Fig. 1b). The wetting phase (water) occupies narrow parts of the crack, whereas the non-wetting phase (gas or oil) occupies wide (central) parts of the crack. At the equilibrium condition, the fluid pressure in each fluid phase is different due to interface tension between immiscible fluids. It is not difficult to understand that the low frequency shear modulus of isotropic partially saturated rock would depend on the fluid saturation due to physical phenomenon: hysteresis of liquid bridges (e.g.: De Souza et al. 2008;Chen, Amirfazli and Tang 2013;Zhang 2016;Shi et al. 2018). Due to this phenomenon, the positive change of pore volume in one direction will not be compensated by a negative change of pore volume in the perpendicular direction and hence the shear modulus will depend on the fluid content. Figure 2 shows a typical experimental setup used to investigate the hysteresis of liquid bridges between parallel plates (e.g. De Souza et al. 2008;Chen et al. 2013;Zhang 2016;Shi et al. 2018). The aperture between two plates is w, the diameter of wetted area is l and F is the force between two plates. A liquid drop is placed between two plates. The upper and the lower contact angles are θ a and θ r (i.e. advancing and receding angles). Here the difference between advancing and receding contact angles is called a contact angle hysteresis. It depends on the contact line motion velocity. At zero velocity, a spectrum of static contact angles is observed (Bormashenko 2013a,b), which can be as large as tens of degrees (Ethington 1990). This effect is called a static contact angle hysteresis. The advancing contact angle increases with the contact line advancing velocity, whereas the receding contact angle decreases with the contact line receding velocity. This effect is called a dynamic contact angle hysteresis. In this paper, we address the Gassmann's theory which is a quasi-static theory. Thus, small dynamic corrections to the contact angle hysteresis during a quasi-static contact-line motion can be neglected. The change of aperture within a certain frequency ω is w. Aperture changes induce changes of the force F and the diameter of wetted area l. Experiments are typically conducted under conditions when frequency ω is low enough so that the viscous forces can be neglected. Figure 3 shows schematically typical experimental results

Figure 3
Effect of oscillation amplitude on hysteresis of liquid bridges (e.g. De Souza et al. 2008;Chen et al. 2013;Zhang 2016;Shi et al. 2018). If the oscillation amplitude is small, there is no hysteresis (a), whereas if the oscillation amplitude is large, there are four stages of hysteresis (b). (a1) and (b1) capillary force versus aperture (between two parallel plates); (a2) and (b2) contact angles versus aperture; (a3) and (b3) diameter of the wetted area (contact line displacement) versus aperture. Note here that Fig. 3 showing the hysteresis of liquid bridges (e.g. De Souza et al. 2008;Chen et al. 2013;Zhang 2016;Shi et al. 2018).  and (b1-b3) shows typical experimental results, obtained for different amplitudes of the aperture change w. Figure 3(a1,b1) shows capillary force versus aperture between parallel plates; Fig. 3(a2,b2) shows contact angles versus aperture; and Fig. 3(a3,b3) shows the diameter of the wetted area versus aperture. When the amplitude of the periodic aperture change is small, the contact line is pinned (i.e. l = 0, see Fig. 3a3) and the contact angle is changing within the range θ r < θ i + θ < θ a (Fig. 3a2), θ i is the initial contact angle, which will be introduced later in this section. The change of the capillary force shows no hysteresis (Fig. 3a1), when the periodic amplitude of aperture change is small. The loading path (Fig. 3a1) coincides with unloading path. If the amplitude of the aperture change is increased, the amplitude of the contact angle change will increase as well. The amplitude of the contact angle change cannot exceed this range θ r ≤ θ i + θ ≤ θ a , thus when the aperture change amplitude is sufficiently large, the contact angle will change between advancing and receding angles, as shown in Fig. 3(b2). At the time moment when the contact angle reaches advancing or receding angle, the contact line will slip (Fig. 3b3). Overall, the hysteresis of liquid bridges can be divided into four continuous stages, when the amplitude of the aperture change is sufficiently large, as shown in Fig. 3(b1-b3): 1. Pinning (stretching) 2. Slipping (receding) 3. Pinning (compression) 4. Slipping (advancing) When the contact angle is greater than the receding angle, the aperture increase will result in an increase in the force due to pinning (Stage 1). This will be accompanied by a reduction in the contact angle until the receding angle is achieved when the force starts to decrease and the contact line starts to slip inward (Stage 2). If the aperture starts to decrease, the contact angle begins to increase until it reaches the advancing angle (Stage 3). In this stage, the pinning stage occurs again, which correspond to the reduction of capillary force. If the aperture keeps decreasing to the initial aperture, the contact line will slip outward with the contact angle equal to the advancing angle (Stage 4). The hysteresis of the capillary force leads to the energy dissipation due to the contact line friction mechanism. The contact line friction is the dominant mechanism of energy dissipation in these experiments. The intermolecular forces acting between molecules of the solid and those of the liquid, which pin the contact line to the substrate, are responsible for the contact line friction, which occurs not over the entire solid-liquid interface, but only at the three-phase line (Yaminsky 2000;Bormashenko 2013a,b). Contact line friction occurs due to slippage and rolling of fluid molecules over the surface of the solid at the contact line location, that is in one dimensional (Ren and Weinan 2007). It is interesting to note here that viscous dissipation occurs in three dimensional due to interaction of fluid molecules moving with different velocities, whereas frictional dissipation between two solids in a contact occurs over a certain surface area in two dimensional. Thus, fundamentally, the contact line friction is different from other dissipation mechanisms.
Next, let us discuss the choice of initial contact angle (θ i ). In calculations, we assume that the initial contact angle is equal to the most stable Young's angle θ i = θ Y . It is important to note here that any contact angles within this range θ r < θ < θ a are equilibrium contact angles; however, there is a most stable configuration of the equilibrium contact angle (e.g. Ruiz-Cabello, Rodríguez-Valverde and Cabrerizo-Vílchez 2014). The Young's contact angle is related to three coefficients of interfacial tension formed by the fluid-fluid interface with the solid surface (de Gennes, Brochard-Wyart and Quéré 2013). The equilibrium Young's contact angle can also be calculated from the advancing and receding contact angles, as was shown theoretically by Tadmor (2004) and confirmed experimentally by Chibowski (2008).
Next, after introducing the hysteresis of liquid bridges phenomenon, we come back to the Fig. 1(b), showing two orthogonal partially saturated cracks under pure shear stress perturbation. It is clear without any mathematical calculations that due to hysteresis of the capillary forces ( Fig. 3b1) the net change of pore volume in isotropic partially saturated rock will not be equal to 0. Hence, the shear modulus depends on the fluid saturation. The net change of the pore volume will be equal to 0 only if there is no hysteresis of capillary forces, that is when the contact line is pinned for small w. Hence, the net volume change of the partially saturated rock depends on the wave amplitude and will vanish in the limit when the wave amplitude is 0. It will be no dependence of shear modulus on the fluid if the surface Figure 4 (a) Initial equilibrium configuration of two identical (plane-strain) cracks: σ initial hydrostatic far-field compressive stress; p we and p nw initial pressure in the wetting and non-wetting fluid phases; θ initial contact angle for the wetting phase; c initial contact line location. (b) Perturbation of the initial equilibrium configuration by the pure shear stress perturbation ( τ ). Perturbations of p we and p nw are identical in both (interconnected) cracks (in the low the frequency limit), whereas perturbations of contact angles θ i and contact line locations c i are different in both cracks. tension between immiscible fluid is 0, thus the change of the capillary force will be 0, that is F = 0 for any w.

M A T H E M A T I C A L F O R M U L A T I O N
In this section, we develop mathematical model, describing the hysteresis of liquid bridges in the rock containing two partially saturated cracks, shown in Fig. 4(a). The aim of this paper to demonstrate that the shear softening is caused by the hysteresis of liquid bridges in compliant cracks; that is why we are focusing only on the minimum number of input parameters required to describe this effect. Thus, we neglect porosity of channels which connect two cracks together, because otherwise we would also need to consider the distribution of pore throats and compressibilities of pore throats. Furthermore, real rocks contain cracks-like pores of different aspect ratio which are closed at different confining stress (e.g. Zimmerman 1990), while for simplicity we are considering only two identical cracks, which is sufficient for explanation of the shear softening effect by the hysteresis of liquid bridges without making the model too complex. Two identical plane-strain cracks (Crack # 1 and Crack #2) are aligned along x and y directions, respectively. Initial geometries of two identical cracks are described by deformable elliptical cavities with semi-axes a and b. In calculations, we consider very narrow cracks with a â b. The distance between two cracks is sufficiently large, so that the deformation and stress concentrations around one crack is not affected by deformation and stress concentrations around another crack (e.g. Zimmerman 1990). The deformation and stress concentration around cracks are affected only by fluid pressure inside cracks and stresses on the external boundary of the representative elementary volume (REV). The wetting fluid phase occupies thin parts of the crack (i.e. tips), whereas the non-wetting fluid phase occupies wide parts of the crack (centre). Pressures in the non-wetting ( p nw ) and wetting ( p we ) fluid phases are different due to interfacial tension and wettability effects and denoted here as capillary pressure ( p cap ) (e.g. Barenblatt, Entov and Ryzhik 1990): Initial saturation of two cracks is identical, as shown in Fig. 4(a). The fluid pressure acting on the surface of crack #1 is p nw if |x| ≤ c and p we if c < |x| ≤ a. Analogically, the fluid pressure acting on the surface of crack #2 is p nw if |y| ≤ c and p we if c < |y| ≤ a (Fig. 4a). The coordinate c defines the location of the contact line, as shown in Fig. 4(a). We consider a uniform initial far-field confining stress σ , acting on the external boundary of the REV, as shown in Fig. 4(a). Nonuniform initial far-field confining stress will result in different compressibilities along x and y directions, but we are not interested in this anisotropic case, because it contradicts to assumptions of Gassmann's theory, as it was discussed around Fig. 1(a) in previous section. Cracks are assumed to be connected hydraulically, thus initial fluid pressure in each fluid phase is identical in two cracks. Furthermore, we consider a quasi-static deformation of the external boundary of REV, so that there is enough time for fluid pressure in each fluid phase to equilibrate between two cracks. The initial contact angle is equal to Young's angle θ = θ Y in both cracks. Because the initial far-field normal stress, the initial fluid pressure and initial contact angles are identical in two identical cracks, the initial saturation degree is also identical in two cracks. Rozhko (2016), Rozhko and Bauer (2018) and Rozhko (2019) derived equations, describing the initial equilibrium saturation and the capillary pressure of partially saturated crack. The capillary pressure inside the equilibrium crack is calculated as follows (Rozhko 2016): where p cl is the crack closure pressure, calculated as (Rozhko 2016): where μ and v are the Shear modulus and Poisson's ratio of the rock mineral. The capillary pressure in equation (3) is twoway coupled to the deformation of the crack aperture. This equation follows from the Laplace equation and the analytical expression for the crack aperture at the contact line location (Rozhko 2016). In equation (3), σ n is the far-field normal stress (i.e. normal to the crack surface), equal to σ at the initial condition, shown in Fig. 4(a). Here we are using a sign convention when compressive stresses, strains and displacements are negative, whereas compressive pressure is positive. In equation (3), γ is the surface tension between immiscible fluids and the angle β defines the location of the contact line: The total volume of partially saturated crack (V tot ) is calculated as (Rozhko 2016): Although the volume of the wetting fluid phase inside the partially saturated crack (V we ) is calculated as (Rozhko and Bauer 2018;Rozhko 2019): And More details about derivation of above equations can be found in Rozhko (2016) and Rozhko and Bauer (2018).
The volume of the non-wetting phase inside the crack is calculated as The wetting-phase saturation of the crack (S we ) is defined as Above equations are applicable for the following saturations range of the crack (see details in Rozhko and Bauer 2018 In calculations, we consider very narrow cracks, a b, thus above equations are applicable almost for the entire range of saturations. Next, following Mavko et al. (2014, p. 138), we consider the perturbations of the equilibrium state caused by a change in a far-field stress by 'pure shear' stress perturbation ( τ ) (Fig. 4.b). The amplitude of stress perturbation induced by seismic waves is very small (comparing to initial confining stress) typically around τ ∼ 10 2 − 10 4 Pa. This pure shear stress perturbation induces perturbations of normal stresses: σ n,1 = τ in the crack #1 and σ n,2 = − τ in the crack #2. Here by sub-indexes i = 1 or i = 2, we will denote cracks #1 or crack #2, respectively. According to Fig. 4(b), the far-field 'pure shear' stress perturbation ( τ ) will also induce perturbations to the following parameters: (1) perturbation of the fluid pressure in the wetting phase p we , (2) perturbation of the fluid pressure in the non-wetting phase p nw , (3) perturbation of the contact line location c i and (4) perturbation of the contact angle θ i . Here, we are considering the low frequency limit, when the fluid can flow between two cracks to equilibrate any gradients of fluid pressure in each fluid phase. Thus, the perturbation of fluid pressure in the wetting phase would be p we and is identical in both cracks, as well as perturbation of the fluid pressure in the nonwetting phase would be p nw and is identical in both cracks. Thus, perturbation of p cap = p nw − p we is also identical in two cracks. Although perturbations of other parameters, shown in Fig. 4(b), are not identical, that is c 1 = c 2 , θ 1 = θ 2 and σ n,1 = σ n,2 . Here, to be more specific, we can say that σ n,1 = − σ n,2 due to pure-shear perturbation of boundary conditions, whereas c 1 = − c 2 and θ 1 = − θ 2 due to non-linearity effects, as it will be discussed below.
To find perturbations of these parameters we need to consider perturbations of equilibrium parameters, such as perturbations of capillary pressure (in each crack), perturbations of crack volumes and perturbations of wetting phase volumes in both cracks. Because the hysteresis of liquid bridges is a highly non-linear effect, small stress perturbations of amplitude τ ∼ 10 2 − 10 4 Pa are enough to cause this non-linearity due to deformation of highly compliant cracks. For example, recent publication of Rozhko and Bauer (2018) have shown that a small stress perturbation (∼ 10 2 − 10 4 Pa) are sufficient to cause large changes in the contact angle θ of the order of a few degrees or even larger. Thus, we cannot use linear (Taylor's) expansions of equations (3), (6) and (7), because perturbations of θ i are large. However, we can consider a much smaller increment of the wave amplitude T: where τ is the wave induced pure shear stress perturbation, defined in Fig. 4(b), whereas N iter is the number of iteration steps (N iter 1). The number N iter can be sufficiently large so that changes of all parameters, including θ i will be small. Typically, N iter ∼ 10a/b is sufficient to insure small deformations around crack aperture per stress increment T.
For crack # 1, Taylor's expansion of equilibrium equations (3), (6) and (7) becomes and Similarly, for crack #2, Taylor's expansion of equilibrium equations (3), (6) and (7) becomes and In these equations, changes of normal stresses in two cracks are σ n,1 = T and σ n,2 = − T. All partial derivatives are calculated analytically and can be found in our previous publication Rozhko (2019).
Next, we are considering undrained boundary conditions, when there is no flow of the wetting and non-wetting phases in or out of the external boundary of REV, while fluids can flow freely inside REV between two cracks. Thus, in our model, the masses of the wetting fluid and non-wetting fluid remain the same in the REV. In this case, changes of volumes of the wetting and non-wetting fluid phases are related to changes of pressures in those phases via the bulk moduli of the wetting (K we ) and non-wetting (K nw ) fluid phases as follows, where those changes are driven by small increment of the wave amplitude T: and It is interesting to note here that if we set K we = 0 and K nw = 0 in equations (20) and (21), we will get drained boundary conditions for REV, because perturbations of fluid pressure in each fluid phase will be equal to 0. Furthermore, one of the fluid phase can be under drained boundary conditions, whereas another fluid phase can be under undrained boundary condition. It can be modelled by setting corresponding fluid bulk modulus to 0 in equations (20) or (21). Thus, equations (20) and (21) embrace different types of boundary conditions for two immiscible fluids in REV: undrained, drained and mixed-mode boundary conditions.
Equations (14)- (21) are eight equations with 10 unknown parameters: p cap , p we , θ 1 , β 1 , V tot,1 , V we,1 , θ 2 , β 2 , V tot,2 and V we,2 . However, the number of unknown parameters will be reduced to eight, if we are considering either contact line pinning or the contact line motion conditions in two cracks. For example, the change of angle β i in equations (14)- (19) is related to the change of the contact line location, according to equation (5). If the contact line is pinned, we have β i = 0 and θ i = 0; whereas if the contact line moving, we have β i = 0 and θ i = 0. Thus, the system of equations (14)-(21) will have a unique solution if we are considering either the contact line pinning or the contact line motion conditions. In order to calculate response, produced by the stress perturbation τ , we need to solve the system of equations (14)-(21) N iter times (iterative solution). After each step, the following variable parameters need to be updated: In above equation, the super-script (k) denotes the iteration step number, where k = 0 denotes the initial conditions. We do not need to update other parameters, such as p (k) cap , i , according to equations (3), (6) and (7), and thus can be calculated analytically. After applying equation (3) for both cracks, we will get p (k) cap,1 = p (k) cap,2 , as expected. After applying equations (6) and (7), we will get we,2 . After each step, all partial derivatives which appear in equations (14)-(21) must be re-calculated as well, because they depend on p (k) we , σ (k) n,i , β (k) i and θ (k) i . These parameters are identical for both cracks only at the initial conditions when (k = 0). Furthermore, volumes of the wetting phase V we,i and crack volumes V tot,i should be updated after each step, in equations (20) and (21).
Next, we discuss more details about modelling the conditions when the contact line is pinned or moving. First, we define the contact line displacement, using equation (5) as follows: According to Fig. 4, the contact line moves to the receding direction if c i > 0 and to advancing direction if c i < 0. Thus, the contact line in crack # i is pinned at the iteration step k if one of the following three conditions is taking place: . The contact line in crack # i at the iteration step k will be moving in advancing direction when (θ (k−1) i = θ a and c i < 0) and will be moving in receding direction when (θ (k−1) i = θ r and c i > 0). Next, we can derive analytical solutions to the system of linear equations (14)-(21) under four different conditions. Case 1: The contact line is pinned in both cracks: Here for simplicity we defined volumes of the wetting phase inside two cracks as and volumes of the non-wetting phase as Case 2: The contact line is moving in crack # 1, but pinned in crack # 2: Case 3: The contact line is moving in crack # 2, but pinned in crack # 1: Case 4, the contact line is moving in both cracks: Next, after solving the system of equations (14)-(21) N iter times, we will find the response induced by the stress perturbation τ , that is p There is no explicit time dependence in this algorithm; however, the wave induced stress perturbation τ can be a function of time. When we discretize τ on much smaller stress increments, using equation (13) the time axis can also be discretized in the same manner accordingly. Thus, implicitly we can simulate any time-dependent signal using this algorithm by discretizing the wave-form on subintervals with subsequent numerical integration of equations of this paper.

S H E A R M O D U L U S A N D A T T E N U A T I O N I N P A R T I A L L Y -S A T U R A T E D R O C K
In this section, a set of equations for the relaxed moduli and attenuation factors for a material containing two partially saturated cracks (Fig. 1b) will be given. Analytical solutions from the previous section describing wave-induced perturbations of fluid pressure and crack volumes will be utilized.
The effective elastic moduli of a material containing arbitrary inclusions can be obtained by using the Betti-Rayleigh-reciprocity theorem (Schmeling 1985;Mavko and Jizba 1991). It connects two elastic states of equilibrium of a linear elastic body. If these two states 1 and 2 are represented by the surface displacements u 1 and u 2 and the surface stress vectors or tractions T 1 and T 2 acting on the body, the theorem can be written as where the integration must be carried out over the total outer and inner surfaces F (e.g. Schmeling 1985). The shear modulus of the rock containing two partially saturated cracks (Fig. 1b) can be calculated by the following approach. Figure 5 shows two states of stresses and displacements, which can be connected using reciprocity theorem. Figure 5(a) shows the body containing cracks loaded by an external 'pure shear' tractions τ and internal waveinduced perturbations of fluid pressure p f l,i inside crack #i. The change of fluid pressure in each crack depends on p we and p cap , which are identical in both cracks and also depends on the change of the contact line location c i , which is different in both cracks, thus volume-average changes of fluid pressure are not identical in both cracks. p f l,i is not uniform along the crack length, as it is shown in Fig. 5(a) by p f l,1 (x) and p f l,2 (y) for cracks # 1 and #2, respectively. U and u shows body surface and pore wall displacements. In thin cracks, it is sufficient to consider only normal components of the crack wall displacements ( u y,1 (x) and u x,2 (y)), because tangential components are negligibly small. In Fig. 5(b), the body containing cracks is loaded by an external pure 'shear tractions' τ , the same as in Fig. 5(a). In Fig. 5(b), the thin cracks are loaded with tractions, the same as external tractions. Thus, the fluid pressure in crack # 1 is p f l,1 = − τ and inside crack #2 is p f l,2 = τ . Here we remind for clarity that we are using a sign convention when compressive stresses are negative while compressive pressure is positive. The resulting displacements U 0 and u 0 (Fig. 5b) are the same if cracks would contain solid matrix material (e.g. Schmeling 1985). Two states of stress in Fig. 5 can be combined by the reciprocity theorem giving: Where, according to Figure 5, shear tractions are given by n is the surface normal unit vector, F ext is the total outer surface of representative elementary volume (REV) and F cr,i is the surface of crack #i. In equation (37), the summation is taken over two cracks. The specific strain energy represented by the first integrals on either side of equation (37) can be taken to define the undisturbed shear modulus of rock mineral and effective shear modulus of REV, (e.g. Schmeling 1985): and where V REV is the volume of REV, which can be related to the initial crack porosity (n c ) as follows: where the initial crack porosity is the porosity at zero effective stress. In equation (37), displacements of the crack wall u 0,i are the same if the crack would contain the solid rock mineral, whereas displacements u i represent displacements of open crack that is partially saturated with water and gas. Because partially saturated crack is much more compressible than crack 'filled' with the solid rock mineral, we have | u i | | u 0,i |, hence the integration term 2 i=1 F cr,i p f l,i (x)n · u 0,i dF cr in equation (37) can be neglected. Similar simplifications were also used in literature for more arbitrary pore geometries (e.g. Schmeling 1985).
In equation (40), the displacement on the outer boundary of REV (U) is unknown, however, using equations (37)-(41), we can still calculate the effective shear modulus of REV as follows: In equation (42) the traction τ is uniform inside the crack, thus the following simplification can be made: where the integral F cr,i u i dF cr represents the change in volume of crack # I due to deformation of crack aperture. The initial volume of crack # i is calculated using equation (6) and input parameters: at the iteration step k = 0, whereas the finale crack volume is calculated using input parameters in equation (6) at the iteration step k = N iter . Note here that the capillary pressure, which appears in equation (6), is calculated using equation (3) and the same input parameters: p (k) we,i , σ (k) n,i , β (k) i and θ (k) i . Thus, we obtain: where Hence, equation (45) combined with equation (42) becomes: The above equation is very general and applicable for calculation of relaxed and unrelaxed shear modulus (in low and high frequency limits) of partially saturated, fully saturated and dry rocks. This equation depends on the change of crack volumes with applied stress perturbations, whereas the change of crack volumes depends on the fluid content of cracks. Furthermore, this equation can be applied to calculate effective shear modulus at any iteration step (k), that is where Thus, the shear strain of REV can be calculated at any iteration step k using equation: Using equation (49), we can calculate the elastic strain energy of REV, induced by stress perturbation τ as follows: Note here that equation (50) does not represent the absolute elastic energy, but only the wave energy, caused by transient stress τ .
The attenuation of elastic wave energy can be calculated by the integration of the elastic energy over period of the elastic wave, that is Equation (51) suggests that elastic wave energy will be attenuated if there is a hysteresis of shear strain during loading and unloading cycles. The attenuation factor 1 Q is defined as a fraction of energy dissipated over wave period (Mavko et al. 2009;Müller et al. 2010).
Next, it will be interesting to compare results, derived for partially saturated rock with results, calculated for dry rock or for fully saturated rock with a single fluid phase. If we are considering a dry rock with two cracks or fully saturated rock (single phase) in the low frequency limit (relaxed modulus), equation (46) yields Equation (53) directly follows from equation (6) in which p we = p f l = const and p cap = 0 due to full saturation of cracks. Equation (53) does not depend explicitly on the effective stress σ = σ + p f l (negative in compression); however, this equation is applicable as long as cracks are open, when p cl + σ + p f l > 0. If cracks are closed, the rock will deform as if there are no cracks, thus effective rock moduli will be equal to moduli of the rock mineral, that is μ dry = μ min . If the rock contains many cracks of different aspect ratio, those cracks will be closed at different effective confining stress (e.g. Zimmerman 1990). Thus, apparent shear modulus of dry rocks containing cracks of different aspect ratio will depend on the effective confining stress (e.g. Zimmerman 1990). In this paper, we are considering only two identical cracks for simplicity. However, the model can be extended to more realistic case of many cracks with different length and aspect ratio. This, more realistic case will be considered at the end of this paper. Furthermore, equation (53) does not depend on the fluid compressibility, because in the low frequency limit the change of fluid pressure (induced by pure shear stress perturbation) will be 0 in isotropic rock. Thus, in equation (53), we re-derived conclusion obtained by Gassmann's that μ sat = μ dry , which is valid for isotropic rock at low frequency when the pore fluid does not interact with the rock frame.

N U M E R I C A L R E S U L T S A N D D I M E N S I O N A L A N A L Y S I S
In this section, we present numerical results for shear moduli (μ eff ) and seismic attenuation ( 1 Q ), calculated for representative elementary volume with two partially saturated cracks. In calculations, we will apply a periodic 'pure-shear'stress perturbation to the external boundary of representative elementary volume (REV), as shown in Fig. (6a). The first cycle is shown by continuous black curve, whereas subsequent cycles are shown by the red dashed curve in Fig. 6(a). Figure 6(b) shows calculated stress-strain hysteresis of REV with two cracks, partially saturated with water and gas system with 50% of water saturation (S we = 0.5). Other input parameters of the model are given in Table 1. This saturation degree (and input parameters from Table 1) corresponds to initial capillary pressure of p cap = 0.56 MPa, calculated using equations (3) and (11). Similar to Fig. 6(a), the first cycle of the stress-strain hysteresis is shown by continuous black curve in Fig. 6(b), whereas subsequent cycles are shown by the red dashed curve. Calculations show that subsequent cycles do not return to its initial value with γ REV = 0 and τ = 0. According to equation (49), the effective shear modulus of REV is inversely proportional to the slope of largest diagonal of the tetragon, shown in Fig. 6(a), whereas the attenuated energy is proportional to the area of the tetragon, according to equation (51). Figure 7 shows calculated changes of other parameters of the model. Alteration of fluid pressure in the wetting and non-wetting phases are shown in Fig. 7(a,b), respectively. Calculations show that amplitudes of the wave-induced fluid pressure perturbations are negligibly small, comparing to the wave amplitude. Similar to Fig. 6, the first cycle is shown by continuous black curve, whereas subsequent cycles are shown by the red dashed curve in Fig. 7. Figure 7(c,d) shows alterations of the contact line location in the first and second cracks, whereas Fig. 7(e,f) Figure 6 (a) Applied a periodic 'pure-shear' stress perturbation to the external boundary of REV. (b) Calculated stress-strain hysteresis of REV with two cracks, partially saturated with water and gas system with 50% water saturation (S we = 0.5). Other input parameters of the model are given in Table 1. In (a) and (b), the first cycle is shown by continuous black curve, whereas subsequent cycles are shown by the red dashed curve. See Fig. 7 for changes of other parameters of the model. shows alterations of the contact angles in the first and second cracks. Overall calculations of Figs 6 and 7 show that the hysteresis of liquid bridges is highly non-linear process, leading to the attenuation of the wave energy in the low (static) frequency limit. This non-linearity is highly sensitive to the wave amplitude, as it was discussed around Fig. 3. Next, we will investigate how the wave amplitude affect this nonlinearity and how this will affect the effective shear modulus and attenuation in partially saturated rocks. We will consider both gas-water and oil-water systems. Assumed input parameters for both systems are given in Table 1, showing that the surface tension between gas and water interface is about 2.5 times higher than the surface tension for oil and water interface. Also, the gas phase is less wetting phase comparing to oil phase, as it is reflected in assumed input contact angles for the water phase. Also, the assumed bulk modulus of oil is much greater than the bulk modulus of gas, as shown in 1 .2 GPa Surface tension between water and gas (γ ) 0 .073 Pa × m Surface tension between water and oil (γ ) 0 .029 Pa × m Advancing contact angle for wetting phase in water-gas system (θ a ) 5°R eceding contact angle for wetting phase in water-gas system (θ r ) 1°A dvancing contact angle for wetting phase in water-oil system (θ a ) 37°R eceding contact angle for wetting phase in water-oil system (θ r ) 35°I nitial crack porosity (n c ) 1 0 −3 Major semi-axis (a) 1 0 −3 m Initial aspect ratio (b/a) 1 0 −4 Effective stress (σ + p we ) −15 MPa Table 1. Figure 8(a) shows calculated effective shear modulus as a function of the wave-amplitude calculated for partially saturated rock with gas-water and oil-water systems, at water saturations S we = 0.3 and S we = 0.6. Dashed black curve in Fig. 8(a) shows the shear modulus of dry rock. Calculations show that the shear softening effect depends on the wave amplitude and on the type of fluid and saturation degree. When the wave amplitude is small, the shear modulus of partially saturated rock coincides with shear modulus of dry rock. At small wave amplitude, the contact line is pinned, whereas if the contact line is slipping, the shear modulus of partially saturated rock is decreased. Wang, Schmitt and Wang (2015) argued that, due to the slippage at the solid-fluid interface, the stiffness of the rock will be lower than predicted by Biot-Gassmann models, which assume Stoke's no-slip boundary conditions between fluid and solid. Our theoretical calculations support Wang et al. (2015) explanation, because the Stoke's no-slip boundary conditions are not applicable at the contact line location (Ren and Weinan 2007). The contact line motion occurs by slippage and rolling of fluid molecules over the surface of the solid at the contact line location (Ren and Weinan 2007), this leads to the shear softening, predicted only when the contact line is moving. Calculations of Fig. 8(a) show that the shear softening occurs at certain critical wave amplitude τ c . At this wave amplitude, the contact line starts to move causing both shear softening and attenuation of elastic wave energy by the contact line friction mechanism. Figure 8(a) shows that the shear softening is stronger for gaswater system as compared to oil-water system. Furthermore, Fig. 8(a) shows that the shear softening occurs at lower wave amplitude for gas-water system as compared to oil-water system. Figure 8(b) shows calculated attenuation factors for corresponding cases shown in Fig. 8(a). Calculations show that the attenuation is higher for gas-water system as compared to oil-water system. Next, we will perform the dimensional analysis to understand better the impact of different input parameters on the shear softening and on the wave attenuation in partially saturated rocks. The shear softening is highly non-linear process described by the system of coupled equations, which are solved iteratively. However, to understand the physics, we will neglect the non-linearity and iterations. First, we will estimate the critical stress perturbation required for the onset of the contact line motion τ c . To do so, we will use the system of linear equations (30) in which all coefficients are calculated at the initial conditions, that is at zero interaction step (k = 0), whereas the change of the contact angle, required for the onset of the contact line motion, is selected (approximately) as | θ 1 | = θ a −θ r 2 and θ 1 = − θ 2 . The initial contact angle is equal to Young's angle and calculated with equations (1), which predicts approximately intermediate angle between advancing and receding angles, that is θ Y ≈ θ a +θ r 2 . Thus, the change of initial angle by θ 1 or by θ 2 will initiate the contact line slippage. Because we neglect the non-linearity and iterations, the system of equations (30) yields the following simplified equations, for the critical wave amplitude, required for the onset of the contact line motion: Partial derivatives in this equation are calculated at initial conditions using equations (A37) and (A39) from Rozhko (2019).
Next, we will estimate the amplitude of shear softening using the system of equations (35) when the contact line is moving in two cracks. In this case, we will assume sufficiently large wave amplitudes ( τ c ) when the contact line motion is much larger than the amplitude of the meniscus bending.  (b) show calculated effective shear modulus and attenuation factor as a function of the wave-amplitude for the rock that is partially saturated with gas-water and oil-water systems, at water saturations of S we = 0.3 and S we = 0.6. Dashed black curve in (a) shows the shear modulus of dry rock. Other input parameters of the model are given in Table 1.
After neglecting the non-linear effects in the system of equations (35), we can derive the expression for the difference between V tot,1 and V tot,2 , which will be used in equation (46) to derive the effective shear modulus of partially saturated rock. After all simplifications using equations (35), (41), (46), (53) and equations (A40) and (A41) from Rozhko (2019), the following expression for the effective shear modulus of partially saturated rock is derived: All partial derivatives in equation (55)  ∂β ) has less transparent physical meaning, because the angle β is related somehow to the contact line location inside the crack, according to equation (5). Let us understand better the physical interpretation of equation (55). To do so, let us first consider the equation (6). From this equation, we can see that the total volume of partially saturated crack is controlled by the following effective stress (e.g. Santos et al. 1990): This effective stress is sometimes referred in literature to as the generalized effective stress or the Bishop's effective stress (e.g. Fjaer et al. 2008).
The effective stress coefficient, which appears in equation (56) is calculated using equation (5) as follows:  (55) to derive the following expression for the effective shear modulus of partially saturated rock: This equation (58) is more understandable than equation (55) from the physical point of view; however, it is still difficult to see how the capillary pressure is changing with the effective stress coefficient while all other parameters are kept constant. In the recent publication of Rozhko (2016), it was demonstrated that the effective stress coefficient, controlling the volume of partially-saturated crack can be well approximated by the wetting-phase saturation: Thus, the equation (58) can be simplified further as Figure 9 show collapse of data curves of Fig. 8 after scaling of vertical and horizontal axes (all curves follow the same trend). See text for details.
In this form, all input parameters of equation (60) have quite clear physical interpretation. Equation (60) shows no dependence on the fluid bulk moduli. Also, the critical stress amplitude ( τ c ) neither depends on the fluid bulk moduli, according to equation (54). Thus, the shear softening effect does not depend on the fluid bulk moduli. We discussed previously after equation (21) that the current model is applicable both for drained and undrained boundary conditions for both fluids, where the drained boundary conditions can be set by using zero values for fluid bulk moduli. Because equations (54), (55), (60) does not depend on the fluid bulk moduli, it implies that the shear softening does not depend on the fluid boundary conditions during pure shear stress-perturbation experiments, that is the same shear softening effect will be observed for drained and undrained boundary conditions. Let us now apply equations (54) and (55) to explain the diversity of experimental data, shown in Figure 8. We will modify the vertical and horizonal axis of Fig. 8 to demonstrate the 'data collapse'. The horizontal axis of Fig. 8(a,b) will be scaled by the critical wave amplitude, calculated with equation (54), that is τ → τ τ c . The vertical axis of Fig. 8(a) will be scaled by , whereas the vertical axis of Fig. 8 Figure 9 shows calculated results, shown in Fig. 8 after scaling of vertical and horizontal axes. Figure 9 shows the ''data collapse' when all curves follows the same trend. The data collapse curves do not coincide exactly because the non-linear effects and iterations are neglected during derivation of equations (54) and (55). However, these equations capture the most important physics and neglect input parameters, which are not important (i.e. bulk moduli of pore fluids). Equations (55) and ( (60)) show that the shear softening effect is strong in rock with compliant pores in which the capillary pressure is highly sensitive to the change of the confining stress, while if the rock material does not contain compliant pores or cracks the partial derivative, ∂p cap ∂σ will be small in this rock, and thus the shear softening effect will be small as well.

A P P L I C A T I O N O F T H E M O D E L T O R E S E R V O I R R O C K S
In this section, we discuss application of the model to reservoir rocks. The reservoir rocks may contain many cracks of different length and aspect ratio (e.g. , Zimmerman 1990;Anders, Laubach and Scholz 2014). Hooker et al. (2009) found that the aperture distribution of cracks between grains is best described by log-normal distribution, suggesting that a lower cut-off is attributed to the grain scale. Pruess and Tsang (1990) argued that for the log-normal aperture distribution, a simple approximation to capillary pressure of microcrack (or fracture) system is obtained in closed form that resemble the typical shape of Leverett's J-function, which is described in reservoir engineering as follows (e.g. Barenblatt et al. 1990): where J (S we ) is the dimensionless Leverett's J-function and κ c is the permeability of microcrack system and n c is the porosity of microcrack system. The Leverett's J-function is a generalized approach, valid for capillary pressure description both in consolidated and unconsolidated rocks and for description of capillary pressure in fractures and micro-cracks systems. It must be noted that we are focusing only on microcrack porosity and microcrack permeability in a dual-porosity and dualpermeability rocks. Such differentiation of pore system on compliant pores (cracks or fractures) and stiff pores (matrix) is quite common both in geophysics and reservoir engineering (e.g. Pride and Berryman 2003;Gerke and van Genuchten 1993). There Leverett J-function can be approximated in its simplest form by the power-law function (Brooks and Corey 1966): where λ > 0 and s 0 > 0 are empirical dimensionless coefficients, λ is the crack aperture distribution index and s 0 is related to the crack capillary entry pressure, as follows: Here for simplicity we consider the following relationships between permeability and porosity of the crack system: This relationship follows from the assumption that the crack permeability is proportional to the cube of crack aperture, whereas the crack porosity is linearly proportion to the crack aperture. Due to tortuosity effects, the cubic law scaling relationships between crack porosity and crack permeability may not be always applicable, in this case a different power-law scaling exponent needs to be used, that is κ c = κ 0 n q c , with q > 3. The porosity (and permeability) of the crack system is highly sensitive to the effective stress, which can be described by crack compressibility (e.g . Zimmerman 1990): where C pc is the crack compressibility under constant fluid pressure (e.g. Zimmerman 1990). In the simplest case, if C pc = const, equation (65) can be integrated, yielding: n c = n 0 exp(C pc (σ + p f l )) (e.g. Zimmerman 1990). It must be noted here again that we are using a sign convention when compressive stresses are negative, whereas in engineering and geosciences it is more common to use the sign convention when compressive stresses are positive. Derived equations (54) and (60) (and equation (65)) require that the compressive stress is negative, otherwise these equations must be modified by ∂ ∂σ → − ∂ ∂σ . Walsh (1965) derived the following equation for the crack: Note that according to equation (4), the crack compressibility is reciprocal to crack closure stress, whereas, according to equation (53) it can be related to the crack porosity, dry shear modulus and shear modulus of rock mineral. In equation (66), m is the aspect ratio of the crack, defined as: Here W is the crack aperture at given effective stress and L is the crack length. Experimental data, outlined in Zimmerman (1990, p. 133) shows that the typical values for the aspect ratio distribution function of the sandstone are in the range 10 −5 ࣠ m ࣠ 10 −3 with the most frequent value around m ∼ 10 −4 . Similar range for the crack aspect ratio scaling was suggested for carbonate rocks and granite by (Cheng and Toksöz 1979).
Next, we substitute equations (63)-(66), describing capillary pressure, permeability, porosity and compressibility of microcrack system in natural rocks, to equations (54) and (60), describing the critical wave amplitude ( τ c ) and the shear modulus of partially saturated rock (μ p.sat ). After simplification, we derived the following expressions for τ c and μ p.sat : and where the Young's angle (θ Y ) is calculated using equation (1), which is approximately equal to θ Y ≈ θ a +θ r 2 . The contact angles and interfacial tension can be measured from standard laboratory tests, whereas it extremely difficult to measure other input parameters, such as crack aspect ratio, crack porosity, crack permeability and crack aperture distribution. It is a general problem of any dual-porosity and dual-permeability models, as it is extremely difficult to constrain these input parameters, while at the same time these models can provide physical insights when single-porosity and single-permeability models are not able to.
Next, we consider saturation limits in equation (69). At zero saturation limit (S we → 0), we have as expected μ p.sat → μ dry if λ > 1 and μ p.sat → 0 if λ < 1. In granular material, the typical values for the pore size distribution index λ are in the range 1 < λ <4.2 (Brooks and Corey 1966;Laliberte, Corey and Brooks 2007), that is all reported values for λ index are higher than 1. However, those reported values include both pore systems: matrix and cracks, thus one may argue that obtained bounds for λ may not be representative for cracks. It is not difficult to demonstrate that there is a theoretical reason for λ > 1. The value λ < 1 would imply that the crack porosity is dominated by cracks with smallest possible aperture that is not possible for consolidated rocks, because it requires infinitely large crack area if the crack porosity is finite. In opposite, values λ > 1 imply that the crack porosity is dominated by cracks with highest possible aperture, that is more representative for consolidated rocks. Thus, in the limit S we → 0, equation (69) predicts the expected result: μ p.sat → μ dry . The case of full saturation must be discussed with more details. When S we → 1, equation (69) predicts μ p.sat < μ dry ; however, in the case of full saturation, we would also expect μ p.sat = μ dry , according to discussions around Fig. 1(a). Indeed, in the case of full saturation, the capillary pressure is 0 and after substitution of p cap = 0 into equation (60), we will get the expected result μ p.sat = μ dry when S we = 1. The reason why equation (69) predicts μ p.sat < μ dry when S we = 1 is due to Brooks and Corey (1966) approximation for the capillary pressure. Due to this approximation, the capillary pressure is not 0 in the case of full saturation, according to equation (62). The capillary pressure of rocks drops rapidly to 0 when S we → 1. Thus, Brooks and Corey (1966) approximation is not applicable to explain this endtail behaviour of the capillary pressure at S we → 1. Thus, according to equation (60), we would expect μ p.sat → μ dry in the same manner as p cap → 0, when S we → 1. Thus, if the rock is literally 100% saturated with a single fluid phase, there would not be any liquid bridges and thus no shear softening effect, according to our model. Several authors, who reported the shear softening effect for saturated rocks argued that it was extremely difficult to achieve full saturation of rock sample, especially if the rock sample is tight (e.g. Murphy 1984;Verwer et al. 2010;Li et al. 2017). Even after applying different advanced saturation techniques, there might be some pores that are not completely saturated, containing entrapped oil or gas bubbles. For example, Li et al. (2017) argued that the maximum wetting-phase saturation which was possible to achieve during experiments was around 98%, thus Brooks and Corey approximation is relevant for this case. Furthermore, published literature outlined in Introduction showed that the shear softening effect was not observed on all tested samples of the same rock. There might be several explanations for why some samples do not reveal the shear softening behaviour. One of the possible explanations is that those samples were fully (100%) saturated, thus capillary pressure in those samples were 0. At the reservoir conditions, when the fluid pressure is above bubble point, we do not expect to have any free gas in the formation. However, in the oil reservoir, we still have the interface menisci between oil and water which reduce the shear modulus. Thus, pore-scale interface phenomena effects are relevant both at the laboratory conditions and at the reservoir conditions.

D I S C U S S I O N S
In this section, we would like to discuss the important application of this theory: effect of pore fluids and saturation on mechanical strength. According to linear elastic fracture mechanics, the mechanical stiffness (shear modulus) and strength (tensile or uniaxial compressional) are interrelated properties (e.g. Paterson and Wong 2005). The tensile rock strength and uniaxial compression strength are linearly proportional to the fracture toughness (K Ic ), whereas the fracture toughness is related to the shear modulus (μ), Poisson's ratio and specific fracture energy (G c ) required to brake molecular bonds during propagation of fracture, as follows: Thus, the change of shear modulus by pore fluids will affect the mechanical strength of rocks. In equation (70), we consider specific fracture energy (G c ) as an intrinsic material property, independent on the fluid content. However, there are some publications where authors assume the change of apparent fracture energy G c by pore fluids to explain the water weakening effect (e.g. Bergsaker et al. 2016). In this paper, we assume that G c is independent on the fluid content, while the weakening is due to apparent change of the shear modulus. Furthermore, because bulk modulus is highly sensitive to pore fluids (Rozhko 2019), the Poisson's ratio in equation (70) is also affected by pore fluids. However, the fracture toughness in equation (70) is much more sensitive to the change of shear modulus, compared to Poisson's ratio, thus at the first-order approximation we can assume ν and G c to be independent on fluids in equation (70). Thus, using equations (69) and (70), we can estimate the 'fluid weakening' effect as K Ic p.sat   Table 2 shows published data by different authors for contact angles of water/air, ethylene-glycol/air and methanol/air interfaces on the calcite (Okayama, Keller and Luner 1997;Ethington 1990). Reported contact angles are quite different for the same substrate mineral and the same fluids. It will of course affect predictions by our model significantly. Let us consider the example of input parameters, given in Table 3. Those parameters are typical for chalk and taken from different publications given in Table 3. We assumed that the crack aspect ratio, reported for the limestone and dolomite by Cheng and Toksöz (1979) is also representative for chalk. The mercury injection capillary entry pressure into cracks was assumed to be similar to entry pressure, measured on chalk sample reported by Fabricius (2007), whereas the crack aperture distribution index for microcracks in chalk was assumed to be similar to the pore size distribution index for chalk, reported by Christoffersen and Whitson (1995). Figure 10 shows (a) calculated shear modulus, (b) fracture toughness and (c) capillary pressure of chalk, partially saturated with different liquid and gas system. As it was discussed in the pre-  (2007) Crack aspect ratio 0.5 × 10 −4 Cheng and Toksöz (1979) Crack aperture distribution index (λ) 3.4 Christoffersen and Whitson (1995) vious section, due to Brooks and Corey (1966) approximation, the capillary pressure is not vanishing to zero at full saturation. Thus, after substitution of p cap = 0 into equation (60), we will get the expected results μ p.sat = μ dry and K Ic p.sat = K Icdry . Figure 10(a,b) shows that shear modulus and fracture toughness are highly sensitive both to saturation, surface tension and contact angle measurement. The reduction of compressive strength by the degree of water saturation is a well-known effect in rock mechanics (Papamichos, Brignoli and Santarelli 1997;Fjaer et al. 2008). Figure 10(b) can be compared with experimental data for compressive strength of chalk, reported by Risnes et al. (2005), who demonstrated that methanolsaturated chalk sample is stronger than water-saturated chalk sample, but weaker that Ethylene-Glycol saturated chalk sample, while dry rock sample were the strongest. Dashed curves of Fig. 10(b) show the same trend, where the dashed curves are based on contact angle measurements of corresponding liquids on calcite samples from Mexico (Ethington 1990). Using equation (68), we can also estimate the critical wave amplitude when the shear softening is expected. Thus, for the water saturated sample and the contact angle, reported by Okayama et al. (1997), we have τ c = 0. This result indicates that for the complete wetting surface the slippage condition (and shear softening) will not depend on the wave amplitude. Results would be quite different if we would use contact angles reported by Ethington (1990). Ethington (1990) reported averaged advancing and receding contact angles for water air interface on calcite as θ a ≈ 47 0 and θ r ≈ 7 0 , those angles are quite different from Sessile drop contact angle θ Y ≈ 72 0 ,while the reason for this difference was not discussed by Ethington (1990). Thus, if we would use advancing and receding angles in equations (1) and (68) we will get τ c = 0.852 MPa. The stress amplitude of 0.852 MPa is Figure 10 Effect of pore liquids (water, ethylene glycol and methanol) and saturation on (a) shear modulus; (b) fracture toughness and (c) capillary pressure of chalk (calculated results). Horizontal axis in (a), (b) and (c) shows liquid saturation of microcracks in liquid-air system. Continuous and dashed curves correspond to contact angles measurements on calcite by Okayama et al. (1997) and Ethington (1990), respectively. See Tables 2 and 3  much larger than stress amplitude typical for seismic waves (10 2 − 10 4 Pa). Thus, the shear softening at seismic wave amplitudes may not always be observed, as it is highly sensitive towards wettability of the rock.

C O N C L U S I O N S
r The current paper suggests explanation to the effect of pore fluid on quasi-static shear modulus of isotropic partially saturated rock. It is demonstrated that the shear softening is related to the pore-scale interfacial phenomena effects, specifically, to the slippage and rolling of fluid molecules over the surface of the solid at the contact line location during the contact line motion. Classical Biot-Gassmann's models assume Stoke's no-slip boundary conditions between fluid and solid, while due to slippage at the solidfluid interface the shear modulus of partially saturated rock can be lower than shear modulus of dry rocks. The slippage condition depends not only on material and fluid properties, but also on the wave amplitude. Only for the complete wetting surface, the slippage condition will not depend on the wave amplitude. It is demonstrated that the quasi-static shear modulus of partially saturated rock is sensitive to the following fluid properties: interfacial tension between immiscible fluids, fluid saturation and rock wettability (contact angles) towards saturating fluids and independent on bulk moduli of pore fluids.
r Developed closed form equation allowing to estimate the shear modulus of a partially saturated rock and considering the following additional parameters: capillary pressure, stress sensitivity of capillary pressure and saturation sensitivity of capillary pressure. The shear softening is demonstrated to be possible only in rocks with dual porosity and dual permeability. Also, we predicted the critical wave amplitude defining the onset of shear softening. If the wave amplitude is lower than the critical amplitude, the shear modulus of partially saturated rock will be equal to shear modulus of dry rock, whereas if the wave amplitude is higher than its critical value, then the shear softening effect will take place. This derived equation depends on advancing and receding contact angles, wettability sensitivity of capillary pressure and stress sensitivity of capillary pressure.
r Suggested model was applied to estimate the shear softening and shear weakening effects of pore fluids in natural rocks with log-normal distribution of micro-crack aperture. In such rocks, the capillary pressure of microcrack system was described by Leverett's J-function with Brook-Corey's saturation dependency. Derived closed-form analytical solutions explain alteration of shear modulus and strength by pore fluids and saturation.

A C K N O W L E D G E M E N T S
I would like to thank Mahyar Madadi and anonymous reviewer for comments and suggestions improving the manuscript.