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

Santiago G Solazzi, Damien Jougnot, J Germán Rubino, Klaus Holliger

► To cite this version:

Santiago G Solazzi, Damien Jougnot, J Germán Rubino, Klaus Holliger. Dynamic Relative Permeabilities for Partially Saturated Porous Media Accounting for Viscous Coupling Effects: An Analytical Solution. Transport in Porous Media, 2023, 10.1007/s11242-023-01920-w. http://dx.doi.org/10.1007/s11242-023-01920-w

HAL Id: hal-04030527

https://hal.sorbonne-universite.fr/hal-04030527

Submitted on 15 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamic Relative Permeabilities for Partially Saturated Porous Media Accounting for Viscous Coupling Effects: An Analytical Solution

Santiago G. Solazzi<fnref n="1" /> ⒹD成长· Damien Jougnot<fnref n="2" /> ⒹD· J. Germán Rubino<fnref n="3" /> ⒹD· Klaus Holliger<fnref n="1" /> ⒹD

Received: 19 July 2022 / Accepted: 18 February 2023
© The Author(s) 2023

## Abstract

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.

**Keywords** Viscous coupling · Dynamic relative permeability · Frequency-dependent effects · Porous media · Analytical solution · Partial saturation

<figure><img src="image_2.png" /></figure>

<footnote>Santiago G. Solazzi
santiago.solazzi@unil.ch</footnote>

<footnote n="1">¹ Institute of Earth Sciences, University of Lausanne, 1015 Lausanne, Vaud, Switzerland</footnote>

<footnote n="2">² Sorbonne Université, CNRS, EPHE, UMR7619, METIS, 75005 Paris, Île-de-France, France</footnote>

<footnote n="3">³ CONICET, Centro Atómico Bariloche – CNEA, R8402AGP San Carlos de Bariloche, Rio Negro, Argentina</footnote>
# 1 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 $\kappa$ 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 1987, 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 1987, 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 $\omega$. 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* $\kappa(\omega)$, 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 $\kappa(\omega)$ 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 two-phase 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.
## 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.

### 2.1 Classic Two-Phase Flow Equations in Porous Media

The traditional two-phase flow formulation of Darcy's law is given by (Bear 1972)

$$v_i = -\lambda_i \nabla p_i, \quad \text{for } i = [w, n], \qquad (1)$$

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 $\lambda_i$ respond to

$$\lambda_i = \frac{\kappa_{ri} \kappa}{\eta_i}, \quad \text{for } i = [w, n], \qquad (2)$$

where $\eta_i$ is the dynamic viscosity of the fluid, $\kappa$ the intrinsic permeability determined by the structure of the porous medium alone, and $\kappa_{ri}$ 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)

$$\begin{aligned} v_i &= \frac{1}{V_t} \int_{\Omega_p} w_i \, d\Omega_p, \quad \text{for } i = [w, n], \\ &= \phi \langle w_i \rangle, \end{aligned} \qquad (3)$$

where $V_t$ is the total volume of a representative elementary volume (REV) of the porous medium of interest, $\Omega_p$ is the dominium associated with the pore space, $\phi$ is the porosity, and $\langle \cdot \rangle$ denotes the average over the pore space of the REV. Note that $w_i$ is zero for the $\Omega_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.

### 2.2 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 1987, 1990). The generalized two-phase flow equations in porous media respond to (e.g., Dullien and Dong 1996)

$$v_w = -\lambda_{ww} \nabla p_w - \lambda_{wn} \nabla p_n, \qquad (4)$$

$$v_n = -\lambda_{nw} \nabla p_n - \lambda_{nn} \nabla p_n, \qquad (5)$$
<figure><img src="image_3.png" /><figcaption>Fig. 1 Schematic illustration of two-phase flow in a porous medium and blowup of a possible flow patterns. Dark red and blue arrows denote the flow patterns at the macroscopic scale (left panel) and at the pore scale (blowup)</figcaption></figure>

where $\lambda_{ij}$ are the generalized mobilities, which are given by

$$ \lambda_{ij} = \frac{\kappa_{r_j} \kappa}{\eta_j}, \quad \text{for } i, j = [w, n], \qquad (6) $$

with $\kappa_{r_j}$ denoting the generalized relative permeability coefficients. The off-diagonal coefficients $\kappa_{r_j^w}$ and $\kappa_{r_j^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 $\kappa_{r_j}$ 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 ($\nabla p_j$) which is modulated by the corresponding term. For instance, $\kappa_{r_j^n}$ modulates the effect of $\nabla 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 Avraam and Payatakes 1995, 1999; Roman et al. 2016, 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 = \eta_n/\eta_w$ leads to a significant increase in $\kappa_{r_j^n}$ and $\kappa_{r_j^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 $\kappa_{r_j^n}$ and $\kappa_{r_j^n}$ are virtually insensitive to variations in the viscosity ratio $M$. Interestingly, some works indicate that the cross-term mobilities should be identical $\kappa_{r_j^n}/\eta_w \equiv \kappa_{r_j^n}/\eta_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; Avraam and 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 $\kappa(\omega)$.
## 2.3 Dynamic Permeability

In most studies on Darcy flow, the permeability $\kappa$ 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 frequency-dependent permeability (e.g., Johnson et al. 1987; Sheng and Zhou 1988; Zhou and Sheng 1989)

$$v(\omega) = \phi\langle w \rangle = -\frac{\kappa(\omega)}{\eta} \nabla p(\omega), \quad (7)$$

where $\omega = 2\pi f$ is the angular frequency and $f$ denotes the linear frequency.

The behavior of $\kappa(\omega)$ 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 $\omega_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)

$$\omega_c \approx \frac{2v}{R_p^2}, \quad (8)$$

with $v = \eta/\rho$ denoting the kinematic viscosity and $\rho$ the fluid density (Johnson et al. 1987). Please note that the dynamic permeability $\kappa(\omega)$ does not only depend on the properties of the porous medium alone, but also depend on the physical properties of the saturating fluid ($\eta$ and $\rho$), and on the angular frequency $\omega$.

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 frequency-dependent relative permeability functions $\kappa_{ri}(\omega, S_w)$. For the latter to be valid, the derived $\kappa_{ri}(\omega, S_w)$ functions need to converge to their real-valued counterparts for sufficiently low frequencies.

## 3 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.
## 3.1 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 \le r \le R_1$), while the wetting phase saturates the annular space between the non-wetting phase and the pore wall ($R_1 < r \le 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 $\phi$, and the permeability $\kappa$ respond to

$$S_n = \frac{R_1^2}{R^2}, \qquad (9)$$

$$S_w = \frac{R^2 - R_1^2}{R^2}, \qquad (10)$$

$$\phi = \frac{N_p \tau \pi R^2}{L^2}, \qquad (11)$$

<figure><img src="image_4.png" /><figcaption>Fig. 2 Illustration of the considered analog pore-scale model. a Representative elementary volume (REV) of the porous medium of interest and b a blowup of the annular saturation pattern. Within the capillaries, light blue denotes the wetting and dark red the non-wetting fluid phase</figcaption></figure>
$$ \kappa = \frac{\phi R^2}{\tau^2 8}, \qquad (12) $$

respectively. In Eqs. (11) and (12), L is the sidelength of the cubic REV, $\tau = 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.

## 3.2 Mathematical Formulation of the Problem

### 3.2.1 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)

$$ \rho \left[ \frac{\partial \mathbf{w}}{\partial t} + (\mathbf{w} \cdot \nabla) \mathbf{w} \right] = -\nabla p + \eta \Delta \mathbf{w} + \left( \zeta + \frac{1}{3} \eta \right) \nabla (\nabla \cdot \mathbf{w}) + \mathbf{f}_{\sigma}, \quad (13) $$

where $\mathbf{w}$ is the fluid velocity at the microscopic scale, $p$ the fluid pressure, $t$ the time, and $\zeta$ the secondary viscosity coefficient. Here, $\mathbf{f}_{\sigma}$ denotes the surface tension force, which, assuming constant interfacial tension $\sigma$, responds to (e.g., Tryggvason et al. 2011)

$$ \mathbf{f}_{\sigma} = \sigma \gamma \mathbf{n} \delta_{S}, \qquad (14) $$

where $\mathbf{n}$ is the unit normal vector with respect to the interface, $\gamma$ is the mean curvature of the interface, and $\delta_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 ($\nabla \cdot \mathbf{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

$$ \rho \frac{\partial \mathbf{w}}{\partial t} + \rho(\mathbf{w} \cdot \nabla)\mathbf{w} = -\nabla p + \eta \Delta \mathbf{w} + \mathbf{f}_{\sigma}. \qquad (15) $$

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, $\mathbf{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)

$$ \frac{\partial p}{\partial r} = -\sigma \gamma \delta_S, \qquad (16) $$
$$ \frac{\rho}{\eta} \frac{\partial w(r)}{\partial t} = -\frac{1}{\eta} \frac{\partial p}{\partial z} + \frac{\partial^2 w(r)}{\partial r^2} + \frac{1}{r} \frac{\partial w(r)}{\partial r}. \quad (17) $$

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, $p_n = p_w + \frac{\sigma}{R_1}$ (e.g., Coward et al. 1995). Note that this condition implies the equality $\nabla p_w = \nabla 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

$$ \mathcal{X} = \frac{\partial p}{\partial z}, \quad (18) $$

and consider that both the relative fluid velocity and the forcing term present an oscillatory form given by $w = \hat{w} e^{-i\omega t}$ and $\mathcal{X} = \hat{\mathcal{X}} e^{-i\omega t}$, respectively, with i denoting the imaginary unit. Then, Eq. (17) can be expressed in the space frequency domain as

$$ \left[ \frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial}{\partial r} + k^2 \right] \hat{w}(r) = \frac{\hat{\mathcal{X}}}{\eta}, \quad (19) $$

where $k = \sqrt{i\omega\rho/\eta}$ can be conceptualized as the wave number of the oscillatory process (Johnson et al. 1987). By taking $\bar{w}(r) = \hat{w} - \frac{\hat{\mathcal{X}}}{k^2\eta}$, the corresponding equation reduces to a Bessel-type differential equation

$$ \left[ \frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial}{\partial r} + k^2 \right] \bar{w}(r) = 0. \quad (20) $$

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

$$ \hat{w}(r, \omega) = \begin{cases} \hat{w}_n(r), & \text{for } r = [0, R_1], \\ \hat{w}_w(r), & \text{for } r = [R_1, R]. \end{cases} \quad (21) $$

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)

$$ \hat{w}_n(r, \omega) = \frac{\hat{\mathcal{X}}_n}{\eta_n k_n^2} + \mathcal{C}_1^n J_0(k_n r) + \mathcal{C}_2^n Y_0(k_n r), \quad \text{for } r = [0, R_1], \quad (22) $$

$$ \hat{w}_w(r, \omega) = \frac{\hat{\mathcal{X}}_w}{\eta_w k_w^2} + \mathcal{C}_1^w J_0(k_w r) + \mathcal{C}_2^w Y_0(k_w r), \quad \text{for } r = [R_1, R], \quad (23) $$

where $\mathcal{C}_1^j$ and $\mathcal{C}_2^j$, 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.
### 3.2.2 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 $\kappa_{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 $\hat{w}_n$ value at $r=0$ must be finite. Given that second-kind Bessel functions present a singularity at the origin, that is, $\lim_{r\to 0} Y_0(k_n r) = \infty$, the $\mathcal{C}_2^n$ term has to be null for $\hat{w}_n$ to be non-diverging in Eq. (22). Conversely, the no-slip assumption at the pore walls implies that $\hat{w}_n|_{(r=R)} = 0$, which, considering Eq. (22), reduces to

$$ \mathcal{C}_1^w J_0(k_w R) + \mathcal{C}_2^w Y_0(k_w R) = -\frac{\hat{\chi}_w}{\eta_w k_w^2}. \quad (24) $$

On the other hand, we have the continuity of the fluid shear stress $\tau_j = \eta_j d\hat{w}_j/dr$ at the interface $r = R_1$, which results in

$$ \eta_w k_w \mathcal{C}_1^w J_1(k_w R_1) + \eta_w k_w \mathcal{C}_2^w Y_1(k_w R_1) - \eta_n k_n \mathcal{C}_1^n J_1(k_n R_1) = 0, \quad (25) $$

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 $(\hat{w}_w - \hat{w}_n)|_{(r=R_1)} = 0$, which, using Eqs. (22) and (23), yields

$$ \mathcal{C}_1^w J_0(k_w R_1) + \mathcal{C}_2^w Y_0(k_w R_1) - \mathcal{C}_1^n J_0(k_n R_1) = \left( \frac{\hat{\chi}_n}{\eta_n k_n^2} - \frac{\hat{\chi}_w}{\eta_w k_w^2} \right). \quad (26) $$

Equations (24), (25), and (26) constitute a linear system of equations that can be solved analytically to obtain $\mathcal{C}_1^w$, $\mathcal{C}_2^w$, and $\mathcal{C}_1^n$. 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).

## 3.3 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 $\langle \hat{w}_j \rangle$ over the cross section of the capillaries and multiplying the latter with the porosity $\phi$

$$ \hat{v}_w(\omega) = \frac{\phi}{\pi R^2} 2\pi \int_{R_1}^{R} r \hat{w}_w(r, \omega) dr, \quad \hat{v}_n(\omega) = \frac{\phi}{\pi R^2} 2\pi \int_{0}^{R_1} r \hat{w}_n(r, \omega) dr. \quad (27) $$

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
$$ \lambda_{ww}(\omega) = -\frac{\phi M_0}{\pi R^2 M_1}, \qquad (28) $$

$$ \lambda_{wn}(\omega) = -\frac{\phi S}{\pi R^2 M_1}, \qquad (29) $$

$$ \lambda_{nw}(\omega) = -\frac{\phi (\mathcal{L}_1 - \mathcal{L}_2)}{\pi R^2 M_2}, \qquad (30) $$

$$ \lambda_{nn}(\omega) = -\frac{\phi}{\pi R^2} \left( \frac{\mathcal{L}_3}{M_3} + 1 \right) \mathcal{L}_4. \qquad (31) $$

The expressions for the coefficients $S$, $M_j$, and $\mathcal{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

$$ \kappa_{jw}(\omega) = \frac{\eta_w}{\kappa} \lambda_{ww}(\omega), \qquad (32) $$

$$ \kappa_{rw}(\omega) = \frac{\eta_n}{\kappa} \lambda_{wn}(\omega), \qquad (33) $$

$$ \kappa_{r_n}(\omega) = \frac{\eta_w}{\kappa} \lambda_{nw}(\omega), \qquad (34) $$

$$ \kappa_{r_n}(\omega) = \frac{\eta_n}{\kappa} \lambda_{nn}(\omega). \qquad (35) $$

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 = \sqrt{8\tau^2\kappa/\phi}$, and thus, the corresponding equations can alternatively be expressed in terms of $\kappa$ and $\phi$. It is important to remark here that, as evidenced in Eqs. (50) to (70), the kinematic viscosities of the fluids $v_n = \eta_n/\rho_n$ and $v_w = \eta_w/\rho_w$ play a key role in determining the dynamic relative permeabilities.

# 4 Results

## 4.1 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 $\kappa_j(\omega)$ (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 $\tau = 1$ and $R = 100 \ \mu\text{m}$. We take water and oil as the saturating fluids, the former having a density $\rho_w = 1000 \ \text{kg m}^{-3}$, and the latter one of $\rho_n = 700 \ \text{kg m}^{-3}$, thus resulting in a density ratio of
$N = 0.7$. In the following, the water viscosity is taken as $\eta_w = 0.001$ Pa s. For oil, we consider two viscosities: (i) $\eta_n = 0.01$ Pa s and (ii) $\eta_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 $\kappa_{r_i}(\omega)$ 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 $\kappa_{r_i}^w$ and $\kappa_{r_n}^w$, which modulate the Darcy velocity of water $v_w$. The bottom row (Fig. 3c, d) corresponds to the relative permeabilities $\kappa_{r_i}^w$ and $\kappa_{r_n}^w$, which modulate the Darcy velocity of oil $v_\eta$. 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 $\nabla p_w$ and $\nabla 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 $\kappa_{r_i}^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, $\kappa_{r_i}^w$ is not affected by viscous coupling

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Real part of the dynamic relative permeabilities as functions of saturation $S_w$ for frequencies of $f = 1$ Hz (black lines), $f = 100$ Hz (blue lines), and $f = 1000$ Hz (red lines). 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 squared markers) and $M = 10$ (solid lines). The medium is characterized by a tortuosity and a pore radius of $\tau = 1$ and $R = 100$ µm, respectively</figcaption></figure>
effects. Conversely, $\kappa_{r^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 $\kappa_{r^n}$. Also for $f = 1$ Hz (black lines), we note a large variation in the crossed term $\kappa_{r^0}$ (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 $\kappa_{r^0}$ 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 \neq 1$, the relation $\kappa_{r^0}/\eta_w \equiv \kappa_{r^0}/\eta_n$ holds in all cases.

Let us now analyze the response of $\Re\{\kappa_{r^j}(S_w, \omega)\}$ 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 $\Re\{\kappa_{r^j}(S_w, \omega)\}$ 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 $\omega_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 $\Re\{\kappa_{r^j}\}$ 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 $\Re\{\kappa_{r^0}\}$ (Fig. 3b) and $\Re\{\kappa_{r^0}\}$ (Fig. 3c) show negative values for $f = 1000$Hz and $M = 10$ (red solid lines) for $0.6^0 < 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 $\kappa_{r^j}(S_w, \omega)$ (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, $\mathfrak{T}\{\kappa_{r^j}(S_w, \omega)\}$ 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), $\mathfrak{T}\{\kappa_{r^j}\}$ 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 $\mathfrak{T}\{\kappa_{r^0}\}$ 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 $\mathfrak{T}\{\kappa_{r^0}\}$ 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.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Imaginary part of the dynamic relative permeabilities as functions of saturation for frequencies of $f = 1$ Hz (black lines), $f = 100$ Hz (blue lines), and $f = 1000$ Hz (red lines). 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 $\tau = 1$ and $R = 100$ µm, respectively</figcaption></figure>

## 4.2 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 $\mathfrak{R}\{\kappa_{rw}\}$ (Fig. 5a) increases with $S_w$ while $\mathfrak{R}\{\kappa_{rn}\}$ (Fig. 5d) decreases with $S_w$, for both $M = 1$ (dashed lines) and $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 $\mathfrak{R}\{\kappa_{rw}\}$ and $\mathfrak{R}\{\kappa_{rn}\}$, which are the terms that modulate flow induced by the water pressure gradient $\nabla p_w$. Conversely, we evidence a large increase in the values of $\mathfrak{R}\{\kappa_{rn}\}$ (Fig. 5b) and $\mathfrak{R}\{\kappa_{rw}\}$ (Fig. 5d) with increasing $M$. These terms modulate flow by $\nabla 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 $\omega_c$. For increasing values of $\eta_n$ and, thus, of $M$, $\omega_c$ moves toward to lower frequencies. Figure 5a shows that,
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 5** 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 $\tau = 1$ and $R = 100$ µm, respectively</figcaption></figure>

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 non-diagonal 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 relative 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 $\mathfrak{T}\{\kappa_{r^n}\}$ (Fig. 6b) and $\mathfrak{T}\{\kappa_{r^n}\}$ (Fig. 6d) for $M = 10$ (solid
<figure><img src="image_11.png" /><figcaption>**Fig. 6** Imaginary 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 $\tau = 1$ and $R = 100$ µm, respectively</figcaption></figure>

lines) compared to $M = 1$ (dashed lines), which is expected as the viscosity of the non-wetting phase is increased tenfold to increase $M$.

# 5 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 $\kappa_{r_n}$ and $\kappa_{r_n^2}$, which modulate the effects of the non-wetting fluid pressure gradient $\nabla p_n$ on the flow, are very sensitive to the viscosity ratio, exhibiting variations with regard to the amplitudes and the critical
frequencies. On the other hand, the dynamic relative permeability terms $\kappa_{r,w}^n$ and $\kappa_{r,g}^n$ that modulate the effects of the wetting phase pressure gradient $\nabla 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 $\kappa_{r,w}^n$ and $\kappa_{r,g}^n$ 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 (1956a, 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 do consider 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.

## Appendix 1: Steady-State Flow: Yuster's Problem

Yuster (1951) and Bacri et al. (1990) proposed an analogous model consisting of a porous structure composed of capillary tubes of a single characteristic radius saturated concentrically with two different fluid phases as the one described in Fig. 2. Recall that the non-wetting fluid saturates the center of the capillary ($0 \le r \le R_1$) while the wetting phase saturates the annular space between the non-wetting phase and the pore wall ($R_1 < r \le R$). Under the assumption that the relative fluid velocities are dependent upon the radial coordinate only, the steady-state flow responds to (Bacri et al. 1990)

$$w_i = \frac{1}{4\eta_i} \chi_i r^2 + A_i \log(r) + B_i, \quad \text{with } i = w, n, \qquad (36)$$

where $A_i$ and $B_i$ are integration constants and are $\chi_i = \frac{\partial p_i}{\partial z}$. The boundary conditions are the same as the ones imposed in our approach, that is,

$$w_n|_{r=R_1} = w_n|_{r=R_1}, \qquad (37)$$

$$w_w|_{r=R} = 0, \qquad (38)$$

$$\eta_n \frac{\partial w_n}{\partial r} |_{r=R_1} = \eta_w \frac{\partial w_w}{\partial r} |_{r=R_1}, \qquad (39)$$

$$\frac{\partial w_n}{\partial r} |_{r=0} = 0. \qquad (40)$$

These equations, in combination with Eqs. (4), (5), and (27), allow us to determine the steady-state Darcy velocities

$$\begin{aligned} v_w &= -\frac{K}{\eta_w} [S_w - 3S_w S_n - 2S_n^2 \log(S_n)] \chi_w - \frac{K}{\eta_n} 2MS_n [S_w + S_n \log(S_n)] \chi_n, \\ v_n &= -\frac{K}{\eta_w} 2S_n [S_w + S_n \log(S_n)] \chi_w - \frac{K}{\eta_n} S_n^2 [1 - 2M \log(S_n)] \chi_n, \end{aligned} \qquad (41)$$

where the non-wetting phase saturation, the wetting phase saturation, the fluid viscosity ratio, and the permeability are given by
$$S_n = \frac{R_1^2}{R^2}, \qquad (42)$$

$$S_w = \frac{R^2 - R_1^2}{R^2}, \qquad (43)$$

$$M = \frac{\eta_n}{\eta_w}, \qquad (44)$$

$$\kappa = \frac{\phi R^2}{8}, \qquad (45)$$

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)

$$\kappa_{rw} = [S_w - 3S_w S_n - 2S_n^2 \log(S_n)], \qquad (46)$$

$$\kappa_{rw^n} = 2MS_n[S_w + S_n \log(S_n)], \qquad (47)$$

$$\kappa_{rw^n} = S_n^2(1 - 2M \log(S_n)), \qquad (48)$$

$$\kappa_{rw} = 2S_n(S_w + S_n \log(S_n)). \qquad (49)$$

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 = \sqrt{8\tau^2\kappa/\phi}$, and hence, the coefficients can alternatively be expressed in terms of $\kappa$ and $\phi$. The coefficients respond to

$$M_0 = A + B + C + D + E + F, \qquad (50)$$

$$A = \eta_w k_n k_w^2 \pi J_0[k_n S_n^{1/2} R](A_1 + A_2), \qquad (51)$$
<figure><img src="image_12.png" /><figcaption>**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$</figcaption></figure>

$$A_1 = J_1[k_w S_n^{1/2} R](k_w R^2 S_w Y_0[k_w R] - 2RY_1[k_w R]), \quad (52)$$

$$A_2 = Y_1[k_w S_n^{1/2} R](-k_w R^2 S_w Y_0[k_w R] + 2RJ_1[k_w R]), \quad (53)$$

$$B = \eta_n k_n^2 k_w \pi (k_w R^2 S_w Y_0[k_w R] - 2RJ_1[k_w R]) J_1[k_n S_n^{1/2} R] Y_0[k_n S_n^{1/2} R], \quad (54)$$

$$C = \eta_n k_n^2 k_w \pi J_0[k_w S_n^{1/2} R] J_1[k_n S_n^{1/2} R] (-k_w R^2 S_w Y_0[k_w R] + 2RY_1[k_w R]), \quad (55)$$

$$D = -2\eta_n k_n^2 k_w \pi S_n^{1/2} R J_1[k_n S_n^{1/2} R] J_1[k_w S_n^{1/2} R] Y_0[k_w R], \quad (56)$$

$$\mathcal{E} = 2\eta_n k_n^2 k_w \pi S_n^{1/2} R J_1[k_n S_n^{1/2} R] J_0[k_w R] Y_1[k_w S_n^{1/2} R], \quad (57)$$

$$\mathcal{F} = 8\eta_n k_n^2 J_1[k_n S_n^{1/2} R], \quad (58)$$

$$S = S_1 + S_2, \quad (59)$$

$$S_1 = -4\eta_w k_w^2 J_1[k_n S_n^{1/2} R], \quad (60)$$

$$S_2 = 2\eta_w k_w^3 \pi S_n^{1/2} R J_1[k_n S_n^{1/2} R] (J_1[k_n S_n^{1/2} R] Y_0[k_n R] - J_0[k_n R] Y_1[k_n S_n^{1/2} R]), \quad (61)$$

$$\mathcal{L}_1 = \frac{4S_n^{1/2} R}{(k_w k_n)} J_1[k_n S_n^{1/2} R], \quad (62)$$
$$ \mathcal{L}_2 = \frac{2\pi S_n R^2}{k_n} J_1[k_n S_n^{1/2} R] (J_1[k_w S_n^{1/2} R] Y_0[k_w R] - J_0[k_w R] Y_1[k_w S_n^{1/2} R]), \quad (63) $$

$$ \mathcal{L}_3 = 2\eta_w k_w J_1[k_n S_n^{1/2} R] (J_1[k_w S_n^{1/2} R] Y_0[k_w R] - J_0[k_w R] Y_1[k_w S_n^{1/2} R]), \quad (64) $$

$$ \mathcal{L}_4 = \frac{\pi S_n R^2}{(\eta_n k_n^2)}, \quad (65) $$

and

$$ \mathcal{M}_1 = - \eta_w k_n k_w^4 (Q + \mathcal{R}), \quad (66) $$

$$ \mathcal{M}_2 = k_w S_n^{1/2} R (Q + \mathcal{R}), \quad (67) $$

$$ \mathcal{M}_3 = k_n S_n^{1/2} R (Q + \mathcal{R}), \quad (68) $$

$$ Q = \eta_n k_n J_1[k_n S_n^{1/2} R] (J_0[k_w S_n^{1/2} R] Y_0[k_w R] - J_0[k_w R] Y_0[k_w S_n^{1/2} R]), \quad (69) $$

$$ \mathcal{R} = \eta_w k_w J_0[k_n S_n^{1/2} R] (J_0[k_w R] Y_1[k_w S_n^{1/2} R] - J_1[k_w S_n^{1/2} R] Y_0[k_w R]), \quad (70) $$

where, once again, $J_m$ and $Y_m$ the mth-order Bessel functions of the first and second kinds, respectively.

## Appendix 3: Notation

<table><tr><td>f</td><td>Linear frequency</td></tr><tr><td>f<sub>σ</sub></td><td>Surface tension force</td></tr><tr><td>J<sub>m</sub></td><td>Bessel function of the first kind of order m</td></tr><tr><td>k<sub>i</sub></td><td>Wave number of the oscillatory process</td></tr><tr><td>L</td><td>Sidelength of the REV</td></tr><tr><td>L<sub>p</sub></td><td>Pore length</td></tr><tr><td>M</td><td>Viscosity ratio between the pore fluids</td></tr><tr><td>N</td><td>Density ratio between the pore fluids</td></tr><tr><td>N<sub>p</sub></td><td>Number of pores in the REV</td></tr><tr><td>p<sub>i</sub></td><td>Pressure of the fluid phase i</td></tr><tr><td>S<sub>i</sub></td><td>Saturation of the fluid phase i</td></tr><tr><td>r</td><td>Radial coordinate</td></tr><tr><td>R</td><td>Pore radius for the annular flow model</td></tr><tr><td>R<sub>p</sub></td><td>Characteristic pore radius</td></tr><tr><td>R<sub>l</sub></td><td>Inner radius for the annular flow model</td></tr><tr><td>v<sub>i</sub></td><td>Darcy velocity of the fluid phase i</td></tr><tr><td>θ<sub>i</sub></td><td>Fourier transform of the Darcy velocity of phase i</td></tr><tr><td>V<sub>i</sub></td><td>Volume of the REV</td></tr></table>
<table><tr><td>w<sub>i</sub></td><td>Pore-scale velocity of the fluid phase i</td></tr><tr><td>ŵ<sub>i</sub></td><td>Fourier transform of the pore-scale velocity of phase i</td></tr><tr><td>Y<sub>m</sub></td><td>Bessel function of the second kind of order m</td></tr><tr><td>z</td><td>Axial coordinate</td></tr><tr><td>γ</td><td>Mean curvature of the fluid-fluid interface</td></tr><tr><td>δ<sub>S</sub></td><td>Delta function defined at the fluid-fluid interface</td></tr><tr><td>ζ</td><td>Secondary viscosity coefficient</td></tr><tr><td>η<sub>i</sub></td><td>Dynamic viscosity of the fluid phase i</td></tr><tr><td>κ</td><td>Absolute permeability</td></tr><tr><td>κ<sub>ri</sub></td><td>Relative permeability of fluid phase i</td></tr><tr><td>κ<sub>j</sub></td><td>Relative permeability of fluid phase i with respect to phase j</td></tr><tr><td>λ<sub>i</sub></td><td>Mobility of the fluid phase i</td></tr><tr><td>λ<sub>ij</sub></td><td>Generalized mobility of fluid phase i with respect to phase j</td></tr><tr><td>ν<sub>i</sub></td><td>Kinematic viscosity of the fluid phase i</td></tr><tr><td>ρ<sub>i</sub></td><td>Density of the fluid phase i</td></tr><tr><td>σ</td><td>Capillary tension</td></tr><tr><td>τ</td><td>Tortuosity</td></tr><tr><td>φ</td><td>Porosity</td></tr><tr><td>ω</td><td>Angular frequency</td></tr><tr><td>ω<sub>c</sub></td><td>Biot's characteristic frequency</td></tr><tr><td>Ω<sub>p</sub></td><td>Dominium associated with the pore space</td></tr><tr><td>Γ<sup>w</sup><sub>1</sub>, Γ<sup>w</sup><sub>2</sub>, Γ<sup>n</sup><sub>1</sub>, Γ<sup>n</sup><sub>2</sub></td><td>Integration constants</td></tr><tr><td>L<sub>1</sub>, L<sub>2</sub>, L<sub>3</sub>, L<sub>4</sub></td><td>Coefficients of the analytical model</td></tr><tr><td>M<sub>0</sub>, M<sub>1</sub>, M<sub>2</sub>, M<sub>3</sub></td><td>Coefficients of the analytical model</td></tr><tr><td>S</td><td>Coefficient of the analytical model</td></tr><tr><td>χ<sub>i</sub></td><td>Pressure gradient of phase i with respect to the axial direction</td></tr><tr><td>χ̂<sub>i</sub></td><td>Fourier transform of pressure gradient of phase i</td></tr><tr><td>w</td><td>Subindex indicating wetting fluid association</td></tr><tr><td>n</td><td>Subindex indicating non-wetting fluid association</td></tr></table>

**Acknowledgements** The authors thank Steffen Berg and three anonymous reviewers for lucid comments and helpful suggestions that greatly improved this manuscript.

**Funding** Open access funding provided by University of Lausanne. This project has received funding from the European Union's HORIZON 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 101007851. Experimental and observational data were neither used nor created for this research. Mathematical expressions to replicate the results of this paper are available in the paper.

**Declarations**

**Competing interests** The authors declare no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
References

Abramowitz, M., Stegun, I.A.: Handbook of Mathematical Functions: With Formulas, Graphs and Math-ematical Tables. Dover Publications, Mineola (1965)

Auriault, J.L., Borne, L., Chambon, R.: Dynamics of porous saturated media, checking of the generalized law of darcy. J. Acoust. Soc. Am. **77**(5), 1641–1650 (1985)

Auriault, J.L., Lebaigue, O., Bonnet, G.: Dynamics of two immiscible fluids flowing through deformable porous media. Transp. Porous Media **4**(2), 105–128 (1989)

Avraam, D., Payatakes, A.: Generalized relative permeability coefficients during steady-state two-phase flow in porous media, and correlation with the flow mechanisms. Transp. Porous Media **20**(1), 135–168 (1995)

Avraam, D., Payatakes, A.: Flow mechanisms, relative permeabilities, and coupling effects in steady-state two-phase flow through porous media. The case of strong wettability. Ind Eng Chem Res **38**(3), 778–786 (1999)

Ayub, M., Bentsen, R.G.: Interfacial viscous coupling: a myth or reality? J. Pet. Sci. Eng. **23**(1), 13–26 (1999)

Bachmat, Y., Bear, J.: Macroscopic modelling of transport phenomena in porous media. 1: The continuum approach. Transp. Porous Media **1**(3), 213–240 (1986)

Bacri, J.C., Chaouche, M., Salin, D.: Modèle simple de perméabilités relatives croisées. Comptes rendus de l'Académie des sciences Série 2, Mécanique, Physique, Chimie, Sciences de l'univers, Sciences de la Terre **311**(6):591–597 (1990)

Badchin, A.J., Yuan, J.Y.: On the capillary coupling between two phases in a droplet train model. Transp. Porous Media **26**(2), 225–228 (1997)

Barrière, J., Bordes, C., Brito, D., et al.: Laboratory monitoring of P waves in partially saturated sand. Geo-phys. J. Int. **191**(3), 1152–1170 (2012)

Bear, J.: Dynamics of Fluids in Porous Media. American Elsevier, New York (1972)

Bentsen, R.G.: Effect of momentum transfer between fluid phases on effective mobility. J. Pet. Sci. Eng. **21**(1–2), 27–42 (1998)

Bentsen, R.G.: The physical origin of interfacial coupling in two-phase flow through porous media. Transp. Porous Media **44**(1), 109–122 (2001)

Biot, M.A.: Theory of propagation of elastic waves in a fluid saturated porous solid. I. Low frequency range. J. Acoust. Soc. Am. **28**, 168–178 (1956a). https://doi.org/10.1063/1.1721956

Biot, M.A.: Theory of propagation of elastic waves in a fluid saturated porous solid. II. Higher frequency range. J. Acoust. Soc. Am. **28**, 179–191 (1956b)

Biot, M.A.: Mechanics of deformation and acoustic propagation in porous media. J. Appl. Phys. **33**(4), 1482–1498 (1962). https://doi.org/10.1063/1.1728759

Blunt, M.J.: Multiphase Flow in Permeable Media: A Pore-Scale Perspective. Cambridge University Press, Cambridge (2017)

Bolt, G., Groenevelt, P.: Coupling phenomena as a possible cause of “non-darcian” behaviour of water in soil. Hydrol Sci J **14**(2), 17–28 (1969)

Bordes, C., Sénéchal, P., Barrière, J., et al.: Impact of water saturation on seismoelectric transfer functions: a laboratory study of coseismic phenomenon. Geophys. J. Int. **200**(3), 1317–1335 (2015)

Burdine, N.: Relative permeability calculations from pore size distribution data. J. Petrol. Tech. **5**(03), 71–78 (1953)

Chapman, A., Higdon, J.: Oscillatory stokes flow in periodic porous media. Phys. Fluids A Fluid Dyn. **4**(10), 2099–2116 (1992)

Charlaix, E., Kushnick, A., Stokes, J.: Experimental study of dynamic permeability in porous media. Phys. Rev. Lett. **61**(14), 1595 (1988)

Coward, A.V., Papageorgiou, D.T., Smyrlis, Y.S.: Nonlinear stability of oscillatory core-annular flow: a generalized Kuramoto-Sivashinsky equation with time periodic coefficients. Z. Angew. Math. Phys. **46**(1), 1–39 (1995)

de Gennes, P.G.: Theory of slow biphasic flows in porous-media. Phys. Chem. Hydrol. **4**(2), 175–185 (1983)

de la Cruz, V., Spanos, T.: Mobilization of oil ganglia. AIChE J. **29**(5), 854–858 (1983)

Dullien, F., Dong, M.: Experimental determination of the flow transport coefficients in the coupled equa-tions of two-phase flow in porous media. Transp. Porous Med. **25**(1), 97–120 (1996)

Ehrlich, R.: Viscous coupling in two-phase flow in porous media and its effect on relative permeabilities. Transp. Porous Media **11**(3), 201–218 (1993)

Guarracino L (2007) Estimation of saturated hydraulic conductivity $k_s$ from the van genuchten shape param-eter $\alpha$. Water Resour Res 43(11)
Gunstensen, A.K., Rothman, D.H.: Lattice-Boltzmann studies of immiscible two-phase flow through porous media. J. Geophys. Res. Solid Earth **98**(B4), 6431–6441 (1993)

Huang, H., Lu, X.y.: Relative permeabilities and coupling effects in steady-state gas-liquid flow in porous media: a lattice boltzmann study. Phys. Fluids 21(9):092,104 (2009)

Huber, C., Su, Y.: A pore-scale investigation of the dynamic response of saturated porous media to transient stresses. Geofluids **15**(1–2), 11–23 (2015)

Jardani, A., Revil, A.: Seismoelectric couplings in a poroelastic material containing two immiscible fluid phases. Geophys. J. Int. **202**(2), 850–870 (2015)

Johnson, D.L., Koplik, J., Dashen, R.: Theory of dynamic permeability and tortuosity in fluid-saturated porous media. J. Fluid Mech. **176**, 379–402 (1987)

Kalaydjian, F.: A macroscopic description of multiphase flow in porous media involving spacetime evolution of fluid/fluid interface. Transp. Porous Media **2**(6), 537–552 (1987)

Kalaydjian, F.: Origin and quantification of coupling between relative permeabilities for two-phase flows in porous media. Transp. Porous Media **5**(3), 215–229 (1990)

Kozeny, J.: Über Kapillare Leitung des Wassers im Boden. Wien Akad Wiss **136**(2a), 271–306 (1927)

Landau, L.D., Lifshitz, E.M.: Course of Theoretical Physics, vol. 6 Fluid Mechanics, vol. 6. Pergamon Press, Exeter (1959)

Li, H., Pan, C., Miller, C.T.: Pore-scale investigation of viscous coupling effects for two-phase flow in porous media. Phys. Rev. E **72**(2):026,705 (2005)

Li, J.X., Rezaee, R., Müller, T.M. et al: Pore size distribution controls dynamic permeability. Geophys. Res. Lett. **48**(5):e2020GL090,558 (2021)

Lo, W.C., Sposito, G., Majer, E.: Wave propagation through elastic porous media containing two immiscible fluids. Water Resour. Res. (2005). https://doi.org/10.1029/2004WR003162

Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. **12**(3), 513–522 (1976)

Müller, T.M., Sahay, P.N.: Stochastic theory of dynamic permeability in poroelastic media. Phys Rev E **84**(2):026,329 (2011)

Nguyen, N.A., Jougnot, D., Thanh, L.D., et al.: Predicting water flow in fully and partially saturated porous media: a new fractal-based permeability model. Hydrogeol. J. **29**(6), 2017–2031 (2021)

Paceman, D.W.: Fundamentals of Numerical Reservoir Simulation. Elsevier, New York (1977)

Picchi, D., Battiato, I.: The impact of pore-scale flow regimes on upscaling of immiscible two-phase flow in porous media. Water Resour. Res. **54**(9), 6683–6707 (2018)

Pride, S.R.: Governing equations for the coupled electromagnetics and acoustics of porous media. Phys. Rev. B **50**(21):15,678 (1994)

Pride, S.R.: Relationships between seismic and hydrological properties. In: Rubin Y, Hubbard S (eds) Hydrogeophysics, chap 9, pp. 253–290. Springer, Berlin (2005)

Pride, S.R., Morgan, F.D., Gangi, A.F.: Drag forces of porous-medium acoustics. Phys. Rev. B **47**(9), 4964 (1993)

Rakotomalala, N., Salin, D., Yortsos, Y.: Viscous coupling in a model porous medium geometry: effect of fluid contact area. Appl. Sci. Res. **55**(2), 155–169 (1995)

Revil, A., Mahardika, H.: Coupled hydromechanical and electromagnetic disturbances in unsaturated porous materials. Water Resour. Res. **49**(2), 744–766 (2013)

Roman, S., Soulaine, C., AlSaud, M.A., et al.: Particle velocimetry analysis of immiscible two-phase flow in micromodels. Adv. Water Resour. **95**, 199–211 (2016)

Roman, S., Soulaine, C., Kovscek, A.R.: Pore-scale visualization and characterization of viscous dissipation in porous media. J. Colloid Interface Sci. **558**, 269–279 (2020)

Rose, W.: Coupling coefficients for two-phase flow in pore spaces of simple geometry. Transp. Porous Media **5**(1), 97–102 (1990)

Rothman, D.H.: Macroscopic laws for immiscible two-phase flow in porous media: results from numerical experiments. J. Geophys. Res. Solid Earth **95**(B6), 8663–8674 (1990)

Rubino, J.G., Holliger, K.: Seismic attenuation and velocity dispersion in heterogeneous partially saturated porous rocks. Geophys. J. Int. **188**, 1088–1102 (2012). https://doi.org/10.1111/j.1365-246X.2011.05291.x

Santos, J.E., Douglas, J., Jr., Corberó, J., et al.: A model for wave propagation in a porous medium saturated by a two-phase fluid. J. Acoust. Soc. Am. **87**(4), 1439–1448 (1990)

Scott, P., Rose, W.: An explanation of the Yuster effect. J. Pet. Technol. **5**(11), 19–20 (1953)

Sheng, P., Zhou, M.Y.: Dynamic permeability in porous media. Phys. Rev. Lett. **61**(14), 1591 (1988). https://doi.org/10.1103/PhysRevLett.61.1591

Smeulders, D., Eggels, R., Van Dongen, M.: Dynamic permeability: reformulation of theory and new experimental and numerical data. J. Fluid Mech. **245**, 211–227 (1992)
Solazzi, S.G., Guarracino, L., Rubino, J.G., et al.: Modeling forced imbibition processes and the associated seismic attenuation in heterogeneous porous rocks. J. Geophys. Res. Solid Earth **122**(11), 9031–9049 (2017)

Solazzi, S.G., Rubino, J.G., Jougnot, D., et al.: Dynamic permeability functions for partially saturated porous media. Geophys. J. Int. **221**(2), 1182–1189 (2020)

Tryggvason, G., Scardovelli, R., Zaleski, S.: Direct Numerical Simulations of Gas-Liquid Multiphase Flows. Cambridge University Press, Cambridge (2011)

Tuncay, K., Corapcioglu, M.Y.: Wave propagation in poroelastic media saturated by two fluids. J. Appl. Mech. **64**(2), 313–320 (1997), https://doi.org/10.1115/1.2787309

Whitaker, S.: Flow in porous media II: the governing equations for immiscible, two-phase flow. Transp. Porous Media **1**(2), 105–125 (1986)

Wolfram Research, Inc. (2021) Mathematica, Version 12.3.1, Champaign, Illinois. https://www.wolfram.com/mathematica

Yuster, S.: Theoretical considerations of multiphase flow in idealized capillary systems. In: 3rd World Petroleum Congress, OnePetro (1951)

Zhou, M.Y., Sheng, P.: First-principles calculations of dynamic permeability in porous media. Phys. Rev. B **39**(16):12,027 (1989)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.