Simulation of Counter-Current Spontaneous Imbibition Based on Momentum Equations with Viscous Coupling, Brinkman Terms and Compressible Fluids

Abstract

A numerical model is investigated representing counter-current spontaneous imbibition of water to displace oil or gas from a core plug. The model is based on mass and momentum conservation equations in the framework of the theory of mixtures. We extend a previous imbibition model that included fluid–rock friction and fluid–fluid drag interaction (viscous coupling) by including fluid compressibility and Brinkman viscous terms. Gas compressibility accelerated recovery due to gas expansion from high initial non-wetting pressure to ambient pressure at typical lab conditions. Gas compressibility gave a recovery profile with two characteristic linear sections against square root of time which could match tight rock literature experiments. Brinkman terms decelerated recovery and delayed onset of imbibition. Experiments where this was prominent were successfully matched. Both compressibility and Brinkman terms caused recovery deviation from classical linearity with the square root of time. Scaling yielded dimensionless numbers when Brinkman term effects were significant.

Article Highlights

• Spontaneous imbibition with viscous coupling, compressibility and Brinkman terms.

• Viscous coupling reduces spontaneous imbibition rate by fluid–fluid friction.

• Brinkman terms delay early recovery and explain seen delayed onset of imbibition.

• Gas compressibility accelerates recovery and can be significant at lab conditions.

• Gas compressibility gives recovery with two root of time lines as seen for shale.

1 Introduction

Naturally fractured reservoirs contribute to around 20% of the hydrocarbon reserves discovered worldwide. These reservoirs are usually considered to have dual porosity as given by disconnected matrix blocks and connected fractures. The fractures have high permeability compared to the matrix, but much less volume (Warren and Root 1963). Oil or gas recovery depends on matrix transfer which can occur through fluid expansion, gravity drainage and capillary imbibition. Water injection has been applied successfully in naturally fractured reservoirs, but is most effective when the matrix is water-wet and capillary forces can take up the injected water spontaneously (from the fractures). This process is referred as spontaneous imbibition and can recover as much oil or gas at matrix level as by forced imbibition if the matrix is strongly water-wet, but less under other wetting conditions (Zhou et al. 2000). Imbibition of fluids carrying wettability alternating chemicals can further improve recovery (Zhang et al. 2006; Mamonov et al. 2019; Andersen and Ahmed 2021). Spontaneous imbibition can occur counter-currently, where the wetting and non-wetting fluids flow in opposite direction. This usually happens when all open sides of the matrix are exposed to wetting phase and gravity is negligible (Morrow and Mason 2001). If wetting and non-wetting phases cover different parts of the matrix, the non-wetting phase will be produced predominantly out of the sides covered by non-wetting phase, while wetting phase enters at the sides open to wetting phase (Bourbiaux and Kalaydjian 1990; Andersen et al. 2019a, b). Spontaneous imbibition is regarded as a crucial driving mechanism for oil recovery from naturally fractured reservoirs (Morrow and Mason 2001; Andersen et al. 2014; Abd et al. 2019). Numerous works have modeled this phenomenon with analytical and numerical approaches (Mattax and Kyte 1962; Ma et al. 1997; Mason et al. 2012; Schmid and Geiger 2012).

An established assumption in modeling of single and multiphase flow in porous media is the Darcy equation (Darcy 1856) and its extension with relative permeabilities. This equation states that the flux of a fluid is proportional to only its own pressure gradient and accounts for the presence of other fluids via the saturation dependent relative permeability factor. However, there are restrictions on the validity of the Darcy equation for modeling some porous medium flows; that is, in closely packed media, saturated fluid flows at slow velocity but with relatively large Reynolds numbers. The flows in such closely packed medium behave nonlinearly and cannot be modeled accurately by the linear Darcy equation (Skrzypacz and Wei 2017).

A more general approach is the theory of mixtures where both the solid and fluid materials are considered continua and each spatial point can be occupied by a fraction of fluid and solid particles (Munaf et al. 1993). The flux relations of each phase are described by momentum equations accounting for more mechanisms. A full hierarchy of flow models for porous media can be derived, with Darcy’s equation as a special case. Different assumptions of which mechanisms are significant lead to equations of flow through solids which can be permitted to interact with the flow process. Examples of such interactions include deformation, swelling due to fluid adsorption, diffusion of solvents through polymeric solids and diffusion of biological fluids through biological solids (Rajagopal 2007).

Some advantages of generalizing the conventional Darcy’s law are the ability to capture experimental observations which may contradict its assumptions or predictions. Some works have argued that during multiphase flow the flux of a phase depends on both its own pressure gradient and that of the other phase as well. The importance of these so-called cross-mobilities has been studied in several works (Bourbiaux and Kalaydjian 1990; Bentsen and Manai 1993; Avraam and Payatakes 1995; Standnes et al. 2017) and can explain why measured relative permeabilities should be reduced during counter-current flow compared to co-current flow. This phenomenon has been studied in the context of spontaneous and forced imbibition by Qiao et al. (2018), Andersen et al. (2019a, 2020) and demonstrates that such mechanisms should be accounted for in order to not predict too optimistic imbibition rates. Mathematical analysis on the behavior of two phase flow equations with viscous coupling and compressibility and associated numerical schemes were studied by Qiao et al. (2019b).

In addition, it is well-known that, there is difficulty when applying Darcy’s law for a viscous fluid (Deng and Martinez 2005), especially when the internally frictional resistance in the fluid is greater than the frictional force between the fluid and the solid surface at the porous boundaries of the solid. An easy way to resolve this difficulty is to modify the Darcy equation by including a second-order viscous term. Brinkman firstly proposed this modification (Brinkman, 1949) and the corresponding equation is called as the Brinkman–Darcy equation. There are many investigations related to this formulation. For example, Coclite et al. (2014) mathematically analyzed the Brinkman regularization of the two-phase flow equations and proved existence of weak solutions for such equations. Qiao et al. (2019a) considered the effects of both Brinkman terms and fluid compressibility in a two-phase flow model and found that the injected fluid through a horizontal core has a slow displacement speed when it is compressible and that a high front saturation of injected fluid can be formed when the fluids have an extremely large viscous effect or Brinkman terms. Varsakelis and Papalexandris (2020) dealt with the derivation of tortuosity estimates based on the Darcy–Brinkman for a polymer-filled system.

In this work we consider a system with 1D counter-current spontaneous imbibition of two immiscible fluids, water (wetting fluid) and oil or gas (non-wetting fluid). The flux relations account for viscous coupling and Brinkman terms and fluid compressibility is included. We study the influence of these mechanisms on recovery of non-wetting phase, fluid distributions and deviations from behavior under standard assumptions (relative permeability formulations). Some research questions we wish to address are:

• Can compressibility influence spontaneous imbibition experiments at typical lab conditions and what is their impact?

• Can viscous coupling, compressibility or Brinkman terms explain the delayed onset of imbibition sometimes seen experimentally or lead to changes in the shape of the recovery profile?

Regarding the first question we note that many works measure spontaneous imbibition of water displacing air, but model the system assuming the fluids are incompressible (Akin et al. 2000; Rangel-German and Kovscek 2002; Li and Horne 2006; Le Guen and Kovscek 2006). Many of these experiments are conducted at ambient conditions, but in systems with high capillary pressure which can allow significant gas expansion. Andersen et al. (2019b) noted a better fit of co-current imbibition experiments accounting for gas compressibility. For the second question we note that counter-current spontaneous imbibition in 1D linear systems theoretically should follow a square root of time profile (McWhorter and Sunada 1990). Several imbibition experiments, however, show that the expected flow regime is not established until after a certain delay time (Akin et al. 2000; Føyen et al. 2019). Unconventional porous media can also deviate from the square root of time profile and be better described by other exponents.

The model based on mixture theory describing 1D spontaneous imbibition with relevant mechanisms and boundary conditions is presented in Sect. 2. Results follow in Sect. 3 with sensitivity analysis and interpretation of how the mechanisms fit into experimental setups and observations. Finally, we summarize the paper with conclusions in Sect. 4.

2 Theory

2.1 Theory of Mixtures—Key Points

The model utilized in the present paper is based on the theory of mixtures and a general formulation is hence first presented. In the theory of mixtures, both solid ($s$) and fluid ($f$) materials are idealized as continua. A fundamental assumption is that each spatial point can be occupied by particles of both fluid and solid continua that, regarded as resulting from a local averaging process (Munaf et al. 1993). The mass conservation equation for each species is given by:

$$\left( {\rho_{i} } \right)_{t} + \nabla \cdot \left( {\rho_{i} u_{i} } \right) = q_{i} , \left( {i = s,f} \right)$$

(1)

where ${\rho }_{i}$ is the bulk density, ${u}_{i}$ is the velocity vector and ${q}_{i}$ is the mass source term. The equation of momentum balance for solid and fluid has the form:

$$\nabla \cdot \left( {{\mathbf{T}}_{i} } \right) + {\Pi }_{i} + \rho_{i} b_{i} = \rho_{i} a_{i} ,{ }\left( {i = s,f} \right)$$

(2)

where ${\mathbf{T}}_{i}$ is the partial stress tensor, ${\Pi }_{i}$ interaction body force, ${b}_{i}$ external body force and ${a}_{i}$ acceleration. Note that this is a vector equation to account for momentum balance in 3 directions. Studies show that each partial stress ${\mathbf{T}}^{s}$, ${\mathbf{T}}^{f}$ can depend on both densities ${\rho }_{s}$, ${\rho }_{f}$, deformation gradient of the solid that is elastic or swells due to fluid absorption, etc., and the rate of deformation of the fluid. This allows the effects of viscosity, e.g., shear thickening, shear thinning (Munaf et al. 1993; Rajagopal 2007). The interaction body force ${\Pi }_{i}$ can depend on the densities and their gradients, the solid deformation measure and its gradient, the fluid–solid relative velocity and a measure of a relative acceleration. This allows for effects due to difference in particles acceleration of the constituents (Munaf et al. 1993).

2.2 System and Geometry

We consider two immiscible fluids (wetting and non-wetting) that flow inside porous rock in 1D linear horizontal direction due to spontaneous imbibition, see Fig. 1. The wetting fluid ($w$) is water, while the non-wetting ($nw$) fluid is either oil or gas. Oil, water and the matrix are assumed incompressible, while gas is assumed compressible. The left end of the matrix $x=0$ is open, while the right end $x=L$ is closed. The matrix is initially filled with non-wetting phase.

Fig. 1

Schematic representation of the counter-current spontaneous imbibition system

2.3 Formulation of the Two-Phase Flow Model Based on the Theory of Mixtures

We now formulate our flow model based on the above formulations of mixture theory. The following assumptions are made:

• (i) The solid is a rigid porous body and thus the mass conservation and momentum balance equations of the solid can be ignored.

• (ii) Two immiscible fluids $i=w,nw$ each occupies part of the porous space.

• (iii) The interaction forces that come into play include the frictional forces that fluids encounter at the boundaries of the pores as well as the viscous coupling forces that one fluid exerts on another. This can be captured by a “drag like” term that is proportional to the difference in the velocity between the two constituents. The drag coefficients being a constant.

• (iv) The frictional effects within the fluid due to its viscosity are reflected in the expression for the partial stress tensor for the fluid phase.

• (v) The flow is sufficiently slow that inertial nonlinearities can be neglected.

• (vi) All fluids can be modeled compressible, but in our examples only gas is compressible, while oil and water are incompressible.

2.3.1 Transport Equations

Since motion and compressibility of the solid are ignored, only the mass and momentum conservation equations for the fluids are considered. We replace the subscript $f$ with the relevant phase considered ($w,nw$). The pore space is fully occupied by fluids; therefore, the total saturation equals unity:

$$s_{w} + s_{nw} = 1$$

(3)

where $s_{w}$, $s_{nw}$ are wetting and non-wetting fluid saturations, respectively. We will have two mass balance equations, one for each fluid:

$$\left( {\phi m} \right)_{t} { } + { }\left( {\phi mu_{w} } \right)_{x} { } = { }0$$

(4)

$$\left( {\phi n} \right)_{t} { } + { }\left( {\phi nu_{nw} } \right)_{x} { } = { }0$$

(5)

Source terms are not included ($q_{i} = 0$). $\phi$ is the rock porosity, and $m$ and $n$ are the mass per pore volume of wetting and non-wetting fluids, respectively, given by:

$$m = s_{w} \rho_{w} ,{ }n = s_{nw} \rho_{nw} .$$

(6)

$u_{w}$ and $u_{nw}$ are wetting and non-wetting fluids interstitial velocities. Compared to the general theory of mixtures, we have made the following transformation:

$$\rho_{i} \to \phi s_{i} \rho_{i} u_{i} ,{ }\left( {i = w,nw} \right),$$

(7)

where the phase density $\rho_{i}$ now is per volume phase rather than per bulk volume. Also, Qiao et al. (2018) introduced an effective porosity $\phi_{e}$ which combines rock porosity and the phase residual saturations $s_{wr} ,s_{nwr}$ (where the respective phases do not flow) as follows:

$$\phi_{e} = \phi \left( {1 - s_{wr} - s_{nwr} } \right)$$

(8)

We can then make the transformations:

$$\phi \to \phi_{e} ,{ }s_{w} \to S_{w} = \frac{{s_{w} - s_{wr} }}{{1 - s_{nwr} - s_{wr} }}$$

(9)

(where ${S}_{w}$ is normalized water saturation between 0 and 1) to account for residual saturations in a way consistent with traditional modeling of porous media flow. In the following we will thus consider effective porosities and normalized saturations.

2.3.2 Momentum Equations

For the momentum balance equations, see (2), the inertia of the fluids is neglected, and we assume no external body forces (such as gravity) and negligible acceleration of fluids:

$$a_{i} = b_{i} = 0, \left( {i = nw,w} \right).$$

(10)

Equation (2) reduces to:

$$\nabla \cdot \left( {{\mathbf{T}}_{i} } \right) + {\Pi }_{i} = 0, \left( {i = nw,w} \right).$$

(11)

The interaction body force ${\Pi }_{i}$ depends on the relative velocity between solid and fluid and on the relative velocity between fluids, which in 1D can be expressed as (Qiao et al. 2018; Qiao and Evje 2020):

$${\Pi }_{i} = P_{i} \nabla S_{i} - \hat{k}_{i} u_{i} - \hat{k}\left( {u_{i} - u_{j} } \right),{ }\left( {i,j = nw,w;i \ne j} \right),$$

(12)

where $\hat{k}_{i}$ and $\hat{k}$ are the fluid–solid and fluid–fluid interaction coefficients respectively. The term $P_{i} \nabla S_{i}$ is related to interfacial force imposed by other phase on phase $i$, arising from averaging the mathematical equations (Drew and Segel 1971). The partial stress has the form:

$${\mathbf{T}}_{i} = - S_{i} P_{i} {\mathbf{I}} + \varepsilon_{i} S_{i} \rho_{i} \nabla u_{i} , \left( {i = nw,w} \right).$$

(13)

${p}_{i}$ is the fluid pressure, $\mathbf{I}$ is the identity tensor, and ${\varepsilon }_{i}$ is the coefficient of the Brinkman term and has unit m²/s. Taking (11) and (13) to 1D and combining with (12), we obtain the momentum balance equations for both fluids:

$$S_{w} \left( {P_{w} } \right)_{x} = - \hat{k}_{w} u_{w} - \hat{k}\left( {u_{w} - u_{nw} } \right) + \varepsilon_{w} \left( {m\left( {u_{w} } \right)_{x} } \right)_{x} ,$$

(14)

$$S_{nw} \left( {P_{nw} } \right)_{x} = - \hat{k}_{nw} u_{nw} + \hat{k}\left( {u_{w} - u_{nw} } \right) + \varepsilon_{nw} \left( {n\left( {u_{nw} } \right)_{x} } \right)_{x} .$$

(15)

Thus, we get modified Darcy–Brinkman equations which are able to account for fluid–fluid interaction (viscous coupling). The Darcy–Brinkman equation is a governing equation for flow through a porous medium with an extra Laplacian (viscous) term (Brinkman term) added to the classical Darcy equation. The equation has been widely used to analyze flows in highly porous media (Deng and Martinez 2005). There has been some literature discussion on what values are representative of the Brinkman coefficients. Although in bulk phase it can be compared to a kinematic fluid viscosity, research suggests that it can be greater than that in porous media flow. Valdes-Parada et al. (2007) showed that the Brinkman coefficient, being greater than the fluid viscosity, should be a decreasing function of the porous medium porosity. Several authors have also pointed out the deviation from a kinematic viscosity value and we refer to Kim and Russel (1985), Kolodziej (1988) and Martys et al. (1994).

2.3.3 Final Set of Equations

The capillary pressure is defined as a difference between the non-wetting phase pressure ${P}_{nw}$ and wetting phase pressure ${P}_{w}$.

$$P_{c} = { }P_{nw} {-}P_{w} .$$

(16)

For fixed rock and fluid properties, capillary pressure is considered a function of saturation only as long as the saturation changes monotonously. Variations in rock and fluid properties are accounted for by assuming an invariant $J$-function (Bear 2013):

$$P_{c} = \sigma \sqrt {\frac{\phi }{K}} { }J\left( {S_{w} } \right),$$

(17)

where $\sigma$ is oil–water or water–gas interfacial tension. The fluid compressibility is accounted for by letting the densities to vary according to the relations (Qiao et al. 2019a; Qiao and Evje 2020):

$$\rho_{w} - \tilde{\rho }_{w0} = \frac{{P_{w} }}{{C_{w} }},$$

(18)

$$\rho_{nw} - \tilde{\rho }_{nw0} = \frac{{P_{nw} }}{{C_{nw} }}.$$

(19)

For weakly compressible fluids (such as liquids), $\frac{{P_{i} }}{{C_{i} }} < \tilde{\rho }_{i0}$ and $\tilde{\rho }_{i0}$ correspond to the density at a low pressure. $C_{i}$ can then be considered the inverse compressibility. Incompressible fluids are obtained by letting $C_{i} \to \infty$. For highly compressible fluids (such as gas) we have $\frac{{P_{i} }}{{C_{i} }} > \tilde{\rho }_{i0}$. In fact, setting $\tilde{\rho }_{i0} = 0$ results in the ideal gas law for the appropriate choice of $C_{i}$. In our examples involving gas, we set ${C}_{nw}={10}^{5} \mathrm{Pa}\frac{{m}^{3}}{\mathrm{kg}}$ to obtain a gas density of 1 kg/m³ at standard conditions. For direct comparison with an incompressible case we will set ${\stackrel{\sim }{\rho }}_{i0}=1$ kg/m³ and $C_{nw} = 10^{10} {\text{Pa}}\frac{{m^{3} }}{{{\text{kg}}}}$.

We make use of the two mass equations, (4) and (5), and the above-defined relations, (16), (18) and (19), to eliminate non-wetting phase pressure. A wetting phase pressure evolution equation can be obtained by summing up the two mass equations, (4) and (5), after pre-multiplying them with ${\rho }_{nw}$ and ${\rho }_{w}$, respectively:

$$(P_{w} )_{t} + \tilde{\eta }\rho_{w} \left( {nu_{nw} } \right)_{x} + \tilde{\eta }\tilde{\rho }_{nw} \left( {mu_{w} } \right)_{x} = 0$$

(20)

where we have introduced the dynamic coefficients:

$$\tilde{\eta } = \frac{{C_{w} C_{nw} }}{{S_{nw} \rho_{w} C_{w} + S_{w} \tilde{\rho }_{nw} C_{nw} }}$$

(21)

$$\tilde{\rho }_{nw} = \rho_{nw} - \frac{{S_{nw} P_{c}^{^{\prime}} }}{{C_{nw} }}$$

(22)

To summarize, the model consists of 4 equations that should be solved:

• 1 mass conservation Eq. (23),

• 1 pressure evolution Eq. (24),

• 2 momentum balance Eqs. (25) and (26) which include viscous coupling terms and Brinkman terms.

  $$\left( m \right)_{t} + \left( {mu_{w} } \right)_{x} = 0,$$

  (23)
  $$(P_{w} )_{t} + \tilde{\eta }\rho_{w} \left( {nu_{nw} } \right)_{x} + \tilde{\eta }\tilde{\rho }_{nw} \left( {mu_{w} } \right)_{x} = 0,$$

  (24)
  $$S_{nw} (P_{w} + P_{c} )_{x} = - \hat{k}_{nw} u_{nw} + \hat{k}\left( {u_{w} - u_{nw} } \right) + \varepsilon_{nw} \left( {n\left( {u_{nw} } \right)_{x} } \right)_{x} ,$$

  (25)
  $$S_{w} (P_{w} )_{x} = - \hat{k}_{w} u_{w} + \hat{k}\left( {u_{nw} - u_{w} } \right) + \varepsilon_{w} \left( {m\left( {u_{w} } \right)_{x} } \right)_{x} .$$

  (26)

2.3.4 Alternative Formulation of the Water Mass BALANCE

In this section we make an alternative representation of the water mass balance (23) using capillary pressure, fractional flow function and total Darcy velocity. We may express the interstitial phase velocities ${u}_{w}$, ${u}_{nw}$ by solving the two momentum Eqs. (14) and (15)

$$\begin{aligned} u_{w} = & - \frac{{S_{w} \left( {\hat{k}_{nw} + \hat{k}} \right)}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}(P_{w} )_{x} \\ & \quad - \frac{{S_{nw} \hat{k}}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}\left( {P_{nw} } \right)_{x} \\ & \quad \quad + \varepsilon_{w} \frac{{\hat{k}_{nw} + \hat{k}}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}\left( {m\left( {u_{w} } \right)_{x} } \right)_{x} \\ & \quad \quad + \varepsilon_{nw} \frac{{\hat{k}}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}\left( {n\left( {u_{nw} } \right)_{x} } \right)_{x} , \\ \end{aligned}$$

(27)

$$\begin{aligned} u_{nw} = & - \frac{{S_{nw} \left( {\hat{k} + \hat{k}_{w} } \right)}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}(P_{nw} )_{x} \\ & \quad - \frac{{S_{w} \hat{k}}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}(P_{w} )_{x} \\ & \quad \quad + \varepsilon_{nw} \frac{{\hat{k}_{w} + \hat{k}}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}\left( {n\left( {u_{nw} } \right)_{x} } \right)_{x} \\ & \quad \quad + \varepsilon_{w} \frac{{\hat{k}}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}\left( {m\left( {u_{w} } \right)_{x} } \right)_{x} \\ \end{aligned}$$

(28)

By introducing the notations of cross-term mobilities:

$$\hat{\lambda }_{ww} = \frac{{S_{w}^{2} \left( {\hat{k}_{nw} + \hat{k}} \right)}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}\phi_{e}$$

(29)

$$\hat{\lambda }_{nn} = \frac{{S_{nw}^{2} \left( {\hat{k}_{w} + \hat{k}} \right)}}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}\phi_{e}$$

(30)

$$\hat{\lambda } = \frac{{\hat{k}S_{w} S_{nw} }}{{\hat{k}_{nw} \hat{k}_{w} + \hat{k}\left( {\hat{k}_{nw} + \hat{k}_{w} } \right)}}\phi_{e}$$

(31)

the Darcy phase velocities ${U}_{i}={u}_{i}{S}_{i}{\phi }_{e}, \left(i=nw,w\right)$ are given as:

$$U_{w} = - \hat{\lambda }_{ww} (P_{w} )_{x} - \hat{\lambda }(P_{nw} )_{x} + \varepsilon_{w} \frac{{\hat{\lambda }_{ww} }}{{S_{w} }}\left( {m\left( {u_{w} } \right)_{x} } \right)_{x} + \varepsilon_{nw} \frac{{\hat{\lambda }}}{{S_{nw} }}\left( {n\left( {u_{nw} } \right)_{x} } \right)_{x}$$

(32)

$$U_{nw} = - \hat{\lambda }_{nn} (P_{nw} )_{x} - \hat{\lambda }(P_{w} )_{x} + \varepsilon_{nw} \frac{{\hat{\lambda }_{nn} }}{{S_{nw} }}\left( {n\left( {u_{nw} } \right)_{x} } \right)_{x} + \varepsilon_{w} \frac{{\hat{\lambda }}}{{S_{w} }}\left( {m\left( {u_{w} } \right)_{x} } \right)_{x}$$

(33)

Based on total Darcy phase velocity $U_{T} = U_{nw} + U_{w}$ and phase pressure difference defined in (16), we have the water phase Darcy velocity in terms of advection, diffusion and Brinkman effects:

$$U_{w} = \hat{f}_{w} U_{T} + W\left( {P_{c} } \right)_{x} - W\frac{{\varepsilon_{nw} }}{{S_{nw} }}\left( {n\left( {u_{nw} } \right)_{x} } \right)_{x} - \frac{{\varepsilon_{w} }}{{S_{w} }}W_{w} \left( {m\left( {u_{w} } \right)_{x} } \right)_{x}$$

(34)

where the following relations are used:

$$\hat{\lambda }_{w} = \hat{\lambda }_{ww} + \hat{\lambda }, \;\hat{\lambda }_{nw} = \hat{\lambda }_{nn} + \hat{\lambda },$$

(35)

$$\hat{\lambda }_{T} = \hat{\lambda }_{nw} + \hat{\lambda }_{w} ,\quad \hat{f}_{w} = \frac{{\hat{\lambda }_{w} }}{{\hat{\lambda }_{T} }}, \quad W = \hat{f}_{w} \hat{\lambda }_{nw} - \hat{\lambda }, W_{w} = \hat{f}_{w} \hat{\lambda }_{w} - \hat{\lambda }_{ww} .$$

(36)

Here ${\widehat{\lambda }}_{w}$, ${\widehat{\lambda }}_{nw}$ and ${\widehat{\lambda }}_{T}$ are called generalized mobilities for water and non-wetting phase and total generalized mobility. ${\widehat{f}}_{w}$ is the water fractional flow function and $W$ and ${W}_{w}$ are mobility coefficients. The water transport equation is then obtained by combining (4) and (34) and using ${\left({\phi }_{e}m{u}_{w}\right)}_{x}={\left({{\rho }_{w}U}_{w}\right)}_{x}$ and $m={S}_{w}{\rho }_{w}$ such that:

$$\left( {\phi_{e} S_{w} \rho_{w} } \right)_{t} = - \left( {\rho_{w} \hat{f}_{w} U_{T} } \right)_{x} - \left[ {\rho_{w} W(P_{c} )_{x} } \right]_{x} - \left( {\rho_{w} W\frac{{\varepsilon_{nw} }}{{S_{nw} }}\left( {n\left( {u_{nw} } \right)_{x} } \right)_{x} + \rho_{w} \frac{{\varepsilon_{w} }}{{S_{w} }}W_{w} \left( {m\left( {u_{w} } \right)_{x} } \right)_{x} } \right)_{x} .$$

(37)

We can observe from (37) that the water mass evolves with changes from an advection term, capillary diffusion term as well as the Brinkman terms. This expression will be useful for understanding the mechanisms driving the transport.

2.4 Boundary Conditions

On the left boundary of the domain $x=0$ we assume a zero capillary pressure and water phase pressure equal to atmospheric pressure, 1 bar:

$$P_{c} \left. \right|_{x = 0} = 0, P_{w} \left. \right|_{x = 0} = 10^{5} {\text{Pa}}.$$

(38)

Boundary values of saturation and densities follow directly from the fixed capillary pressure and absolute pressure in (38):

$$S_{w} |_{{{\text{x}} = 0}} = P_{c}^{ - 1} \left( {P_{c} \left. \right|_{x = 0} } \right),{ }\rho_{w} \left. \right|_{x = 0} = \rho_{w} \left( {P_{w} \left. \right|_{x = 0} } \right),{ }\rho_{nw} \left. \right|_{x = 0} = \rho_{nw} \left( {P_{nw} \left. \right|_{x = 0} } \right),$$

(39)

$$m|_{x = 0} = \rho_{w} |_{x = 0} \cdot S_{w} |_{x = 0} ,{ }n|_{x = 0} { } = { }\rho_{nw} |_{x = 0} \cdot S_{nw} |_{x = 0} .$$

(40)

We assume $m\left( {u_{w} } \right)_{x}$ and $n\left( {u_{nw} } \right)_{x}$ in the viscous terms on the left open boundary equal to 0 in the momentum Eqs. (25) and (26):

$$m\left( {u_{w} } \right)_{x} |_{x = 0} = n\left( {u_{nw} } \right)_{x} |_{x = 0} = 0.$$

(41)

On the right (closed) boundary a zero-flux condition is applied for both phases:

$$u_{w} |_{x = L} = u_{nw} |_{x = L} = 0.$$

(42)

2.5 Initial Conditions

The initial state is given by full saturation of non-wetting phase:

$$S_{w} \left( x \right)\left. \right|_{t = 0} = 0$$

(43)

The initial pressure of the wetting phase is that at the boundary, while that of the non-wetting phase is that of the wetting phase plus the initial capillary pressure:

$$P_{{\text{w}}} \left( x \right)\left. \right|_{t = 0} = 10^{5} {\text{ Pa}},$$

(44)

That also determines the initial phase densities:

$${ }\rho_{{\text{w}}} \left( x \right)\left. \right|_{t = 0} = \rho_{w} \left( {10^{5} {\text{ Pa}}} \right),{ }\rho_{{{\text{nw}}}} \left( x \right)\left. \right|_{t = 0} = \rho_{nw} \left( {P_{nw} = 10^{5} {\text{ Pa}} + P_{c} \left( {S_{w} = 0} \right)} \right)$$

(45)

2.6 Functional Relations

The dimensionless capillary pressure function $J\left({S}_{w}\right)$ takes the form (Andersen et al. 2017):

$$J\left( {S_{w} } \right) = \frac{{J_{1} }}{{\left( {1 + k_{1} S_{w} } \right)^{{n_{1} }} }} - \frac{{J_{2} }}{{\left( {1 + k_{2} \left( {1 - S_{w} } \right)} \right)^{{n_{2} }} }} + J_{3}$$

(46)

${J}_{1},{J}_{2},{J}_{3},{k}_{1},{k}_{2},{n}_{1},{n}_{2}$ are nonnegative curve fitting parameters. The fluid–solid and fluid–fluid interaction coefficients are defined as (Standnes et al. 2017; Qiao et al. 2018):

$$\hat{k}_{w} = I_{w} \frac{{\mu_{w} }}{K}\phi_{e} S_{w}^{\alpha }$$

(47)

$$\hat{k}_{nw} = I_{nw} \frac{{\mu_{nw} }}{K}\phi_{e} \left( {1 - S_{w} } \right)^{\beta }$$

(48)

$$\hat{k} = I\frac{{\mu_{w} \mu_{nw} }}{K}\phi_{e} S_{w} \left( {1 - S_{w} } \right)$$

(49)

Here ${I}_{w}$, ${I}_{nw}$ and $I$ are the wetting fluid–solid, non-wetting fluid–solid and wetting–non-wetting fluids interaction parameters that characterize the strength of resistance forces. $\alpha$ and $\beta$ are exponents that characterize how the interaction terms depend on saturation.

2.7 Co- and Counter-Current Relative Permeabilities

Under the assumption of incompressible co-current flow (both phases have equal pressure gradients) and negligible viscous terms, it was shown in Qiao et al. (2018) that co-current relative permeabilities could be formulated as:

$$k_{rw}^{co} = \frac{{S_{w}^{2 - \alpha } \left( {I_{nw} + IS_{nw}^{1 - \beta } \mu_{w} } \right)}}{{I_{nw} I_{w} + I\left( {I_{nw} S_{nw} S_{w}^{1 - \alpha } \mu_{nw} + I_{w} S_{nw}^{1 - \beta } S_{w} \mu_{w} } \right)}}$$

(50)

$$k_{rnw}^{co} = \frac{{S_{nw}^{2 - \beta } \left( {I_{w} + IS_{w}^{1 - \alpha } \mu_{nw} } \right)}}{{I_{nw} I_{w} + I\left( {I_{nw} S_{nw} S_{w}^{1 - \alpha } \mu_{nw} + I_{w} S_{nw}^{1 - \beta } S_{w} \mu_{w} } \right)}}$$

(51)

where the Eqs. (47), (48) and (49) have been used and we omit the details of derivation here. It was later shown in Andersen et al. (2020) that for incompressible counter-current flow (where the Darcy fluxes are opposite) without viscous terms, counter-current relative permeabilities were related to the co-current relative permeabilities by a saturation dependent coefficient $C$ (which is the same for both phases):

$$C = \frac{{I_{nw} I_{w} + I\left( {I_{nw} S_{nw} S_{w}^{1 - \alpha } \mu_{nw} + I_{w} S_{nw}^{1 - \beta } S_{w} \mu_{w} } \right)}}{{\left( {I_{nw} + IS_{nw}^{1 - \beta } \mu_{w} } \right)\left( {I_{w} + IS_{w}^{1 - \alpha } \mu_{nw} } \right)}}$$

(52)

Taking the co-current relative permeability expressions (50) or (51) factored by $C$ in (52) results in the following counter-current relative permeability expressions:

$$k_{rw}^{cou} = \frac{{S_{w}^{2 - \alpha } }}{{I_{w} + IS_{w}^{1 - \alpha } \mu_{nw} }}, k_{rnw}^{cou} = \frac{{S_{nw}^{2 - \beta } }}{{I_{nw} + IS_{nw}^{1 - \beta } \mu_{w} }}$$

(53)

which can be implemented in a standard porous media simulator for comparison with our momentum equation system under the counter-current flow setup. For the special case where there is no fluid–fluid interaction, the co- and counter-current expressions become identical and equivalent to Corey relative permeabilities:

$$k_{rw} = \frac{{S_{w}^{2 - \alpha } }}{{I_{w} }}, k_{rnw} = \frac{{S_{nw}^{2 - \beta } }}{{I_{nw} }}$$

(54)

where the end points are $1/{I}_{i}$ and the Corey exponents are $2-\alpha$ and $2-\beta$ for wetting and non-wetting phase, respectively.

2.8 Discretization and Grid

The main equations are solved by the finite difference method. The mass balance Eq. (23) is solved by explicit upwind scheme. The system consisting of pressure evolution (24) and two momentum balance Eqs. (25) and (26) is solved implicitly. The physical domain is discretized with uniform step size. Staggered grid is used for discretization of the equations with ${S}_{w}, {P}_{w}, {P}_{c}$ in the center of the nodes and ${u}_{w}, {u}_{o}$ in the cell interfaces. To achieve numerical stability, upwind scheme and a reasonable time interval are utilized, see “Appendix A”. 200 grid cells were used for the simulations, which was found to give sufficient convergence for our purposes. See “Appendix B” for a grid sensitivity analysis and “Appendix C” for analysis of the global mass error evolution of the numerical scheme.@@@

2.9 Recovery Factor

The recovery factor $RF$ is computed as the volume fraction of imbibed wetting phase to displaceable initial volume of non-wetting fluid. This reduces to the average normalized wetting saturation:

$$RF = \frac{1}{L}\mathop \smallint \limits_{0}^{L} S_{w} dx$$

(55)

Given sufficient time, $RF$ will reach ${S}_{weq}$ (the normalized saturation where ${P}_{c}=0$), which is equal to 1 for strongly water-wet media and less than one otherwise.

2.10 Simulation Result Format

We illustrate our results inspired by classical spontaneous imbibition behavior. McWhorter and Sunada (1990) derived exact integral solutions for 1D linear spontaneous imbibition flow of two incompressible fluids. They showed that cumulative water imbibed ${Q}_{w}$ (from a unit area), and equivalently volume oil produced from that area, was given by a constant inflow parameter $A$ and the square root of time:

$$Q_{w} \left( t \right) = \mathop \smallint \limits_{0}^{t} q_{w0} \left( t \right)dt = { }2A\sqrt t$$

(56)

where ${q}_{w0}$ is the water volumetric flux at the inlet. Self-similar solutions were obtained in the sense that the position of each saturation was proportional to the square root of time and the saturation derivative of a function $F$; hence, profiles at different times could be compressed or expanded to completely overlap. The linear dependence with square root of time for position and recovery is valid until the no-flow boundary is reached.

Based on how some parameters appear linearly in the diffusion equation for a 1D system describing linear spontaneous imbibition under standard mechanisms (when using relative permeabilities and incompressible fluids); recovery and saturation profiles appear identically when compared at identical scaled times ${t}_{D}$:

$$t_{D} = \frac{t}{\tau },{ }\tau = \frac{{L^{2} \mu_{ref} }}{\sigma }\sqrt {\frac{\phi }{K}}$$

(57)

where ${\mu }_{ref}$ is a representative viscosity either based on the imbibing fluid or the mean of both fluids. We will use ${\mu }_{ref}={\mu }_{w}$ resulting in a similar scaling as Mattax and Kyte (1962). Scaling with (57) assumes that the capillary diffusion coefficient does not change with variation of the considered parameters; hence, changes of viscosity ratio or saturation functions violates the scaling, but still gives square root of time behavior. Andersen et al. (2020) showed that a model for spontaneous imbibition in terms of momentum equations with viscous coupling could be expressed with effective relative permeabilities and therefore still result in square root of time recovery. Standnes and Andersen (2017) used viscous coupling as a parameter to better explain time scales of recovery during variation of fluid viscosities. Ma et al. (1997) proposed how to include variation in geometry in the scaling by replacing $L$ with a characteristic length ${L}_{c}$, but the change in geometry from linear causes recovery to generally not follow a root of time trend. We will present:

• Recovery factor $RF$ against the square root of scaled or unscaled time

• Saturation profiles at identical scaled or unscaled times.

Scaling recovery allows us to distinguish nonstandard mechanisms by recovery following other trends than the square root of time and saturation profiles to not overlap at identical scaled times.

2.10.1 Scaling Numbers for Relative Importance of Brinkman Terms

We here aim to understand when the Brinkman terms affect the spontaneous imbibition process. We introduce the following scaling:

$$\begin{gathered} x_{D} = \frac{x}{L}, t_{D} = \frac{t}{\tau }, u_{iD} = \frac{{u_{i} }}{{\left( {\frac{L}{{\phi_{e} \tau }}} \right)}}, U_{iD} = \frac{{U_{i} }}{{\left( {\frac{L}{\tau }} \right)}}, \hfill \\ m_{D} = \frac{m}{{\rho_{w} }}, n_{D} = \frac{n}{{\rho_{nw0} }}, W_{D} = \frac{W}{{\left( {\frac{K}{{\mu_{ref} }}} \right)}}, W_{wD} = \frac{{W_{w} }}{{\left( {\frac{K}{{\mu_{ref} }}} \right)}} \hfill \\ \end{gathered}$$

(58)

and use our assumption of constant water density in (37) which yields the following scaled water mass balance equation:

$$\begin{aligned} \phi_{e} \left( {S_{w} } \right)_{T} & = - \left( {\hat{f}_{w} U_{TD} } \right)_{X} - \phi \left[ {W_{D} \left( J \right)_{X} } \right]_{X} \\ & \quad - \left( {\frac{{\varepsilon_{nw} K\rho_{nw0} }}{{\mu_{ref} \phi_{e} L^{2} }}\frac{{W_{D} }}{{S_{o} }}\left( {n_{D} \left( {u_{nwD} } \right)_{X} } \right)_{X} + \frac{{\varepsilon_{w} K\rho_{w} }}{{\mu_{ref} \phi_{e} L^{2} }}\frac{{W_{wD} }}{{S_{w} }}\left( {m_{D} \left( {u_{wD} } \right)_{X} } \right)_{X} } \right)_{X} . \\ \end{aligned}$$

(59)

The first term on the RHS may be nonzero due to compressibility (in gas cases) but will be zero otherwise. The capillary diffusion term is considered the main driving force during counter-current spontaneous imbibition. By dividing the fluxes in the Brinkman terms to the flux in the capillary diffusion term, we obtain:

$$\frac{{{\text{Br}}}}{{{\text{Ca}}}} = \frac{{\varepsilon_{nw} K\rho_{nw0} }}{{\mu_{ref} \phi \phi_{e} L^{2} }}\frac{{\left( {n_{D} \left( {u_{nwD} } \right)_{X} } \right)_{X} }}{{S_{nw} \left( J \right)_{X} }} + \frac{{\varepsilon_{w} K\rho_{w} }}{{\mu_{ref} \phi \phi_{e} L^{2} }}\frac{{W_{wD} \left( {m_{D} \left( {u_{wD} } \right)_{X} } \right)_{X} }}{{S_{w} W_{D} \left( J \right)_{X} }}$$

(60)

All the saturation and gradient-dependent terms are normalized such that the coefficients with constant parameters are dimensionless scaling numbers that reflect the relative magnitude of the two mechanisms. For given input functions we can thus expect the Brinkman terms to become more important by increasing the magnitude of these numbers. Also, by varying the parameters within them in such a way that the numbers remain constant, we can expect the impact to be similar.

3 Results and Discussions

3.1 Input Parameters

The model reference input parameters are given in Table 1. The porous medium is assumed homogeneous in all the examples. The saturation function parameters and oil, water and rock parameters are based on experimental measurements from Bourbiaux and Kalaydjian (1990). The $J$-function parameters in (46) were fitted directly to their measured capillary pressure curve, while fluid–fluid and rock–fluid interaction coefficients and exponents were set to consistently match their co-current relative permeabilities using (50) and (51) and counter-current spontaneous imbibition data. Especially, explaining their imbibition measurements required lower mobilities during counter-current imbibition than the ones measured under steady-state relative permeability tests which was achieved by a nonzero fluid–fluid interaction coefficient $I$. Parameters obtained from those data related to our presented model were obtained in Qiao et al. (2018), where the experimental data were reproduced with a consistent set of parameters, as listed in the table. Their match of the imbibition data is shown in Fig. 2 for demonstration. Brinkman terms were then not accounted for (${\varepsilon }_{i}=0$) and the fluids were assumed incompressible. Water and oil are also here modeled as incompressible, while gas is modeled compressible according to the ideal gas law with a density of 1 kg/m³ at atmospheric pressure. Our base assumption will also be that ${\varepsilon }_{i}=0$ since the experimental data would not allow direct estimation of this parameter and to provide a reference behavior before adding the role of the Brinkman terms. Gas–water capillary pressure is obtained from oil–water capillary pressure by scaling with the IFT according to (17), both shown in Fig. 3 left. Higher initial capillary pressures lead to higher non-wetting phase initial pressures and densities. The pressure–gas density relation is shown with initial density marked for different permeabilities $K$ in Fig. 3 right.

Table 1 The model reference input parameters based on matching experimental data from Bourbiaux and Kalaydjian (1990) in the work by Qiao et al. (2018) not considering fluid compressibility or Brinkman terms

Fig. 2

Experimental validation of the model, modified from Qiao et al (2018), assuming incompressible fluids and zero viscous terms. The model captures the delay caused by viscous coupling on the effective relative permeabilities as multiphase flow changes from co-current (‘conventional’) to counter-current (‘generalized’). The experimental data are from Bourbiaux and Kalaydjian (1990)

Fig. 3

Left—Capillary pressure curves for oil–water (based on Bourbiaux and Kalaydjian 1990) and water–gas (based on scaling with IFT); Right—Pressure–density relation for gas, IFT = 72 mN/m. The maximum gas pressure and its corresponding gas density at initial condition are plotted with two different absolute permeabilities

The interaction coefficients applied in the momentum equations for the oil–water reference case are shown in Fig. 4. Note that they are proportional to viscosity divided by permeability, see Eqs. (47), (48) and (49), and will change from case to case. This ensures consistent behavior with Darcy’s law if the fluid–fluid interaction terms and viscous terms are zero. For example, setting $I=0$ in Eqs. (50) and (51) produces Corey relative permeabilities. To preserve the same amount of viscous coupling and identical counter-current relative permeabilities when defining parameters for the gas–water case the product $I{\mu }_{nw}$ was kept the same as for the oil–water case; hence, the low gas viscosity requires a higher values of $I$, see Eq. (49).

Fig. 4

Water–rock (${\widehat{{\varvec{k}}}}_{{\varvec{w}}}$), oil–rock (${\widehat{{\varvec{k}}}}_{{\varvec{o}}}$) and water–oil interaction ($\widehat{{\varvec{k}}}$) coefficients. These curves are proportional to viscosity divided by permeability and are thus case dependent

3.2 Oil–Water Simulations (Incompressible Fluids)

3.2.1 Variation of the Viscous Coupling Coefficient $I$

In Fig. 5 we have varied the viscous coupling coefficient $I$ for the oil–water reference case. As mentioned, oil and water have been assumed incompressible and we assume the Brinkman term is zero. Higher values of the coefficient $I$ increases the viscous coupling and reduces the effective relative permeabilities of both fluids, see Eq. (53). The imbibition rate is therefore reduced at greater $I$. In all cases, the recovery follows a square root of time profile at early time, as derived by Andersen et al. (2020) under these conditions. By implementing the effective relative permeabilities into a porous media simulator, in this case IORCoreSim (Lohne 2013), the same results were obtained (circles) as with our model (lines), thus validating our numerical model. The comparison was made for $I=0$ (with no viscous coupling) and $I=2300$ (the reference case). Note that the reference case demonstrates a significant reduction of imbibition rate caused by the viscous coupling effect. The saturation profiles, presented at equal times, display less water imbibed with higher viscous coupling since the imbibition rate is reduced.

Fig. 5

Simulation of water–oil displacement for different choices of the fluid–fluid interaction coefficient ${\varvec{I}}$ (in Pa⁻¹ s⁻¹) with comparison against the core scale simulator IORCoreSim (Lohne 2013) using counter-current relative permeability functions

3.2.2 Variation of Viscous Terms

First, we want to estimate the effect of the viscous forces within the wetting fluid by comparing the numerical results obtained by varying the magnitude of the wetting fluid’s effective viscosity. The relevant results are shown in Fig. 6. We can see that the water front is slower and the saturation profile is steeper when its viscosity effect is strong. The reason is that high viscosity coefficient ${\varepsilon }_{w}$ can reduce the water front velocity gradient since it requires more energy to break the front shape of displacing water. We may understand it as a gel if water has a high viscosity coefficient. A close look at the saturation profiles shows that the profiles with high viscous coefficient obtain different shapes with time: the black curves obtain lower front saturations with time (~ 0.99 at 1 h, 0.95 at 4 h, 0.9 at 8 h and 0.85 at 15 h). On the other hand, theory predicts that saturation profiles should be invariant in shape before encountering the no-flow boundary if the viscous terms and compressibility are ignored (Andersen et al. 2020). In other words, each saturation’s distance increases proportionally with the square root of time but their relative position is fixed. We also see that the curves with large viscous coefficients deviate from linear recovery profiles against square root of time, most noticeable for the black curve. Based on our simulations, the viscosity coefficient does not have strong effect on the oil recovery rate before the viscosity coefficient increases to ${\varepsilon }_{w}>{10}^{4}$ m²/s.

Fig. 6

Variation of the Brinkman coefficient of water ${{\varvec{\varepsilon}}}_{{\varvec{w}}}$. Left—wetting fluid saturation profiles at different times. Right—oil recovery profiles

In Fig. 7, we investigate the viscous effect from the non-wetting fluid (oil) by varying ${\varepsilon }_{o}$. Higher values reduce the recovery rate and delays the saturation profiles. In this case, oil moves like a gel for high oil viscosity coefficients; therefore, it is difficult to deform the front shape of displaced oil. This also will make the oil recovery process very slow and the viscous effect significant. Note that water takes significant time to accumulate at the inlet side and maintains a low saturation there until water has managed to spread across the system. The increasing saturation at the inlet is again demonstrating that the profile is not invariant and changes with time. Varying the oil coefficient has a profound impact on the recovery rate and reduces it already at ${\varepsilon }_{o}>{10}^{3}$ m²/s. For moderate values of ${\varepsilon }_{o}$ a delay (very low imbibition rates) is seen at early times and then followed by linear recovery against square root of time. That can explain the observed induction times sometimes reported experimentally before theoretical imbibition behavior is established (Akin et al. 2000; Tang and Firoozabadi 2000; Zhou et al. 2000; Føyen et al. 2019). We see especially that the recovery lines of ${\varepsilon }_{o}={10}^{3}$ and ${10}^{4}$ m²/s are parallel to the one without viscous terms. For even larger values, that does not seem to be the case. Water seems to displace oil in a less efficient way at a given value of ${\varepsilon }_{o}$ compared to if that value was set for ${\varepsilon }_{w}$: see the recovery profiles in Fig. 7 compared with those at the same times in Fig. 6 Note also that the saturation profiles are affected in opposite ways, by becoming flatter for high oil coefficients and steeper for high water coefficients.

Fig. 7

Variation of the Brinkman coefficient of oil ${{\varvec{\varepsilon}}}_{{\varvec{o}}}$. Left—wetting fluid saturation profiles at different times. Right—oil recovery profiles

The previous dominance of the oil viscous coefficient was further investigated in Fig. 8 where we vary the viscosity coefficients of both water and oil with equal values. We see a strong reduction of the recovery rate when both coefficients are increased. Lower coefficient values are required to get similar recovery response as for only varying the oil coefficient (one order of magnitude) or only the water coefficient (two orders of magnitude). The water saturation profiles at different times are similar to those in Fig. 7 where we only increase the oil viscosity effect and have low saturations near the inlet, but are steeper inside the core which may be due to the water coefficient.

Fig. 8

Variation of the Brinkman coefficients of both oil and water (with same value). Left—wetting fluid saturation profiles at different times. Right—oil recovery profiles

3.2.3 Importance of Brinkman Terms Based on Dimensionless Numbers

In this section we explore scaling based on the dimensionless numbers presented in (60) used to indicate the relative importance of the Brinkman terms compared to capillary diffusion. The numbers specifically state that if we increase the Brinkman coefficients, permeability or density or reduce viscosity, porosity or length, that will increase the impact of the Brinkman terms over the capillary diffusion term (for the same saturation functions). Note that this is based on the assumption that the Brinkman coefficients are constant.

In Fig. 9 we show the reference case with zero Brinkman terms (blue curve) and with Brinkman coefficients ${10}^{4}$ m²/s for oil and water (red curve) which has caused some delay in recovery. Compared to this case we make three new cases where either the Brinkman coefficients are increased by a factor 4 (black ‘ + ’ markers), the length is halved (green curve) or permeability is increased by a factor 4 (purple ‘*’ points). This results in the same increased ratio of $\frac{{\varepsilon }_{i}K}{{L}^{2}}$ of 4 compared to the red curve case and all three cases fall on the same line as predicted by the scaling where keeping the dimensionless numbers same should give same behavior.

Fig. 9

Demonstration of scaling the relative importance of Brinkman terms over capillary diffusion. The blue curve has zero Brinkman terms, while the red curve has significant delay by the Brinkman terms. The remaining cases have increased the ratio $\frac{{\varepsilon }_{i}K}{{L}^{2}}$ by a factor 4 by either changing length, permeability or Brinkman coefficients which has given the same result as predicted by the scaling

3.2.4 Match of Experiments with Induction Time

Inspired by the previous observations, we here apply the developed model to explain some experiments from the literature that display nonstandard behavior. Zhou et al. (2000) conducted forced and spontaneous imbibition experiments on cores with different wetting states, as determined by how long they had been aged with a crude oil. Behbahani and Blunt (2005) matched the experiments based on core scale simulation and consistent saturation functions from pore scale. The experimental data did, however, display induction time in some cases, where the imbibition rate stayed low for an early period before developing into more significant magnitudes. Our goal was to improve the match of two spontaneous imbibition experiments by using Brinkman terms to account for the induction time, i.e., the period with low rates.

To model the experiments, the known fluid and rock parameters were input to the model such as oil and water viscosities, IFT, porosity, permeability, core dimensions and initial saturations, see Table 2. Gravity and fluid compressibility effects were ignored. As our model was 1D and the experiments were conducted with all faces open on cylindrical cores, we used the Ma et al. (1997) characteristic length, which for a cylinder with height $H$ and diameter $D$ is:

$$L_{c} = \frac{DH}{{2\left( {D^{2} + 2H^{2} } \right)^{0.5} }}$$

(61)

Table 2 Experimental rock and fluid parameters from Zhou et al. (2000) and momentum equation and capillary pressure parameters used to match two of their experiments

This length is the effective distance to the no-flow boundary. As input parameters for the $J$-function and the momentum equations, we respectively tuned the correlation (46) to the capillary pressure data obtained by Behbahani and Blunt (2005) and also tuned the momentum equation parameters in the counter-current relative permeabilities (53) to their relative permeabilities. As they had used the same relative permeabilities to model co- and counter-current imbibition and not captured the induction time (as is not possible with standard assumptions), these data were just treated as initial guesses. As we focused only on the spontaneous imbibition experiments, there was no information to determine the fluid–fluid interaction parameter $I$ and it was set to 0.

The experimental recovery data were first plotted against the square root of time and a linear trend was identified, see dashed/dotted lines in Fig. 10. As this linear trend did not start at zero recovery at zero time, but after a delay (induction time), the slope of the line was first matched to produce a recovery trend that was parallel to the experimental line (i.e., the shifted full lines in Fig. 10). In other words, during the main imbibition process the imbibition rate was matched by the model. This was done by tuning the magnitude of the $J$-function, but not its shape (especially the point where capillary pressure crosses zero determines end recovery during spontaneous imbibition). This resulted in a recovery curve with identical shape as the experimental one, but shifted on the time axis (full lines in Fig. 11). To match the time shift, the oil Brinkman term coefficient was increased from zero. This shifted the data sufficiently and continuously to provide a successful match of the experiments (see the dashed-dotted lines in Fig. 11). The final parameters used to match the data are listed in Table 2. It could be possible to improve the match of the late time recovery of the 48-h test by modifying the relative permeabilities; however, our main goal was to demonstrate the usefulness of the Brinkman terms to model induction time data.

Fig. 10

Experimental counter-current imbibition data from Zhou et al. (2000) where two cores were aged in crude oil with different aging time (indicated). Linear trends in the data were identified (dashed/dotted lines) between recovery and the square root of time, which did not start at zero time. Shifting the lines to zero time (full lines) indicate the amount of induction time

Fig. 11

Simulation of the experimental counter-current imbibition data from Zhou et al. (2000) where two cores were aged in crude oil with different aging time (indicated). First the slopes of the linear data were matched assuming zero Brinkman terms (full lines), then the original data were matched by including nonzero Brinkman terms to capture the induction time

Spontaneous imbibition has been studied during the induction time with in-situ images (Føyen et al. 2019). They indicate a period where local fronts form and develop before merging to more symmetrical profiles. Induction time has been associated with more oil-wet state or low initial water saturation (Tang and Firoozabadi 2000; Zhou et al. 2000), non-uniform wetting or interaction with epoxy used to seal core surfaces (Føyen et al. 2019); conditions that complicate the development of a continuous water film. As there are many reservoirs with mixed- or oil-wet state and low water saturation, this mechanism can increase the imbibition time. By tracing the straight lines in Fig. 10 to the end of the linear data period we see that in the 72-h test the induction time caused the duration of the imbibition to increase by a factor of $\sim {\left(\frac{110}{65}\right)}^{2}=2.9$, while for the 48-h test it increased by $\sim {\left(\frac{60}{50}\right)}^{2}=1.44$, which both are significant. In many of Tang and Firoozabadi (2000)’s tests the induction time lasted 25–30% of the full imbibition time, i.e., the time increased by a factor ~ 1.3. Although a standard model could capture a longer time scale, it would not capture the variation in imbibition rate from low to high and then low again. How this effect upscales and depends on different parameters is not much explored.

We note that so-called non-equilibrium models (Barenblatt et al. 2002; Silin and Patzek 2004) exist as modifications of classical models which also could offer better adaption to experimental observations. In those models, the relative permeability and capillary pressure functions are expressed using effective saturations as defined by the saturation plus a relaxation time coefficient multiplied by the saturation time derivative. By selecting a proper saturation dependent relaxation time coefficient Silin and Patzek (2004) obtained a recovery solution with a power function relation at early time and square root of time at later time. They were able to match experimental data, but did not demonstrate how significant the improvement was compared to classical approaches. The motivation of these models has been experimental findings where imbibition does not display self-similar behavior, i.e., where recovery does not follow the square root of time and saturation profiles do not overlap when plotted against a similarity variable. We point out that several of these experiments were conducted using air as non-wetting fluid in low permeability media, e.g., by Le Guen and Kovscek (2006). As we will see in the following section, gas compressibility can explain deviation from self-similarity.

3.3 Gas–Water Simulations (Compressible Non-Wetting Fluid)

3.3.1 Impact of Gas Compressibility (by Variation of K) Without Viscous Terms

The effect of gas compressibility is investigated in this section. The gas is modeled as ideal, which is realistic at lab conditions, by setting ${C}_{g}={10}^{5}$ and ${\stackrel{\sim }{\rho }}_{g0}=0$. The gas then has an initial density proportional to the initial non-wetting pressure, which increases with capillary pressure, see (45). By considering a system with lower permeability, high water-wetness, or high IFT the initial gas density will be high. The boundary pressure (1 bar), and gas density at that pressure (1 kg/m³), is the same in all cases. As the system approaches zero capillary pressure where both phases reach ambient pressure, the gas density in the system will approach the one at the boundary (1 kg/m³). We expect gas compressibility to be more significant when the initial density is high compared to the end density, i.e., there is a significant gas expansion from start to finish. We test cases at different absolute permeabilities to obtain different initial densities. As seen in Fig. 3 right, where density is plotted against pressure, the reference permeability of 118 mD gives an initial density almost identical to 1 kg/m³ at which we expect negligible compressibility effects. Reducing the permeability by factors of 10 to 11.8 and 1.18 mD gives an initial density of ~ 1.6 kg/m³ and ~ 5 kg/ m³, respectively. In the latter case the gas thus expands fivefold from its initial state to its end state. For contrast, the examples are also generated where the gas phase is made incompressible case with a constant density of 1 kg/ m³; characterized by setting ${\stackrel{\sim }{\rho }}_{g0}=1$ kg/ m³ and ${C}_{g}={10}^{10}$.

In Fig. 12 the results of the three cases (compressible and incompressible) are shown as water saturation profiles at equal times and recovery against square root of time. We note that reducing permeability has two effects: (1) it reduces the imbibition rate and delays imbibition profiles, and (2) gives greater separation between the ‘compressible’ and ‘incompressible’ cases with same permeability. For the reference permeability there is no real difference, as expected by the negligible density change. For the cases with 10 and 100 times lower permeability, the separation becomes more significant.

Fig. 12

Effect of gas compressibility on imbibition front and recovery factor shown in regular time; interfacial tension is 72 mN/m

To eliminate the effect of permeability on the time scale, the saturation profiles are shown at same scaled times and recovery is plotted against the square root of scaled time in Fig. 13. The scaling causes all incompressible cases and the reference compressible case (118 mD) to overlap. This is expected since when only the rock–fluid and fluid–fluid effects are significant the imbibition rate is proportional to the square root of permeability. The only difference between the curves is then related to the gas compressibility. Inspecting the figures, one can see that the compressibility accelerates the imbibition process. The ‘compressible’ recovery profiles are always higher than the ‘incompressible’ recovery profiles, and the ‘compressible’ saturation profiles have advanced farther than the ‘incompressible’ ones. The effect is most significant at early time when gas has the highest density, but small saturation changes can give large capillary pressure changes and cause gas expansion which contributes to drive gas out. At $t=10\tau$ ($\sqrt{t/\tau }=3.2$) water saturations have increased over the entire core, reducing capillary pressure and non-wetting pressure. Most of the gas expansion has then occurred. There is a visible decline in imbibition rate even earlier, already at $\sqrt{t/\tau }=2$ for the 1.18 mD case with fast water front movement. This case has a recovery factor which is 0.20 higher for the compressible cases compared to the incompressible case for a significant period of time. The imbibition processes, however, appear to end at similar times.

Fig. 13

Effect of gas compressibility on imbibition front and recovery factor shown in scaled time; the interfacial tension is 72 mN/m

3.3.2 Matching Experimental Data

Counter-current 1D spontaneous imbibition experiments were conducted by Roychaudhuri et al. (2013) where deionized water imbibed into cubic shale samples with only one side open for flow, and gas was displaced. Shales have very low permeability and are developed by hydraulic fracturing with mainly water-based fluids. Shales can consequently develop large capillary pressures, which we predict to cause significant gas compression. Imbibition in shale has demonstrated unexpected imbibition behaviors and particularly Roychaudhuri et al. (2013) observed that their imbibition data obtained two linear trends when recovery was plotted against the square root of time. Standard modeling should only display one linear trend, but as seen from our sensitivity analysis in Fig. 13, significant gas compressibility causes two apparent linear slopes. We thus expected the inclusion of compressibility to better match their experiments and did so for their test #17.

Without measured saturation functions for those samples, we tuned the reference saturation functions and compared recovery data scaled from 0 to 1. The time scale of the test depends on the magnitude of the capillary diffusion coefficient, which combines fluid mobilities, permeability, $J$-function, porosity and interfacial tension. The impact of compressibility depends on the initial capillary pressure, which depends on permeability, porosity, interfacial tension and $J$-function, but not fluid mobilities. Representative parameters were assumed for air and deionized water, while arbitrary end saturations were assumed. Porosity, permeability and length were from the experimental measurements. The magnitude of the $J$-function, as well as the magnitude of the relative permeabilities were tuned to match the experiments. The applied parameters are listed in Table 3. Notably, the relative permeabilities required lower magnitudes than the base functions (as seen by the high ${I}_{w},{I}_{nw},I$), while the $J$-function had a similar magnitude as the base function. The best match of the data indicated that the initial gas density was 4.6 kg/m³ which suggests significant expansion toward 1 kg/m³ during the recovery process.

Table 3 Input parameters to match the experimental data from Roychaudhuri et al. (2013)

The experimental data are shown in Fig. 14 as recovery against square root of time. Two straight lines are plotted to indicate the two distinct linear sections. It is seen that the model captures very well the two linear slopes drawn through the experimental data when gas compressibility is accounted for. Treating the gas as an incompressible fluid results in a straight line. Our model hence explains the untraditional observation of two linear slopes of recovery vs square root of time during counter-current imbibition in shale. It is well documented that shales can display nonstandard flow and storage mechanisms such as slip, adsorption, Knudsen diffusion, confinement of nanopores, etc. (Sun et al. 2015; Pitakbunkate et al. 2016; Zhang et al. 2019; Klewiah et al. 2020). However, their importance and to what extent effective Darcy properties can capture them is debated. We here find that the main nonstandard behavior was captured by accounting for the compressibility of the gas, which our model predicts has more impact in tight media with high capillarity.

Fig. 14

Experimental data (triangles) from test #17 in Roychaudhuri et al. (2013) of water–gas counter-current imbibition in shale. Recovery vs square root of time displays two linear slopes which could be matched with our model by accounting for gas compressibility. Treating the gas as incompressible results in a classical straight line

3.3.3 Combined Impact of Compressibility and Viscous Terms

Finally, we are interested in estimating the combined effect of Brinkman viscous terms and gas compressibility. The Brinkman coefficient for both fluids was varied with equal values for both fluids and the absolute permeability was set at 1.18 mD to get a pronounced compressibility effect. Corresponding incompressible cases were also generated. The results are shown with scaled time in Fig. 15. For compressible or incompressible cases, an increase in the Brinkman coefficients $\varepsilon$ reduces recovery rate. Compressible cases have higher recovery than their corresponding incompressible case. It also seems that compressibility plays a dominating role on recovery over the Brinkman coefficient at early times. For $\varepsilon \le {10}^{4}$ m²/s the recovery profiles of the compressible cases are close to identical up to a recovery of $RF=0.5$, but differ at higher recoveries (later times), with high Brinkman coefficients yielding lower recovery.

Fig. 15

Effect of Brinkman coefficient (same value is used for both fluids) during gas–water imbibition. The effect of compressible and incompressible gas is also shown. The absolute permeability and interfacial tension are 1.18 mD and 72 mN/m respectively

4 Conclusions

In this paper we have studied counter-current spontaneous imbibition driven by capillary forces. The model was formulated based on mass balance and momentum balance equations which account for fluid–rock interactions, fluid–fluid interaction, fluid compressibility and Brinkman terms. The model was parameterized based on a previous match of literature experimental data and validated against a commercial software for conditions where that was possible. Sensitivity analyses were carried out for a water–oil and a water–gas system and used to assess the potential role and importance of the different mechanisms.

• (1) Based on water–oil flow simulations, we observe that the impact of increasing the Brinkman term of the displaced fluid (oil) has more impact on the fluid flow than the Brinkman term of the displacing fluid (water). For both terms, a larger value can delay the recovery process.

• (2) With a significant Brinkman term for an incompressible non-wetting phase the recovery curve was delayed at early times, but displayed linear behavior against the square root of time at later times. This could explain and match oil–water spontaneous imbibition experiments with induction time reported in the literature.

• (3) Using scaling a dimensionless number was derived to indicate when Brinkman terms were important relative to capillary diffusion. For other parameters constant, same impact was seen when Brinkman coefficient times permeability divided by length squared was fixed. Stronger impact was seen by increasing this ratio.

• (4) We found that recovery of displaced fluid was affected by its compressibility behavior, here defined as an ideal gas. Incompressible displaced fluid was displaced slower than if it was compressible, but the difference was only significant when the initial non-wetting pressure was so much greater than the external ambient pressure that the gas could expand significantly.

• (5) Conditions where gas compressibility is significant could be typical and our results demonstrated that strongly water-wet media with permeability less than 10 mD will have higher imbibition rates than if the fluids were incompressible, which affects simulation of the experiments. The recovery process seemed to end at comparable times.

• (6) Gas compressibility could explain why water–gas imbibition in shale shows two linear slopes of recovery in a square root of time plot. It is mainly because compressibility is more significant than higher permeabilities. Literature experiments could be explained with this mechanism included. The two-slope observation should be expected for any porous media with high capillarity during water–gas counter-current imbibition.

• (7) Combining the effects from both fluid compressibility and Brinkman viscous terms, the compressibility effect could overcome the recovery delay resulting from the Brinkman viscous terms, leading to a fast gas recovery process.

Appendices

Appendix A

1.1 Discretization

The length scale is normalized by introducing dimensionless coordinate:

$$x_{D} = \frac{x}{L}$$

(62)

The Eqs. (23) to (26) become:

$$\left( m \right)_{t} + \left( {mu_{w} } \right)_{x} /L = 0$$

(63)

$$s_{w} (P_{w} )_{x} /L = - \hat{k}_{w} u_{w} + \hat{k}\left( {u_{nw} - u_{w} } \right) + \varepsilon_{w} \left( {m\left( {u_{w} } \right)_{x} } \right)_{x} /L^{2}$$

(64)

$$s_{o} (P_{w} + P_{c} )_{x} /L = - \hat{k}_{nw} u_{nw} + \hat{k}\left( {u_{w} - u_{nw} } \right) + \varepsilon_{nw} \left( {n\left( {u_{nw} } \right)_{x} } \right)_{x} /L^{2}$$

(65)

$$(P_{w} )_{t} + \tilde{\eta }\rho_{w} \left( {nu_{nw} } \right)_{x} /L + \tilde{\eta }\tilde{\rho }_{nw} \left( {mu_{w} } \right)_{x} /L = 0$$

(66)

where we skipped notation “D” for ${x}_{D}$ for reading simplicity. The spatial domain [0,1] is divided uniformly into M computational cells. A staggered grid (Fig. 

Fig. 16

Staggered grid used for discretization

16) is used where the values for masses ${m}_{i},{n}_{i}$ and pressures ${P}_{w}, {P}_{nw},{P}_{c}$ are located in the center of the cell and velocities ${u}_{w},{u}_{nw}$ on the interfaces.

On the left boundary we set water pressure and capillary pressure as:

$$P_{w} |_{{{\text{j}} = 0}} = 10^{5} {\text{Pa}},{ }P_{c} |_{{{\text{j}} = 0}} = 0$$

(67)

The boundary water saturation is calculated by inverting the capillary pressure function:

$$s_{w} |_{{{\text{j}} = 0}} = P_{c}^{ - 1} \left( {P_{c} \left. \right|_{{{\text{j}} = 0}} } \right)$$

(68)

Since the right end is closed we set the fluids velocities there to zero:

$$u_{w} |_{N + 1/2} = u_{o} |_{N + 1/2} { } = { }0$$

(69)

A time domain $\left[ {0,T} \right]$ is discretized using a logarithmic scale in order to provide sufficient accuracy at the initial stage of the imbibition process which is characterized by large velocity gradients. First the mass conservation equation is discretized:

$$\frac{{m_{j}^{k + 1} - m_{j}^{k} }}{\Delta t} = - \frac{1}{\Delta xL}\left( {\left[ {mu_{w} } \right]_{{j + \frac{1}{2}}}^{k + 1} - \left[ {mu_{w} } \right]_{{j - \frac{1}{2}}}^{k + 1} } \right)$$

(70)

where the upwind scheme is used to discretize the flux:

$$[mu_{w} ]_{{j + \frac{1}{2}}}^{k + 1} = \left\{ {\begin{array}{*{20}c} {m_{j}^{k} u_{w,j + 1/2}^{k} , if u_{w,j + 1/2}^{k} \ge 0;} \\ {m_{j + 1}^{k} u_{w,j + 1/2}^{k} , if u_{w,j + 1/2}^{k} \ge 0} \\ \end{array} } \right.$$

(71)

In order to get coefficients in the next step, we need to update the water saturation:

$$s_{w,j}^{k + 1/2} = \frac{{m_{j}^{k + 1} }}{{\rho_{w} \left( {P_{w,j}^{k} } \right)}}$$

(72)

non-wetting fluid pressure:

$$P_{o,j}^{k + 1/2} = P_{w,j}^{k} + P_{c} \left( {s_{w,j}^{{k + \frac{1}{2}}} } \right)$$

(73)

and the non-wetting fluid mass:

$$n_{j}^{k + 1/2} = \left( {1 - s_{w,j}^{{k + \frac{1}{2}}} } \right)\rho_{o} \left( {P_{o,j}^{{k + \frac{1}{2}}} } \right)$$

(74)

The next step is simultaneous computation of wetting fluid pressure and wetting and non-wetting fluid velocities with pressure evolution and momentum conservation equations discretized as follows:

$$\frac{{m_{j}^{k + 1} - m_{j}^{k} }}{\Delta t} = - \frac{1}{\Delta xL}\left( {\left[ {mu_{w} } \right]_{{j + \frac{1}{2}}}^{k + 1} - \left[ {mu_{w} } \right]_{{j - \frac{1}{2}}}^{k + 1} } \right)$$

(75)

$$\begin{aligned} s_{{o,j + \frac{1}{2}}}^{{k + \frac{1}{2}}} \frac{1}{{\Delta x}}\left( {P_{{w,j + 1}}^{{k + 1}} - P_{{w,j}}^{{k + 1}} } \right) = & - s_{{o,j + \frac{1}{2}}}^{{k + \frac{1}{2}}} \frac{1}{{\Delta x}}\left( {P_{{c,j + 1}}^{{k + \frac{1}{2}}} - P_{{c,j}}^{{k + \frac{1}{2}}} } \right) - L \cdot \hat{k}_{{o,j + \frac{1}{2}}}^{{k + \frac{1}{2}}} u_{{o,j + \frac{1}{2}}}^{{k + 1}} \\ & \quad \quad + L \cdot \hat{k}_{{j + \frac{1}{2}}}^{{k + \frac{1}{2}}} \left( {u_{{w,j + \frac{1}{2}}}^{{k + 1}} - u_{{o,j + \frac{1}{2}}}^{{k + 1}} } \right) \\ & \quad \quad + \varepsilon _{o} \frac{1}{{\Delta x^{2} L}}\left( {n_{{j + 1}}^{{k + \frac{1}{2}}} \left[ {u_{{o,j + \frac{3}{2}}}^{{k + 1}} - u_{{o,j + \frac{1}{2}}}^{{k + 1}} } \right] - n_{j}^{{k + \frac{1}{2}}} \left[ {u_{{o,j + \frac{1}{2}}}^{{k + 1}} - u_{{o,j - \frac{1}{2}}}^{{k + 1}} } \right]} \right) \\ \end{aligned}$$

(76)

$$\begin{aligned} s_{{w,j + \frac{1}{2}}}^{{k + \frac{1}{2}}} \frac{1}{{\Delta x}}\left( {P_{{w,j + 1}}^{{k + 1}} - P_{{w,j}}^{{k + 1}} } \right) = & - L \cdot \hat{k}_{{w,j + \frac{1}{2}}}^{{k + \frac{1}{2}}} u_{{w,j + \frac{1}{2}}}^{{k + 1}} - L \cdot \hat{k}_{{j + \frac{1}{2}}}^{{k + \frac{1}{2}}} \left( {u_{{w,j + \frac{1}{2}}}^{{k + 1}} - u_{{o,j + \frac{1}{2}}}^{{k + 1}} } \right) \\ & \quad + \varepsilon _{w} \frac{1}{{\Delta x^{2} L}}\left( {m_{{j + 1}}^{{k + \frac{1}{2}}} \left[ {u_{{w,j + \frac{3}{2}}}^{{k + 1}} - u_{{w,j + \frac{1}{2}}}^{{k + 1}} } \right] - m_{j}^{{k + \frac{1}{2}}} \left[ {u_{{w,j + \frac{1}{2}}}^{{k + 1}} - u_{{w,j - \frac{1}{2}}}^{{k + 1}} } \right]} \right) \\ \end{aligned}$$

(77)

$$\begin{gathered} L\frac{{P_{w,j}^{k + 1} - P_{w,j}^{k} }}{\Delta t} + [\tilde{\eta }\rho_{w} ]_{j}^{k + 1} \frac{1}{\Delta x}\left( {\left[ {n^{{k + \frac{1}{2}}} u_{o}^{k + 1} } \right]_{{j + \frac{1}{2}}} - \left[ {n^{{k + \frac{1}{2}}} u_{o}^{k + 1} } \right]_{{j - \frac{1}{2}}} } \right) \hfill \\ + \left[ {\tilde{\eta }\tilde{\rho }_{o} } \right]_{j}^{k + 1} \frac{1}{\Delta x}\left( {\left[ {m^{k + 1} u_{w}^{k + 1} } \right]_{{j + \frac{1}{2}}} - \left[ {m^{k + 1} u_{w}^{k + 1} } \right]_{{j - \frac{1}{2}}} } \right) = 0 \hfill \\ \end{gathered}$$

(78)

where

$$s_{{w,j + \frac{1}{2}}}^{k + 1/2} = \left\{ {\begin{array}{*{20}c} {s_{w,j}^{k + 1/2} , if u_{w,j + 1/2}^{k} > 0;} \\ {\left( {s_{w,j}^{{k + \frac{1}{2}}} + s_{w,j + 1}^{{k + \frac{1}{2}}} } \right)/2, if u_{{w,j + \frac{1}{2}}}^{k} = 0;} \\ {s_{w,j + 1}^{k + 1/2} , if u_{w,j + 1/2}^{k} < 0} \\ \end{array} } \right.$$

(79)

Values for ${s}_{w,j+\frac{1}{2}}^{k+1/2}$, ${\widehat{k}}_{w,j+\frac{1}{2}}^{k+1/2}$ and, ${\widehat{k}}_{o,j+\frac{1}{2}}^{k+1/2}$ follow the same logic.

$$\hat{k}_{{j + \frac{1}{2}}}^{k + 1/2} = \left\{ {\begin{array}{*{20}l} {\hat{k}_{j}^{k + 1/2} ,} \hfill & {{\text{if}}\,u_{w,j + 1/2}^{k} > 0 \left| {\& u_{o,j + 1/2}^{k} > 0;} \right.} \hfill \\ {\left( {\hat{k}_{j}^{{k + \frac{1}{2}}} + \hat{k}_{j + 1}^{{k + \frac{1}{2}}} } \right)/2,} \hfill & {{\text{ if}}\,{ }u_{w,j + 1/2}^{k} \cdot u_{o,j + 1/2}^{k} \le 0} \hfill \\ {\hat{k}_{j + 1}^{k + 1/2} ,} \hfill & {{\text{if}}\,{ }u_{w,j + 1/2}^{k} < 0\left| {\& u_{o,j + 1/2}^{k} < 0} \right.{ }} \hfill \\ \end{array} } \right.$$

(80)

There are $3\left(N-1\right)$ equations for $3\left(N-1\right)$ unknowns.

On the left boundary we set the non-wetting fluid pressure: ${P}_{w}{|}_{\mathrm{x}=0}={10}^{5}\mathrm{ Pa}$

We also set the mass flux derivatives for both fluids on the boundary block to be zero:

$$\frac{\partial }{\partial x}\left( {mu_{w} } \right)|_{0} { } = { }0,\frac{\partial }{\partial x}\left( {mu_{o} } \right)|_{0} { } = { }0$$

(81)

$$\begin{aligned} s_{{o,\frac{1}{2}}}^{{k + \frac{1}{2}}} \frac{1}{\Delta x}\left( {P_{w,1}^{k + 1} - 0} \right) = & \quad - s_{{o,\frac{1}{2}}}^{{k + \frac{1}{2}}} \frac{1}{\Delta x}\left( {P_{c,1}^{{k + \frac{1}{2}}} - 0} \right) \\ & \quad - L \cdot \hat{k}_{{o,\frac{1}{2}}}^{{k + \frac{1}{2}}} u_{{o,j + \frac{1}{2}}}^{k + 1} + \hat{k}_{\frac{1}{2}}^{{k + \frac{1}{2}}} \left( {u_{{w,\frac{1}{2}}}^{k + 1} - u_{{o,\frac{1}{2}}}^{k + 1} } \right) \\ & \quad \quad + \varepsilon_{o} \frac{1}{{\Delta x^{2} L}}\left( {n_{1}^{{k + \frac{1}{2}}} \left[ {u_{{o,\frac{3}{2}}}^{k + 1} - u_{{o,\frac{1}{2}}}^{k + 1} } \right]} \right) \\ \end{aligned}$$

(82)

$$\begin{aligned} s_{{w,\frac{1}{2}}}^{{k + \frac{1}{2}}} \frac{1}{\Delta x}\left( {P_{w,1}^{k + 1} - 0} \right) = & \quad - L \cdot \hat{k}_{{w,\frac{1}{2}}}^{{k + \frac{1}{2}}} u_{{w,\frac{1}{2}}}^{k + 1} - L \cdot \hat{k}_{\frac{1}{2}}^{{k + \frac{1}{2}}} \left( {u_{{w,\frac{1}{2}}}^{k + 1} - u_{{o,\frac{1}{2}}}^{k + 1} } \right) \\ & \quad \quad + \varepsilon_{w} \frac{1}{{\Delta x^{2} L}}\left( {m_{1}^{{k + \frac{1}{2}}} \left[ {u_{{w,\frac{3}{2}}}^{k + 1} - u_{{w,\frac{1}{2}}}^{k + 1} } \right]} \right) \\ \end{aligned}$$

(83)

On the right boundary node (no-flux boundary), ${u}_{w}{|}_{N+1/2},{u}_{o}{|}_{N+1/2} = 0$ we have:

$$\frac{{m_{N}^{k + 1} - m_{N}^{k} }}{\Delta t} = - \frac{1}{\Delta x}(0 - \left[ {mu_{w} ]_{N - 1/2}^{k + 1} } \right)$$

(84)

$$\begin{aligned} s_{{o,N - \frac{1}{2}}}^{{k + \frac{1}{2}}} \frac{1}{\Delta x}\left( {P_{w,N}^{k + 1} - P_{w,N - 1}^{k + 1} } \right) & = - s_{{o,N - \frac{1}{2}}}^{{k + \frac{1}{2}}} \frac{1}{\Delta x}\left( {P_{c,N}^{{k + \frac{1}{2}}} - P_{c,N - 1}^{{k + \frac{1}{2}}} } \right) - \hat{k}_{{o,N - \frac{1}{2}}}^{{k + \frac{1}{2}}} u_{{o,N - \frac{1}{2}}}^{k + 1} + \hat{k}_{{N - \frac{1}{2}}}^{{k + \frac{1}{2}}} \left( {u_{{w,N - \frac{1}{2}}}^{k + 1} - u_{{o,N - \frac{1}{2}}}^{k + 1} } \right) \\ & \quad + \varepsilon_{o} \frac{1}{{\Delta x^{2} }}\left( {n_{N}^{{k + \frac{1}{2}}} \left[ {0 - u_{{o,N - \frac{1}{2}}}^{k + 1} } \right] - n_{N - 1}^{{k + \frac{1}{2}}} \left[ {u_{{o,N - \frac{1}{2}}}^{k + 1} - u_{{o,N - 1\frac{1}{2}}}^{k + 1} } \right]} \right) \\ \end{aligned}$$

(85)

$$\begin{aligned} s_{{w,N - \frac{1}{2}}}^{{k + \frac{1}{2}}} \frac{1}{\Delta x}\left( {P_{w,N}^{k + 1} - P_{w,N - 1}^{k + 1} } \right) = & - \hat{k}_{{w,N - \frac{1}{2}}}^{{k + \frac{1}{2}}} u_{{w,N - \frac{1}{2}}}^{k + 1} - \hat{k}_{{N - \frac{1}{2}}}^{{k + \frac{1}{2}}} \left( {u_{{w,N - \frac{1}{2}}}^{k + 1} - u_{{o,N - \frac{1}{2}}}^{k + 1} } \right) \\ & \quad + \varepsilon_{w} \frac{1}{{\Delta x^{2} }}\left( {m_{N}^{{k + \frac{1}{2}}} \left[ {0 - u_{{w,N - \frac{1}{2}}}^{k + 1} } \right] - m_{N - 1}^{{k + \frac{1}{2}}} \left[ {u_{{w,N - \frac{1}{2}}}^{k + 1} - u_{{w,N - 1\frac{1}{2}}}^{k + 1} } \right]} \right) \\ \end{aligned}$$

(86)

$$\frac{{P_{w,N}^{k + 1} - P_{w,N}^{k} }}{\Delta t} + [\tilde{\eta }\rho_{w} ]_{N}^{k + 1} \frac{1}{\Delta x}\left( {0 - \left[ {n^{{k + \frac{1}{2}}} u_{o}^{k + 1} } \right]_{{N - \frac{1}{2}}} } \right) + \left[ {\tilde{\eta }\tilde{\rho }_{o} } \right]_{N}^{k + 1} \frac{1}{\Delta x}\left( {0 - \left[ {m^{{k + \frac{1}{2}}} u_{w}^{k + 1} } \right]_{{N - \frac{1}{2}}} } \right) = 0$$

(87)

Finally, we update the wetting fluid saturation:

$$s_{w,j}^{k + 1} = \frac{{m_{j}^{k + 1} }}{{\rho_{w} \left( {P_{w,j}^{k + 1} } \right)}}$$

(88)

Appendix B

2.1 Grid and Temporal Sensitivity Analysis

A grid sensitivity analysis was performed to check that the numerical solutions were not sensitive to the number of cells. The results are shown in Fig. 

Fig. 17

Oil recovery profiles produced by the script in the present study with the number of nodes in the spatial grids varied. The time step was also varied in the case with $I = 2300$ and 200 nodes

17 for the reference case with two values of $I$ and indicate little change once 200 grid cells have been used. Temporal convergence analysis was also performed by time step variation for the case with $I = 2300$ and 200 grid cells. 0.5 s was approximately the stability limit for that case and two cases with half and a quarter of that time step were simulated. It is indicated that time step size did not affect the solution as long as it was small enough for the simulation to be stable.

Appendix C

Global mass error analysis

Consider the water mass balance Eq. (23) repeated below:

$$\left( m \right)_{t} + \left( {mu_{w} } \right)_{x} = 0,$$

(89)

We integrate it over the core from $x = 0$ to $L$ and divide by $L$.

$$\left( {\overline{m}} \right)_{t} - \frac{1}{{\text{L}}}\left( {mu_{w} } \right)|_{x = 0} = 0, \overline{m} = \frac{1}{L}\mathop \smallint \limits_{x = 0}^{L} mdx.$$

(90)

In the above, $\overline{m}$ is the spatially averaged value of $m$ across the system at a given time. Further, integrate over time to obtain:

$$\left[ {\overline{m}\left( t \right) - \overline{m}_{0} } \right] - \frac{1}{{\text{L}}}\mathop \smallint \limits_{t = 0}^{t} \left( {mu_{w} } \right)|_{x = 0} dt = 0.$$

(91)

This equation tells us that the change in the average conserved variable $\overline{m}$ is due to transport across the open boundary. If the numerical scheme introduces any error, the difference between the two terms might be nonzero. We can therefore define a relative global mass error as follows:

$$E = \frac{{\left[ {\overline{m}\left( t \right) - \overline{m}_{0} } \right] - \frac{1}{{\text{L}}}\mathop \smallint \nolimits_{t = 0}^{t} \left( {mu_{w} } \right)|_{x = 0} dt}}{{\left[ {\overline{m}\left( t \right) - \overline{m}_{0} } \right]}},$$

(92)

To evaluate the magnitude of this error, we ran the reference case (with 200 nodes) and plotted the relative error $E$ in Fig. 

Fig. 18

Relative global water mass error of the reference case as a function of time

18. We clearly see that the water mass is conserved with a rather low value of relative error $E$ defined in (92), more than 12 orders of magnitude lower than one percent of the total mass change.