Dynamic Relative Permeabilities for Partially Saturated Porous Media Accounting for Viscous Coupling Effects: An Analytical Solution

We present an analytical model to compute frequency-dependent relative permeability functions for partially saturated porous media accounting for viscous coupling effects. For this, we consider the oscillatory motion of two immiscible fluid phases and solve the Navier–Stokes equations at the pore scale using suitable interface conditions between fluids. These calculations are combined with the generalized two-phase flow Darcy equations to obtain the corresponding upscaled macroscopic fluxes. By means of an analog pore model consisting of a bundle of cylindrical capillaries in which pore fluids are distributed in a concentric manner, we find closed analytical expressions for the complex-valued and frequency- and saturation-dependent relative permeability functions. These expressions allow for a direct assessment of viscous coupling effects on oscillatory flow for all frequencies and saturations. Our results show that viscous coupling effects significantly affect flow characteristics in the viscous and inertial regimes. Dynamic relative permeabilities are affected by the pore fluid densities and viscosities. Moreover, viscous coupling effects may induce two critical frequencies in the dynamic relative permeability curves, a characteristic that cannot be addressed by extending the classic dynamic permeability definition to partially saturated scenarios using effective fluids. The theoretical derivations and results presented in this work have implications for the estimation and interpretation of seismic and seismoelectric responses of partially saturated porous media.


Introduction
Darcy's law is extensively used by geophysicists, hydrologists, geologists, and reservoir engineers to describe fluid flow in monosaturated porous media at the macroscopic scale. In this context, the permeability is a fundamental hydraulic material property that relates the fluid pressure gradient to the flow rate per unit area (e.g., Bear 1972). Given that many practical problems involve the simultaneous flow of two immiscible pore fluid phases, such as oil and water, or gas and water, Darcy's law has been extended to account for the particular characteristics of this scenario (e.g., Bear 1972).
In the presence of immiscible two-phase flow in porous media, the fluid phases flow through different regions within the pore structure, and thus, each phase effectively senses a different permeability. In this context, relative permeabilities arise as an essential concept for understanding two-phase flow characteristics in porous media, providing valuable information regarding the permeability reduction sensed by a given fluid phase due to the presence of a second phase in the pore space (e.g., Bear 1972;Blunt 2017). In general, relative permeabilities are considered to be (i) functions of the saturation state and (ii) independent of the physical properties of the pore fluids, such as viscosities or densities (Burdine 1953;Mualem 1976). However, it has been classically argued by many authors that the physical properties of flowing pore fluid phases have non-negligible effects on the relative permeability functions (e.g., Yuster 1951;Scott and Rose 1953;Bolt and Groenevelt 1969;de Gennes 1983;Bachmat and Bear 1986;de la Cruz and Spanos 1983;Whitaker 1986;Kalaydjian 1987Kalaydjian , 1990. A particularly interesting phenomenon, which is associated with deviations from the classical extension of Darcy's law to multiphase flow, is the so-called viscous coupling effect. The classic immiscible two-phase flow formulation of Darcy's equations commonly assumes that each fluid phase, wetting or non-wetting, is driven by a unique pressure gradient acting across the corresponding phase (Peaceman 1977;Bear 1972). However, this conceptualization does not account for the effects of momentum transfer across interfaces separating the immiscible pore fluids. To account for the fact that, for instance, a pressure gradient acting on one phase can induce flow of the other phase, and vice versa, the classic two-phase formulation of Darcy's equations needs to be modified. Theoretical developments, based on volume averaging techniques and the principles of thermodynamics, have been used to upscale the Stokes equations in porous media and, thus, to derive novel and more consistent macroscopic two-phase flow formulations (e.g., de Gennes 1983;Bachmat and Bear 1986;de la Cruz and Spanos 1983;Whitaker 1986;Kalaydjian 1987Kalaydjian , 1990. The resulting mathematical expressions, to which we refer in the following as the generalized two-phase flow equations, do not only include the "classic" relative permeability terms but, also, exhibit "cross-terms," which account for interfacial coupling, which include viscous and capillary effects (Ayub and Bentsen 1999). In order to gain insights into the particular effects that viscous coupling has on the relative permeabilities, a number of researchers (e.g., Yuster 1951;Bacri et al. 1990;Rose 1990;Ehrlich 1993) employed simple analogous models, composed of partially saturated capillaries or slates, to obtain analytical solutions to the problem. Later on, interfacial coupling effects on the relative permeabilities were analyzed by means of laboratory experiments (e.g Avraam and Payatakes 1995;Bentsen 1998;Avraam and Payatakes 1999;Roman et al. 2020) and numerical models (e.g, Rothman 1990; Gunstensen and Rothman 1993;Rakotomalala et al. 1995;Huang and Lu 2009;Li et al. 2005). The corresponding results show that such effects produce non-negligible changes in the relative permeability curves of partially saturated porous media and, thus, that disregarding them may lead to interpretation errors.
To our knowledge, all of the above-mentioned studies were performed under steady-state and viscous-dominated flow conditions. However, in many pertinent scenarios in geophysics, such as flow generated by a passing seismic wave, pressure gradients and pore fluid velocities are of oscillatory nature and present inertial effects (e.g., Johnson et al. 1987;Zhou and Sheng 1989). To date, the effects of viscous coupling on the relative permeability curves under oscillatory flow conditions accounting for inertial effects remain largely unexplored.
When studying oscillatory flow scenarios, the viscous skin depth, that is, the spatial extension within the fluid with respect to the pore walls where viscous effects dominate over inertial effects, decreases with the angular frequency . For sufficiently high frequencies, the viscous skin depth becomes smaller than the pore size, and thus, inertia affects the flow characteristics. As such, Darcy's law in its classical form is no longer valid and several researchers have proposed to define a complex-valued and frequency-dependent permeability (e.g., Johnson et al. 1987;Zhou and Sheng 1989). This parameter, usually referred to as dynamic permeability ( ) , correctly accounts for (i) the out-of-phase motion of the fluid with respect to the pressure forcing and (ii) the flux reduction occurring at sufficiently high frequencies due to inertial effects. The dynamic permeability concept is of fundamental importance in the context of seismic (e.g., Biot 1962;Pride 2005) and seismoelectric analyses (e.g., Pride 1994;Jardani and Revil 2015). Even though the behavior of ( ) has been extensively explored for monosaturated porous media (e.g., Johnson et al. 1987;Sheng and Zhou 1988;Charlaix et al. 1988;Zhou and Sheng 1989;Chapman and Higdon 1992;Smeulders et al. 1992;Pride et al. 1993;Müller and Sahay 2011;Huber and Su 2015;Li et al. 2021), the corresponding behavior under partially saturated conditions has not yet been rigorously explored. Revil and Mahardika (2013) proposed an empirical model to account for a saturation dependence in the effective dynamic permeability of the wetting phase in the context of seismoelectric studies. More recently, Solazzi et al. (2020) have proposed an approach to compute the effective dynamic permeability functions of partially saturated media using bundle-of-tube models and the concept of capillary pressure equilibrium. However, both of these works disregard fluid-fluid interactions. As far as the authors know, Auriault et al. (1989) is the only work that accounts for dynamic permeability estimates of partially saturated media using the generalized twophase Darcy equations. The employed homogenization procedure is, however, highly complex and the expressions for the dynamic permeabilities are not explicitly given. Consequently, more work is needed to provide explicit simplified solutions for the dynamic permeability of partially saturated media accounting for interfacial coupling in general and for viscous coupling in particular.
In this work, we explore the effects of viscous coupling on dynamic relative permeability functions of partially saturated porous media under oscillatory flow conditions. To this end, we consider a simple analogous model, inspired by the work of Yuster (1951), in which concentric flow is assumed in cylindrical capillaries. We derive closed analytical expressions for the dynamic relative permeabilities, which are valid across the entire frequency range. Finally, we analyze the behavior of the corresponding functions for different viscosities, densities, frequencies, and saturations.

3 2 Theoretical Background
In this section, we briefly summarize the (i) classic and (ii) generalized/extended Darcy equations for immiscible two-phase flow in porous media for viscous-dominated flow. Based on Darcy's equation for single-phase flow, we introduce the concept of dynamic permeability. The latter will be extended to partially saturated scenarios in Sect. 3 of this paper.

Classic Two-Phase Flow Equations in Porous Media
The traditional two-phase flow formulation of Darcy's law is given by (Bear 1972) where v i denotes the Darcy velocity for phase i, being either the wetting w or non-wetting fluid phase n, and p i the fluid pressure. The mobilities i respond to where i is the dynamic viscosity of the fluid, the intrinsic permeability determined by the structure of the porous medium alone, and r i the relative permeability of phase i. The latter is commonly related to the wetting fluid saturation S w and to the pore space characteristics (e.g., Mualem 1976;Burdine 1953).
Recall that the Darcy velocity for phase i, that is, v i , is a volume average of the corresponding microscopic fluid velocity w i relative to the pore walls (e.g., Auriault et al. 1985) where V t is the total volume of a representative elementary volume (REV) of the porous medium of interest, Ω p is the dominium associated with the pore space, is the porosity, and ⟨⋅⟩ denotes the average over the pore space of the REV. Note that w i is zero for the Ω p -regions in which the i-phase is absent. As previously mentioned, if one wishes to account for the effects of viscous coupling in the two-phase flow equations, Eq. (1) needs to be modified.

Generalized Two-Phase Flow Equations in Porous Media
Given that, at the pore scale, flowing immiscible fluids may interact across their respective interfaces (Fig. 1), several theoretical approaches have been proposed to account for fluid-fluid momentum transfer in multiphase flow in porous media (e.g de Gennes 1983; Bachmat and Bear 1986;de la Cruz and Spanos 1983;Whitaker 1986;Kalaydjian 1987Kalaydjian , 1990. The generalized two-phase flow equations in porous media respond to (e.g., Dullien and Dong 1996) where ij are the generalized mobilities, which are given by with r j i denoting the generalized relative permeability coefficients. The off-diagonal coefficients r n w and r w n account for interfacial coupling between the fluid phases. Please note that, in the context of Eqs. (4) and (5), the indexes i and j in the relative permeabilities r j i are to be interpreted as follows: (i) the subindex i makes reference to the flowing phase that is affected by the corresponding relative permeability and (ii) the subindex j indicates the pressure gradient ( ∇p j ) which is modulated by the corresponding term. For instance, r n w modulates the effect of ∇p n on v w . We remark that Eq. (3) is also behind the definition of v w and v n in this formulation.
As previously stated, the effects of viscous coupling on two-phase flow through porous media have been analyzed through laboratory (e.g Payatakes 1995, 1999;Roman et al. 2016Roman et al. , 2020 and numerical (e.g, Gunstensen and Rothman 1993;Li et al. 2005) experiments. In general, it is observed that an increment in the viscosity ratio M = n ∕ w leads to a significant increase in r n w and r n n . The non-wetting phase experiences an apparent hydraulic slip, usually referred to as lubricating slip, due to wetting films located near the pore walls, which increases with M. Conversely, the terms r w n and r w w are virtually insensitive to variations in the viscosity ratio M. Interestingly, some works indicate that the cross-term mobilities should be identical r w n ∕ w ≡ r n w ∕ n , which is known as the Onsager reciprocity (e.g., Gunstensen and Rothman 1993). It is important to mention that crossed mobility terms may show different values when irreversible processes take place, such as mobilization of isolated fluid ganglia or bubbles (e.g., Li et al. 2005;Payatakes 1995, 1999). The purpose of this paper is to explore the effects of viscous coupling extending Eqs. (4) and (5) to the entire frequency range in the context of oscillatory flow. For this, we will make use of the concept of the frequency-dependent dynamic permeability ( ).
Non-wetting phase Wetting phase Solid matrix

Dynamic Permeability
In most studies on Darcy flow, the permeability is considered to be real-valued and frequency-independent. This is indeed valid as long as the fluid flow through the pore space is viscosity-dominated, that is, when the viscous skin depth is greater than the characteristic pore size (Johnson et al. 1987). However, in the presence of oscillatory pressure forcing, inertial effects increase with frequency and viscous boundary layers may develop along the pore walls. This physical process is evidenced as a decrease in the flow rate and a phase shift between the hydraulic pressure variations and the fluid flow. In this context, the classic formulation of Darcy's law is no longer valid. To circumvent this issue, several investigators have proposed to generalize the classic monosaturated Darcy's law through the introduction of a complex-valued and frequencydependent permeability (e.g., Johnson et al. 1987;Sheng and Zhou 1988;Zhou and Sheng 1989) where = 2 f is the angular frequency and f denotes the linear frequency.
The behavior of ( ) with frequency is controlled by both viscous and inertial forces. The transition from viscous-to inertia-dominated flow occurs at the so-called critical angular frequency c , for which the characteristic pore radius R p of the medium becomes comparable to the viscous skin depth (e.g., Johnson et al. 1987) with = ∕ denoting the kinematic viscosity and the fluid density (Johnson et al. 1987). Please note that the dynamic permeability ( ) does not only depend on the properties of the porous medium alone, but also depend on the physical properties of the saturating fluid ( and ), and on the angular frequency .
In the following, we employ the concept of the dynamic permeability to extend Eqs. (4) and (5) to the entire frequency range. As outlined below, we do so by including all saturation-and frequency-dependent effects in new complex-valued and frequencydependent relative permeability functions r j i ( , S w ) . For the latter to be valid, the derived r j i ( , S w ) functions need to converge to their real-valued counterparts for sufficiently low frequencies.

Model Development
In this section, we derive closed analytical expressions for the dynamic relative permeability functions of partially saturated porous media accounting for viscous coupling effects. For this purpose, we use a simple pore-scale analog model and solve Navier-Stokes equations under oscillatory flow conditions at the microscopic scale, accounting for suitable boundary and interface conditions. Finally, we use the generalized two-phase flow Darcy equations to obtain the extensions of the relative permeability functions to the entire frequency range.

Pore-Scale Characterization
In the presence of a fluid pressure gradient, pore fluids follow preferential flow paths through the pore space. These flow channels are hereby conceptualized using a capillary tube geometry. Even though this pore-scale geometry constitutes an obvious oversimplification, models based on such an approach have a long history (Kozeny 1927) and have proved to be highly effective for the macroscopic description of the hydraulic characteristics of porous media (e.g., Guarracino 2007;Nguyen et al. 2021). In this work, we consider a porous structure composed of capillary tubes of a single characteristic radius R saturated concentrically with two different fluid phases (Fig. 2a). This setup was previously proposed by Yuster (1951) and Bacri et al. (1990) to study the effects of viscous coupling in steady-state Poiseuille-type flow. The non-wetting fluid saturates the center of the capillary ( 0 ≤ r ≤ R 1 ), while the wetting phase saturates the annular space between the non-wetting phase and the pore wall ( R 1 < r ≤ R ) (Fig. 2b), with r being the radial coordinate. The proposed mechanistic description of oscillatory fluid flow processes is then upscaled to the REV scale using the generalized Darcy's law. We therefore wish to extend the steady-state solutions derived by Yuster (1951) and Bacri et al. (1990), which are summarized in "Appendix 1," to the entire frequency range.
In the context of the model described in Fig. 2, the expressions for the non-wetting S n and wetting S w phase saturations, the porosity , and the permeability respond to respectively. In Eqs. (11) and (12), L is the sidelength of the cubic REV, = L p ∕L is the tortuosity, with L p > L denoting the length of the flow channels, and N p is the number of capillaries in the REV.

Oscillatory Flow in a Partially Saturated Capillary
Let us consider a homogeneous, isotropic, porous solid, saturated with a Newtonian fluid. Any process that perturbs the pore fluid pressure, such as, the passage of a seismic wave, will generate fluid displacements relative to the pore walls. Hereafter, without loss of generality for the considered problem, we assume that pore walls are still. The microscopic motion of the fluid phase within the medium is described by the Navier-Stokes equation (e.g., Landau and Lifshitz 1959) where w is the fluid velocity at the microscopic scale, p the fluid pressure, t the time, and the secondary viscosity coefficient. Here, f denotes the surface tension force, which, assuming constant interfacial tension , responds to (e.g., Tryggvason et al. 2011) where n is the unit normal vector with respect to the interface, is the mean curvature of the interface, and S is a delta function that takes a unit value on the interface between fluids and is null elsewhere.
Following Johnson et al. (1987), we shall assume that the solid matrix is rigid and that the pore space is saturated with an incompressible ( ∇ ⋅ w = 0 ) fluid. Please note that the fluid can indeed be regarded as incompressible at the pore scale provided that the prevailing acoustic wavelengths in the fluid are much larger than the typical pore size (Johnson et al. 1987;Charlaix et al. 1988;Zhou and Sheng 1989). We also assume that the fluid flow within the pore space is characterized by a small Reynold's number, such that laminar flow prevails (Johnson et al. 1987;Smeulders et al. 1992). Then, Eq. (13) reduces to Considering a cylindrical capillary whose axis is oriented in the direction z and assuming that (i) the radial and angular components of the fluid velocity are null, that is, w = [0, 0, w(r)] and that (ii) capillary forces act radially, in agreement with the concentric flow model, Eq. (15) can be expressed in cylindrical coordinates as (Landau and Lifshitz 1959) On the one hand, Eq. (16) expresses the continuity of the normal stress at the fluid-fluid interface which, considering the Young-Laplace equation for a cylindrical interface, is satisfied if the capillary force is counteracted by a fluid pressure jump at r = R 1 , that is, Coward et al. 1995). Note that this condition implies the equality ∇p w = ∇p n must be fulfilled for concentric flow to occur (e.g., Picchi and Battiato 2018). On the other hand, Eq (17) expresses the momentum conservation of the system in the axial direction. Let us define the external forcing exerted to the pore fluid in the z direction as and consider that both the relative fluid velocity and the forcing term present an oscillatory form given by w =ŵ e −i t and X =X e −i t , respectively, with i denoting the imaginary unit. Then, Eq. (17) can be expressed in the space frequency domain as where k = √ i ∕ can be conceptualized as the wave number of the oscillatory process (Johnson et al. 1987). By taking w(r) =ŵ −X k 2 , the corresponding equation reduces to a Bessel-type differential equation Please note that this equation holds for two immiscible fluid phases flowing concentrically in a cylindrical capillary (Fig. 2) and is the key expression to be solved in this work. In the context of the proposed model, the relative fluid velocities with respect to the pore walls respond to The general solution for the Bessel equation (Eq. 20) is a combination of zeroth-order Bessel functions J 0 and Y 0 of the first and second kinds, respectively (Abramowitz and Stegun 1965) where ℭ j 1 and ℭ j 2 , with j = n, w , are four integration constants that can be computed from the boundary conditions of the problem, which are detailed in the next section.

Boundary Conditions of the Problem
In order to obtain the relative velocities of both fluid phases (Eqs. 22 and 23) and, through them, derive the expressions for the Darcy velocities (Eq. 3) and dynamic relative permeabilities r j i (Eqs. 4 and 5), we need a set of suitable boundary conditions. To this end, we consider the following: (i) a non-diverging relative velocity at the center of the capillary; (ii) a no-slip condition at the pore walls; continuity of (iii) the shear stresses and of (iv) the velocity at the interface between fluids.
The first assumption implies that the relative velocity ŵ n value at r = 0 must be finite. Given that second-kind Bessel functions present a singularity at the origin, that is, lim r→0 Y 0 (k n r) = ∞ , the ℭ n 2 term has to be null for ŵ n to be non-diverging in Eq. (22). Conversely, the no-slip assumption at the pore walls implies that ŵ w | (r=R) = 0 , which, considering Eq. (22), reduces to On the other hand, we have the continuity of the fluid shear stress j = j dŵ j ∕dr at the interface r = R 1 , which results in where J 1 and Y 1 denote the first-order Bessel functions of the first and second kinds, respectively.
Finally, at the interface between the immiscible fluids, we impose the continuity of fluid relative velocities for the phases ŵ w −ŵ n | (r=R1) = 0 , which, using Eqs. (22) and (23), yields Equations (24), (25), and (26) constitute a linear system of equations that can be solved analytically to obtain ℭ w 1 , ℭ w 2 , and ℭ n 1 . The inferred expressions were verified using the software Mathematica (Wolfram Research, Inc. 2021). Once these constants are obtained, we can compute the relative velocities of both fluid phases (Eqs. 22 and 23).

Analytical Solution of the Problem
The macroscopic Darcy velocities for the wetting and non-wetting phases can be obtained by means of Eq. (3), that is, by averaging the corresponding relative velocities ⟨ŵ j ⟩ over the cross section of the capillaries and multiplying the latter with the porosity Once again, these integrals can be solved analytically in the context of the proposed model. Comparing the corresponding results with Eqs. (4) and (5), we obtain the following expressions for the dynamic mobilities The expressions for the coefficients S , M j , and L j are summarized in "Appendix 2." In order to simplify the interpretation of the derivation and of the corresponding analytical solution, we detail the notation used in this work in "Appendix 3." The corresponding dynamic relative permeabilities are given by These closed analytical expressions are the essential result of this work. We present the equations and coefficients in terms of the characteristic radius R. However, in the context of our model R = √ 8 2 ∕ , and thus, the corresponding equations can alternatively be expressed in terms of κ and . It is important to remark here that, as evidenced in Eqs. (50) to (70), the kinematic viscosities of the fluids n = n ∕ n and w = w ∕ w play a key role in determining the dynamic relative permeabilities.

Dependence on Saturation
Relative permeabilities are commonly analyzed as functions of the wetting phase saturation (e.g., Bear 1972). Let us then analyze the behavior of the above-defined dynamic relative permeabilities r j i ( ) (Eqs. 32 to 35) as functions of S w for different frequencies f. For this, we consider a simple porous medium, as described by Fig. 1, characterized by = 1 and R = 100 μ m. We take water and oil as the saturating fluids, the former having a density w = 1000 kg m −3 , and the latter one of n = 700 kg m −3 , thus resulting in a density ratio of 1 3 N = 0.7 . In the following, the water viscosity is taken as w = 0.001 Pa s. For oil, we consider two viscosities: (i) n = 0.01 Pa s and (ii) n = 0.001 Pa s, which, in turn, correspond to viscosity ratios of M = 10 and M = 1 , respectively. These two scenarios allow us to discern the effects of the viscosity ratio M on the dynamic relative permeability curves. Figure 3 shows the real part of the relative permeabilities r j i ( ) as a function of S w for different frequencies. We consider f = 1 Hz (black lines), f = 100 Hz (blue lines), and f = 1000 Hz (red lines). The top row (Fig. 3a, b) corresponds to the relative permeabilities r w w and r n w , which modulate the Darcy velocity of water v w . The bottom row (Fig. 3c, d) corresponds to the relative permeabilities r w n and r n n , which modulate the Darcy velocity of oil v n . The left (Fig. 3a, c) and right (Fig. 3b, d) columns are associated with the relative permeabilities that modulate the effects of the pressure gradients of the wetting and non-wetting phases ∇p w and ∇p n , respectively, on the flow characteristics. Recall that we consider two scenarios for the oil viscosity: (i) M = 1 (dashed lines with square markers) and (ii) M = 10 (solid lines). Let us first analyze the response of the medium for f = 1 Hz (black lines). In this relatively low-frequency range, and for the considered medium and fluid properties, fluid flow lies in the viscous regime, and thus, it can be shown that the relative permeabilities replicate those predicted by Yuster (1951) and Bacri et al. (1990) ("Appendix 1"). We observe that r w w for different values of M (Fig. 3a) increases with S w and that the behavior is similar for different M. This shows that, for this frequency, r w w is not affected by viscous coupling The medium is characterized by a tortuosity and a pore radius of = 1 and R = 100 μ m, respectively effects. Conversely, r n n (Fig. 3d) monotonically decreases with S w for M = 1 (dashed black line) while it increases ( 0 < S w < 0.4 ), presenting values greater than unity, and then decreases ( 0.4 < S w < 1 ) for M = 10 (solid black line). This behavior is also expected, as the non-wetting phase experiences an apparent hydraulic slip, usually referred to as lubricating slip, due to the presence of a wetting film located near the pore walls (e.g., Li et al. 2005). The greater the value of M is, the greater is the effect of slip on the r n n . Also for f = 1 Hz (black lines), we note a large variation in the crossed term r n w (Fig. 3b) with increasing M (compare dashed and solid lines), indicating that, for M = 10 , the pressure gradient acting on the non-wetting phase increases its capability to drag the wetting phase into flow. Conversely, the remaining crossed term r w n is virtually insensitive to M for this relatively low frequency. Finally, we remark that even though, in the low-frequency range, the crossed terms present identical behaviors for M = 1 and their corresponding responses diverge for M ≠ 1 , the relation r w n ∕ w ≡ r n w ∕ n holds in all cases. Let us now analyze the response of ℜ{ r j i (S w , )} as a function of saturation for frequencies of f = 100 Hz and f = 1000 Hz, for which inertial effects start to prevail. We observe that ℜ{ r j i (S w , )} decreases with frequency, irrespective of the saturation (Fig. 3). Note that a similar behavior was reported by Solazzi et al. (2020), who analyzed the effective dynamic permeabilities disregarding viscous coupling effects. Figure 3 shows that, for M = 1 (dashed lines), the dynamic relative permeabilities do not differ significantly for f = 1 Hz (black lines) and f = 100 Hz (blue lines). However, for M = 10 (solid lines), the dynamic relative permeabilities present a greater change when the frequency increases from 1 Hz (black solid lines) to 100 Hz (blue solid lines). The reason for this is that we are increasing M by increasing the non-wetting phase viscosity. The latter shifts c toward lower frequencies (Eq. 8), and thus, inertial effects arise at lower frequencies. We thus observe the effects of both inertial and viscous coupling on ℜ{ r j i } for M = 10 , while we only observe viscous coupling, but not inertia, in the corresponding curves for M = 1 . At the same time, we note that ℜ{ r n w } (Fig. 3b) and ℜ{ r w n } (Fig. 3c) show negative values for f = 1000 Hz and M = 10 (red solid lines) for 0.6 < S w < 1 . This is a highly interesting characteristic, which is also present in the work of Auriault et al. (1989), that cannot be accounted for without a frequency-dependent formulation accounting for viscous coupling.
For completeness, we illustrate the imaginary part of r j i S w , (Fig. 4) as a function of saturation for frequencies of f = 1 Hz (black lines), f = 100 Hz (blue lines), and f = 1000 Hz (red lines). As previously shown in Fig. 3, we take M = 1 (dashed lines) and M = 10 (solid lines) by considering different non-wetting phase viscosities. We remark that the imaginary part in the dynamic permeability is associated with an out-of-phase motion between the pressure forcing and the resulting flux due to inertial effects (e.g., Johnson et al. 1987). As expected, ℑ{ r j i (S w , )} is virtually null for f = 1 Hz (black lines) for both M = 1 and M = 10 (Fig. 4). Also, we observe that, in the case of M = 1 (dashed lines), ℑ{ r j i } tends to increase with frequency. However, when M = 10 (solid lines), this is not the case. For example, Fig. 4a shows that, for M = 10 and saturations such that S w < 0.9 , the values of ℑ{ r w w } associated with f = 100 Hz (blue solid lines) are greater than those associated with f = 1000 Hz (red solid lines). Conversely, for S w > 0.9 the ℑ{ r w w } values for f = 1000 Hz (blue solid lines) are greater than those associated with f = 100 Hz (red solid lines). These effects are better appreciated analyzing the responses as functions of the frequency, which is done in the following.

Dependence on Frequency
The dynamic permeability in monosaturated porous media is commonly defined as a function of the frequency (e.g., Johnson et al. 1987;Sheng and Zhou 1988;Zhou and Sheng 1989). In the following, we therefore analyze the real and imaginary parts of dynamic relative permeability curves as functions of frequency (Figs. 5 and 6). For this, we consider the fluids and porous media defined in the previous subsection. We depict the behavior for three saturations S w of 0.3 (green lines), 0.6 (blue lines), and 0.9 (red lines). Figure 5 shows the real part of the dynamic relative permeabilities as functions of frequency for different wetting phase saturations. We consider the dynamic relative permeabilities for M = 1 (dashed lines) and M = 10 (solid lines). As expected, we note that ℜ{ r w w } (Fig. 5a) increases with S w while ℜ{ r n n } (Fig. 5d) decreases with S w , for both M = 1 (dashed lines) M = 10 (solid lines). As predicted by the model of Yuster (1951) (Eqs. 46 and 49), we note that changing M has no effects on the low-frequency behavior of ℜ{ r w w } and ℜ{ r w n } , which are the terms that modulate flow induced by the water pressure gradient ∇p w . Conversely, we evidence a large increase in the values of ℜ{ r n w } (Fig. 5b) and ℜ{ r n n } (Fig. 5d) with increasing M. These terms modulate flow by ∇p n . A general characteristic of all dynamic relative permeabilities is that the onset of inertial effects is evidenced by the inflection of the curves occurring at c . For increasing values of n and, thus, of M, c moves toward to lower frequencies. Figure 5a shows that, . The wetting fluid is water and the non-wetting fluid is oil, which results in a density ratio of N = 0.7 . We consider two different oil viscosities, which result in two viscosity ratios: M = 1 (dashed lines with square markers) and M = 10 (solid lines). The medium is characterized by a tortuosity and a pore radius of = 1 and R = 100 μ m, respectively when considering M = 10 , two inflection points arise, which are not present for M = 1 . This implies that the effects of viscous coupling might generate deviations from the seemingly universal scaling law of the monosaturated dynamic permeability, which was proposed in several classical works (Johnson et al. 1987;Sheng and Zhou 1988;Zhou and Sheng 1989). This is an important feature, not present in the effective dynamic permeability curves in partially saturated porous media analyzed by Solazzi et al. (2020), which disregarded viscous coupling effects. Finally, we observe, once again, that nondiagonal relative permeability terms may present small negative values in a restricted frequency range, located approximately between 80 and 200 Hz. Figure 6 shows the imaginary part of the dynamic relative permeabilities as functions of frequency for different wetting phase saturations. We observe that the imaginary part of the dynamic relative permeabilities for M = 1 (dashed lines) shows a peak where the corresponding real parts showed an inflection in Fig. 5. Comparing the imaginary part of the dynamic permeabilities for M = 1 (dashed lines) and M = 10 (solid lines), we see, for all dynamic permeabilities, that the peak value shifts to lower frequencies for increasing M. Interestingly, in Fig. 6a, when considering M = 10 , we observe two peaks, whereas only one is present when disregarding viscous coupling effects ( M = 1 ). This feature was evidently depicted in Fig. 5a as two inflection points. We also evidence a large increase in the values of ℑ{ r n w } (Fig. 6b) and ℑ{ r n n } (Fig. 6d) for M = 10 (solid Real part of the dynamic relative permeabilities as functions of frequency for different wetting phase saturations S w . We consider water and oil as saturating fluids, and saturations of 0.3 (green lines), 0.6 (blue lines), and 0.9 (red lines). We consider two different oil viscosities, which result in two viscosity ratios: M = 1 (dashed lines) and M = 10 (solid lines). The density ratio is taken as N = 0.7 . The medium is characterized by a tortuosity and a pore radius of = 1 and R = 100 μ m, respectively lines) compared to M = 1 (dashed lines), which is expected as the viscosity of the nonwetting phase is increased tenfold to increase M.

Discussion and Conclusions
We have presented an analytical model to compute viscous coupling effects on the relative permeability functions in the presence of oscillatory flow. For this, we have extended the dynamic permeability concept to porous media saturated by two immiscible fluid phases.
The proposed model provides saturation-and frequency-dependent relative permeability functions, which account for momentum transfer between fluids.
Our results show that the real parts of the dynamic relative permeabilities decrease with frequency, with inflection frequencies that depend on the fluid viscosities and densities. Also, non-diagonal relative permeability terms may present small negative values in a restricted frequency range, an observation also made by Auriault et al. (1989). These frequency-dependent effects are also visible in the imaginary parts of the dynamic relative permeabilities. In particular, the relative permeability terms r n w and r n n , which modulate the effects of the non-wetting fluid pressure gradient ∇p n on the flow, are very sensitive to the viscosity ratio, exhibiting variations with regard to the amplitudes and the critical On the other hand, the dynamic relative permeability terms r w w and r w n that modulate the effects of the wetting phase pressure gradient ∇p w on the flow may, on top of the latter effects, exhibit two characteristic frequencies due to viscous coupling effects for sufficiently high M values. When analyzed as a function of saturation S w , we observe that the dynamic relative permeabilities may result in complex curves for increasing frequencies, presenting values larger than unity for r n n and r n w for sufficiently large values of M. A number of uncertainties arise as a consequence of the assumptions made in our analytical model. One obvious simplification adopted in this work is the capillary bundle assumption, which is motivated by the possibility of obtaining closed analytical expressions for the dynamic relative permeabilities. As such, our model does not account for axial variations in pore cross sections associated with pore throats nor for multiple pore interconnections, which are known to produce deviations in monophasic dynamic permeability functions (e.g., Pride et al. 1993;Huber and Su 2015). On the other hand, the concentric flow assumption implies that a (i) continuous interconnection exists between fluid phases and (ii) that capillary tension acts only in the radial direction. In this sense, research shows that the fluid-fluid interfacial area as well as the presence of isolated fluid ganglia is important parameters to account for when computing the relative permeability functions (e.g, Rakotomalala et al. 1995;Li et al. 2005;Roman et al. 2020). When fluid ganglia are present, capillary forces play a predominant role in the flow of the fluid phases, opposing the non-wetting fluid displacement and invalidating the concentric-type flow conditions (e.g., Badchin and Yuan 1997;Ayub and Bentsen 1999;Bentsen 2001). Such capillary coupling effects cannot be accounted for in the proposed analytical model, which focuses on viscous coupling effects alone. Further research is certainly needed to account for both viscous and capillary coupling effects when estimating the dynamic relative permeability estimates. Following classic two-phase flow extensions of Biot's poroelasticity theory (e.g, Santos et al. 1990;Tuncay and Corapcioglu 1997;Lo et al. 2005;Jardani and Revil 2015), the proposed model disregards oscillatory saturation changes. However, we remark that this is an interesting process associated with drainage and imbibition characteristics, which merits further research.
Oscillatory flow characteristics play a key role in determining seismic and seismoelectric signatures of partially saturated media. In this context, most works employ the classic theories of Biot (1956aBiot ( , 1956b and Pride (1994), developed for monosaturated scenarios, using effective fluid properties (e.g., Barrière et al. 2012;Rubino and Holliger 2012;Bordes et al. 2015;Solazzi et al. 2017). As shown by Solazzi et al. (2020), the use of effective fluids in Biot's theory of poroelasticity may lead to errors in the estimation of the energy dissipated due to viscous flow. Our results further stress this point, showing that it is not possible for an effective fluid to account for the complexities associated with viscous coupling effects, which, for example, may induce two inflection points in the dynamic permeability curves. In this context, experimental evidence shows that dissipation due to momentum transfer between immiscible fluid phases at the microscopic scale can have large effects on the macroscopic flow characteristics (Roman et al. 2020). Such effects are further enhanced when relatively small globules of trapped wetting phase fluids interact with a non-wetting phase.
A number of classical papers deal with extensions of Biot's (1956a;b) theory to account for the effects of two immiscible fluid phases (e.g, Santos et al. 1990;Tuncay and Corapcioglu 1997;Lo et al. 2005;Jardani and Revil 2015). Even though some of these works doconsider inertial coupling between saturating fluid phases by means of a coupling parameter (Santos et al. 1990;Lo et al. 2005), these models disregard cross-terms associated with interfacial coupling. As a result, biphasic extensions of Biot's (1956a;b) theory rely on Eq. (1) and, thus, on relative permeabilities which do not account for viscous coupling and, thus, are functions of the saturation only. Further research is needed to combine these formulations with an extended Darcy formulation (Eqs. 4 and 5), where dynamic relative permeability terms accounting for viscous coupling can be employed. In this context, the proposed model may help to develop new and more consistent seismic and seismoelectric theories for partially saturated porous media. This possibility is of particular interest for studying the vadose zone, as unconsolidated soils are characterized by high permeabilities and, thus, present a transition from viscous-to inertia-dominated flow occurring at relatively low frequencies.
respectively. Note that this formulation does not account for the tortousity, which is, however, accounted for in the proposed formulation. The relative permeabilities are therefore given by (Bacri et al. 1990) Figure 7 shows a comparison between Yuster's low-frequency solution (Eqs. 46 to 49) and the proposed frequency-dependent solution (Eqs. 28 to 31) for f = 1 Hz. We consider water as the wetting phase, while the non-wetting phase is chosen such that the viscosity ratio M and the density ratio N are both equal to unity M = N = 1 (Fig. 7a) and M = 1.5 and N = 1 (Fig. 7b). The agreement between the curves proves that the model proposed in this work collapses to that of Yuster's (1951) in the low-frequency viscosity-dominated regime.

Appendix 2: Coefficients of the Proposed Analytical Model
In the following, we give the expressions of the coefficients S , M j , and L j associated with Eqs. (282930) to (31). We do so in terms of the characteristic radius R. Please do, however, note that in our model the characteristic radius is defined as R = √ 8 2 ∕ , and hence, the coefficients can alternatively be expressed in terms of κ and . The coefficients respond to r w n =2S n S w + S n log(S n ) .
(50) M 0 = A + B + C + D + E + F,   Fig. 7 Illustration of the relative permeability terms as functions of the wetting phase saturation S w . We illustrate the predictions of the proposed approach at a frequency of f = 1 Hz (solid lines) and those given by the analytical solution of Yuster (1951) (crosses). We consider a pore radius of 100 μ m, water as the wetting phase, while the non-wetting phase is such that the viscosity ratios M and the density ratios N are a M = N = 1 and b M = 1.5 and N = 1 and where, once again, J m and Y m the mth-order Bessel functions of the first and second kinds, respectively.