Finite-Volume Simulation of Capillary-Dominated Flow in Matrix-Fracture Systems using Interface Conditions

TL;DR

This paper develops a finite-volume simulation approach with discrete interface conditions to accurately model capillary-driven flow in fractured porous media, improving flux estimation without fine grid refinement.

Contribution

It introduces a novel application of interface conditions in finite-volume schemes to enhance flux accuracy in multiphase flow simulations in fractured media.

Findings

Interface conditions improve flux accuracy at coarse grid resolutions.

The method accurately captures capillary-driven imbibition dynamics.

The approach effectively models flow with viscous, buoyancy, and capillary forces.

Abstract

In numerical simulations of multiphase flow and transport in fractured porous media, the estimation of the hydrocarbon recovery requires accurately predicting the capillary-driven imbibition rate of the wetting phase initially present in the fracture into the low-permeability matrix. In the fully implicit finite-volume scheme, this entails a robust methodology that captures the capillary flux at the interface between the matrix and the fracture even when very coarse control volumes are used to discretize the matrix. Here, we investigate the application of discrete interface conditions at the matrix-fracture interface to improve the accuracy of the flux computation without relying on extreme grid refinement. In particular, we study the interaction of the upwinding scheme with the discrete interface conditions. Considering first capillary-dominated spontaneous imbibition and then forced imbibition with viscous, buoyancy, and capillary forces, we illustrate the importance of the interface conditions to accurately capture the matrix-fracture flux and correctly represent the flow dynamics in the problem.

Figures

Equations

$$\varphi \frac{\partial S_{\ell }}{\partial t}+\nabla \cdot {\mathbf{u}}_{\ell }=q_{\ell }\forall \ell \in \{w,n\},$$

$${\mathbf{u}}_{\ell }:=-k{\lambda }_{\ell }(\nabla p_{\ell }-{\rho }_{\ell }g\nabla z)\forall \ell \in \{w,n\}.$$

$$S_n+S_w=1.$$

$$P_c^{(\alpha )}(S_w)=p_n-p_w.$$

$${\mathbf{u}}_T:={\mathbf{u}}_w+{\mathbf{u}}_n=-k{\lambda }_T\nabla p+k({\lambda }_w{\rho }_w+{\lambda }_n{\rho }_n)g\nabla z+k{\lambda }_w\nabla P_c^{(\alpha )},$$

$$\nabla \cdot {\mathbf{u}}_T=q_w+q_n.$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi\frac{\partial S}{\partial t}+\nabla\cdot\Bigg{(}\underbrace{\frac{\lambda_{w}}{\lambda_{T}}\boldsymbol{u}_{T}}_{\begin{subarray}{c}\text{viscous}\\ \text{term}\end{subarray}}+\underbrace{k\frac{\lambda_{w}\lambda_{n}}{\lambda_{T}}(\rho_{w}-\rho_{n})g\nabla z}_{\begin{subarray}{c}\text{buoyancy}\\ \text{term}\end{subarray}}-\underbrace{kD^{(\alpha)}(S)\nabla S}_{\begin{subarray}{c}\text{capillary}\\ \text{term}\end{subarray}}\Bigg{)}=q_{w},}$$

$$D^{(\alpha )}(S):=-\frac{{\lambda }_w{\lambda }_n}{{\lambda }_T}\frac{\partial P_c^{(\alpha )}(S)}{\partial S}.$$

$$P_c^{(m)}(S):=p_e^{(m)}{S}^{-1/{\theta }^{(m)}}-p_e^{(m)}(1-S{)}^{-1/{\theta }^{(m)}},$$

$$P_c^{(f)}(S):=p_{\text{max}}^{(f)}(1-S),$$

$$\sum _{\alpha \in \{m,f\}}{\mathbf{u}}_{\ell }\cdot {\mathbf{n}}_{\alpha }=0\text{along}{\Gamma }_{mf}\text{for}\ell \in \{w,n\},$$

$$\widetilde{P}^{(\beta)}_{c}(S^{(\beta)}):=\min\big{(}P_{c,\textit{max}}^{(\alpha)},\max\big{(}P^{(\beta)}_{c}(S^{(\beta)}),P_{c,\textit{min}}^{(\alpha)}\big{)}\big{)},\quad\forall\,S^{(\beta)}\in[0,1],\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta\in\{m,f\},\beta\neq\alpha,}$$

$${\tilde{P}}_c^{(m)}(S^{(m)})={\tilde{P}}_c^{(f)}(S^{(f)})\text{along}{\Gamma }_{mf}.$$

$$\widetilde{P}^{(m)}_{c}(S^{(m)})=\min\big{(}p_{\textit{max}}^{(f)},\max\big{(}P^{(m)}_{c}(S^{(m)}),0\big{)}\big{)},\qquad\text{and}\qquad\widetilde{P}^{(f)}_{c}(S^{(f)})=P^{(f)}_{c}(S^{(f)}).$$

$$V_i{\varphi }_i\frac{S_{\ell ,i}^{n+1}-S_{\ell ,i}^n}{\Delta t}+\sum _{j\in \text{adj}(i)}{F}_{\ell ,ij}^{n+1}=V_i{q}_{\ell ,i},\forall i\in \{1,\dots ,N\}\forall \ell \in \{w,n\},$$

$$T_{ij}:=T_{ji}:=(\widehat{T}^{-1}_{ij}+\widehat{T}^{-1}_{ji})^{-1}:=\bigg{(}{\displaystyle\frac{d_{i,(ij)}}{k_{i}\mathcal{A}_{ij}}+\displaystyle\frac{d_{j,(ji)}}{k_{j}\mathcal{A}_{ji}}}\bigg{)}^{-1},$$

$$F_{\ell ,ij}^{n+1}:=T_{ij}{\lambda }_{\ell ,ij}^{\text{P}}\Delta {\Phi }_{\ell ,ij},$$

$$\lambda_{\ell,ij}^{P}:=\left\{\begin{array}[]{ll}\lambda_{\ell}(S_{i})&\text{if}\quad\Delta\Phi_{\ell,ij}\geq 0\\ \lambda_{\ell}(S_{j})&\text{otherwise,}\end{array}\right.$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta\Phi_{w,ij}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta\Phi_{n,ij}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}F_{\ell,ij}^{n+1}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}:=\frac{\lambda_{\ell,ij}^{P}}{\lambda_{w,ij}^{P}+\lambda_{n,ij}^{P}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{u}_{T,ij}}+T_{ij}\frac{\lambda_{\ell,ij}^{P}\lambda_{m,ij}^{P}}{\lambda_{w,ij}^{P}+\lambda_{n,ij}^{P}}(\rho_{m}-\rho_{\ell})g\Delta z_{ij}+T_{ij}\frac{\lambda_{\ell,ij}^{P}\lambda_{m,ij}^{P}}{\lambda_{w,ij}^{P}+\lambda_{n,ij}^{P}}(\Delta P_{c,m}-\Delta P_{c,\ell})_{ij},}$$

$$V_{\ell,ij}^{n+1}:=\frac{\lambda_{\ell,ij}^{V}}{\lambda_{w,ij}^{V}+\lambda_{n,ij}^{V}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{u}_{T,ij}},$$

$$\lambda_{\ell,ij}^{V}(S_{i},S_{j}):=\left\{\begin{array}[]{ll}\lambda_{\ell}(S_{i})&\text{if}\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{u}_{T,ij}}>0\\ \lambda_{\ell}(S_{j})&\text{otherwise}.\\ \end{array}\right.$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G_{\ell,ij}^{n+1}:=T_{ij}\frac{\lambda_{\ell,ij}^{G}\lambda_{m,ij}^{G}}{\lambda_{w,ij}^{G}+\lambda_{n,ij}^{G}}(\rho_{m}-\rho_{\ell})g\Delta z_{ij}\quad\text{for}\,\ell\neq m,}$$

$$\lambda_{\ell,ij}^{G}(S_{i},S_{j}):=\left\{\begin{array}[]{ll}\lambda_{\ell}(S_{i})&\text{if}\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\rho_{m}-\rho_{\ell})g\Delta z_{ij}>0}\\ \lambda_{\ell}(S_{j})&\text{otherwise}.\\ \end{array}\right.$$

$$C_{\ell ,ij}^{n+1}:=T_{ij}{D}_{ij}^{(\alpha )}\Delta S_{\ell ,ij},$$

$$D_{ij}^{(\alpha )}(S_i,S_j):=\underset{S\in [\min (S_i,S_j),\max (S_i,S_j)]}{\max }{D}^{(\alpha )}(S).$$

$$\left\{\begin{array}[]{ll}F_{w,ij}^{n+1}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{u}_{T,ij}},S_{i},S_{ij}^{(\alpha)})=-F_{w,ji}^{n+1}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\overline{u}_{T,ij}},S_{j},S_{ij}^{(\beta)})\\[5.0pt] \widetilde{P}^{(\alpha)}_{c}(S_{ij}^{(\alpha)})=\widetilde{P}^{(\beta)}_{c}(S_{ij}^{(\beta)}).\\ \end{array}\right.$$

$$\displaystyle F_{\ell,ij}^{n+1}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{u}_{T,ij}},S_{i},S_{ij}^{(\alpha)}):=V_{\ell,ij}^{n+1}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{u}_{T,ij}},S_{i},S_{ij}^{(\alpha)})+G_{\ell,ij}^{n+1}(S_{i},S_{ij}^{(\alpha)})+C_{\ell,ij}^{n+1}(S_{i},S_{ij}^{(\alpha)}),$$

$$S^{(f)}_{ij}=h(S^{(m)}_{ij}):=(P^{(f)}_{c})^{-1}\big{(}\widetilde{P}^{(m)}_{c}(S^{(m)}_{ij})\big{)},$$

Full text

1 Introduction

Reservoir simulation is an essential tool for a comprehensive understanding of both petroleum recovery and groundwater processes. The modeling of multiphase flow and transport in porous media entails the difficult task of solving the governing partial differential equations (PDEs). In heterogeneous porous media, a major modeling challenge stems from the fact that the saturation-dependent flux coefficients representing the relative permeabilities and capillary pressures are spatially discontinuous between different rock types, leading to complex capillary-driven flow dynamics characterized by discontinuous pressure and saturation profiles.

These discontinuities arise because capillary pressure depends on saturation but also on intrinsic rock properties exhibiting large and abrupt spatial variations, such as the pore throat size and wettability. These contrasts in intrinsic rock properties and saturation functions are particularly severe in fractured porous media. For these problems, the wetting phase propagating into the fracture quickly invades the lower-permeability matrix driven by a large capillary pressure gradient, while the non-wetting phase is displaced into the fractures and can be produced. This counter-current imbibition process at the matrix-fracture interface is a key oil recovery mechanism that must be resolved accurately to obtain a reliable prediction of the flow dynamics.

Therefore, it is essential to design with an accurate and specialized numerical scheme and solution strategy able to capture abrupt spatial jumps in the relative permeability and capillary pressure functions. The accuracy of standard flux approximations applied to fractured systems is often limited by the presence of large cells in the matrix next to the fracture, since this configuration often leads to large errors in the flux computation. This issue is particularly severe for Discrete Fracture Model (DFM) simulations (Karimi-Fard et al.,, 2003) and dual-permeability models. Therefore, we particularly seek robust numerical schemes that produce physically-consistent predictions of the capillary-driven imbibition rate in the context of fractured media on coarse grids.

Following the work of Van Duijn et al., (1995), interface conditions have been used in the finite-volume method to improve the accuracy of the flux in the presence of capillary heterogeneity (Enchéry et al.,, 2006; Cancès,, 2009; Evje and Friis,, 2012; Brenner et al.,, 2013; Hamon et al.,, 2018; Ahmed et al.,, 2018). In this approach, the flux is computed by solving a local nonlinear system that enforces a capillary equilibrium at the interface. Along with the interface conditions, the introduction of additional interface-based degrees of freedom makes it possible to resolve complex interfacial flow dynamics that would otherwise require extreme spatial refinement. Interface conditions have initially been applied to ensure physically consistent simulations of capillary trapping taking place during drainage in the context of DNAPL trapping (Niessner et al.,, 2005; Papafotiou et al.,, 2010), and more recently, to represent capillary imbibition in matrix-fracture systems (Brenner et al.,, 2017, 2018; Aghili et al.,, 2019). For completeness, we mention here the alternative numerical approaches to handle capillary approaches in the mixed finite-element method (Hoteit and Firoozabadi,, 2005) and in discontinuous Galerkin schemes (Evje and Friis,, 2012; Arbogast et al.,, 2013).

In this work, we focus on a widely used fully implicit finite-volume method with cell-centered variables to demonstrate that interface conditions can be used in existing reservoir simulators to improve the accuracy of oil recovery predictions. We show that, in the context of high-contrast matrix-fracture systems with capillary heterogeneity, a straightforward modification to the standard flux computations to better account for the local capillary equilibrium at the interface results in a more accurate imbibition rate.

To do that, we extend the analysis performed in Hamon et al., (2018) to demonstrate the usefulness of the discrete interface conditions for the simulation of capillary-dominated counter-current imbibition. We perform a detailed analysis of the impact of the upwinding of the saturation-dependent coefficients of the flux on the accuracy of a fully implicit finite-volume scheme based on the interface conditions. We consider both the standard Phase-Potential Upwinding (PPU) and the Implicit Hybrid Upwinding (IHU). The latter was initially developed to improve the nonlinear solver convergence for two-phase transport with buoyancy but has been extended to coupled flow and transport with three phases and with capillarity (Hamon and Tchelepi,, 2016; Hamon et al.,, 2016; Lee and Efendiev,, 2016; Hamon et al.,, 2018; Lee and Efendiev,, 2018; Moncorgé et al.,, 2019).

In the first numerical example, we simulate capillary-driven spontaneous imbibition. We demonstrate that the use of interface conditions improves the accuracy of the interfacial flux on practical (coarse) grids, which results in a more accurate prediction of the imbibition rate. We also note that IHU-C, combining the IHU scheme with interface conditions, is more accurate than the PPU-based schemes. Then, we consider the capillary interactions with viscous and buoyancy forces in a forced imbibition case, and focus on the ability of the schemes to capture non-wetting phase trapping in the matrix due to the capillary end-effect. The conclusions are supported by a truncation error analysis.

In the remainder of the paper, we first review the governing equations and the interface conditions in Section 2. Then, we present the fully implicit finite-volume scheme in Section 3, with an emphasis on the approximation of the saturation-dependent coefficients in PPU or IHU. After analyzing the truncation error of the schemes for the homogeneous case in Section 4, we present numerical examples in Section 5.

2 Mathematical model

2.1 Governing equations

We consider two-phase flow in an incompressible porous medium. The incompressible wetting and non-wetting phases are denoted by the subscripts $w$ and $n$, respectively. The mass conservation laws read

[TABLE]

where $\phi(\boldsymbol{x})$ is the porosity of the medium, $S_{\ell}(\boldsymbol{x},t)$ is the saturation of the phase $\ell$, and $t$ is the time. The source/sink term $q_{\ell}$ is positive for injection and negative for production. The phase velocity, $\boldsymbol{u}_{\ell}$, is defined using the multiphase extension of Darcy’s law as

[TABLE]

Here, $k(\boldsymbol{x})$ is the scalar absolute permeability of the medium, $\lambda_{\ell}$ is the phase mobility defined as $k_{r\ell}/\mu_{\ell}$, where $k_{r\ell}(\boldsymbol{x},S_{\ell})$ is the phase relative permeability, $\mu_{\ell}$ is the phase viscosity, and $\rho_{\ell}$ is the phase density. We assume that the relative permeabilities are strictly increasing functions of their own saturation. In this work, we set $\rho_{n}=800\,\text{kg}.\text{m}^{-3}$ and $\rho_{w}=1000\,\text{kg}.\text{m}^{-3}$. We assume the same viscosity $\mu_{w}=\mu_{\textit{n}}=1\,\text{cP}$ for both phases. The gravitational acceleration is denoted by $g$, and the depth by $z$ (positive going downward). Even though studying buoyancy effects is not the main focus of this work, the gravity term is included in (2) for completeness. The saturations of both phases are constrained by the following equation:

[TABLE]

To account for capillary heterogeneity, we divide the domain $\Omega$ into non-overlapping subdomains $\Omega_{\alpha}$ referred to as rock regions such that $\cup_{\alpha}\,\Omega_{\alpha}=\Omega$. That is, capillary pressure and relative permeabilities are functions of saturation that can be spatially discontinuous between rock regions. In each rock region, the capillary pressure constraint relates the two phase pressures with

[TABLE]

We assume that capillary pressure is a strictly decreasing function of the wetting-phase saturation. Capillary pressure functions that do not satisfy this property require a specific treatment in the numerical scheme (see Brenner et al., (2017)) and will be considered in future work.

2.2 Pressure and transport equations

In this work, we consider two types of upwinding methods in the flux approximation required by the finite-volume scheme. Although the Phase-Potential Upwinding (PPU) methodology is applied directly to (1), the Implicit Hybrid Upwinding (IHU) approach is constructed as the approximation of the flux in fractional flow form. This formulation is obtained by decomposing the governing equations into a flow problem and a transport problem for one of the two phases. To obtain this decomposition, we first sum the phase velocities to define the total velocity as

[TABLE]

where $\lambda_{T}$ is the total mobility, defined as the sum of the phase mobilities. The total velocity is a function of space as well as of the reference pressure chosen to be $p=p_{n}$ and the reference saturation denoted by $S=S_{w}$. Summing (1) and then using the saturation constraint results in the elliptic equation governing the temporal evolution of the pressure field:

[TABLE]

Next, the highly nonlinear parabolic transport equation is obtained by using (5) to formally eliminate the pressure variable from (1), yielding

[TABLE]

where we have linearized the capillary pressure term of (2) to introduce a nonlinear capillary diffusion coefficient defined as

[TABLE]

The capillary diffusion coefficient is used to construct a differentiable and bounded capillary flux in the IHU approach. The function $S\mapsto D^{(\alpha)}(S)$ depends on space and saturation, is nonnegative, and is equal to zero only at the saturation endpoints. The maximum value of the capillary diffusion coefficient is denoted by $D^{(\alpha)}_{\textit{max}}$. We consider two rock regions, the matrix and the fracture, but the methodology is applicable to any number of regions. We refer to the matrix region with the superscript $\alpha=m$, and to the fracture region with the superscript $\alpha=f$.

In the matrix, we rely on the capillary pressure model of Skjaeveland et al., (1998). Capillary pressure is written as a function of saturation as

[TABLE]

where $p_{e}^{(m)}$ is the entry pressure in the matrix, and $\theta^{(m)}$ determines the saturation exponent. In this work, we set $\theta^{(m)}=4$. This type of capillary pressure function – used to fit experimental curves in Masalmeh et al., (2005) and described in the context of numerical studies in Tavassoli et al., (2005); Schmid and Geiger, (2013) – is representative of a capillary imbibition process in a mixed-wet system. The wetting-phase saturation $S^{\star}$ that satisfies $P_{c}(S^{\star})=0$ is used to distinguish two imbibition regimes – namely, a spontaneous imbibition regime in which the capillary pressure function is positive ($S<S^{\star}$), and a forced imbibition regime in which the capillary pressure function is negative ($S>S^{\star}$).

We highlight that the capillary pressure curve obtained in (9) is unbounded when $S$ is close to zero or one. Therefore, to avoid numerical problems with the standard Phase-Potential Upwinding (PPU) scheme with this analytical capillary pressure, we modify the formulation given by (9) in a preprocessing step. This procedure aims at enforcing upper and lower bounds on the capillary pressure in the matrix while preserving its monotonicity and differentiability with respect to saturation. The details are given in A. In the fracture, we use the following linear capillary pressure model:

[TABLE]

where the maximum capillary pressure in the fracture is a constant denoted by $p_{\textit{max}}^{(f)}>0$. The relative permeabilities in the fracture are linear. Figure 1 shows the capillary pressure curves used in this work as well as the corresponding capillary diffusion coefficients.

2.3 Interface conditions for local capillary equilibrium

In heterogeneous porous media, interface conditions can be used to enforce a local capillary equilibrium at the interfaces between rock regions with different capillary pressure functions (Van Duijn et al.,, 1995). Once discretized in the numerical scheme, these conditions provide a methodology to better capture the capillary jump in the computation of the flux. Here, we consider these conditions in the context of matrix-fracture systems. Considering an interface $\Gamma_{mf}$ between the matrix, $\Omega_{m}$, and the fracture, $\Omega_{f}$, the two conditions imposed at the discontinuity are given below.

Mass flow rate conservation condition at the discontinuity

The first condition ensures that the mass flow rate of each phase is conserved at the interface between rock regions $\Omega_{m}$ and $\Omega_{f}$. It reads

[TABLE]

where $\boldsymbol{n}_{\alpha}$ is the outward normal to the interface with respect to each rock region.

Extended pressure condition at the discontinuity

We first define the extended pressure as

[TABLE]

where $P^{(\alpha)}_{c,\textit{min}}$ and $P^{(\alpha)}_{c,\textit{max}}$ denote, respectively, the minimum and the maximum of the capillary pressure function on the saturation interval. Assuming that the capillary pressure curves have been bounded using the procedure described in A and therefore have finite endpoints, the quantities $P^{(\alpha)}_{c,\textit{min}}$ and $P^{(\alpha)}_{c,\textit{max}}$ are well defined. Using that, the extended pressure condition enforces a local capillary equilibrium at the interface with the requirement that

[TABLE]

The left-hand side and right-hand side of (13) reduce to, respectively,

[TABLE]

The extended pressure condition states that although the capillary pressure functions are different in $\Omega_{f}$ and in $\Omega_{m}$, the phase pressures are still continuous at the interface whenever S^{(m)}\in[(P^{(m)}_{c})^{-1}\big{(}p^{(f)}_{\textit{max}}\big{)}, (P^{(m)}_{c})^{-1}\big{(}0\big{)}]. This condition limits the possible saturation pairs satisfying this local capillary equilibrium if S^{(m)}\in[(P^{(m)}_{c})^{-1}\big{(}p^{(f)}_{\textit{max}}\big{)}, (P^{(m)}_{c})^{-1}\big{(}0\big{)}], and otherwise enforces a physically consistent discontinuity in the saturation profile at the interface. Figure 2 shows example of the capillary equilibrium in a matrix-fracture system.

3 Fully

implicit finite-volume scheme

The domain is discretized into $N$ cells, with $N^{(m)}$ (respectively, $N^{(f)}$) cells in the matrix (respectively, in the fracture). A fully implicit finite-volume discretization of (1) results in the following equation:

[TABLE]

where $F_{\ell,ij}^{n+1}$ is the numerical flux of phase $\ell$ for interface $(ij)$ between two cells $i$ and $j$, $adj(i)$ is the set of neighbors of cell $i$, $V_{i}$ is the bulk volume of the cell $i$, and $\Delta t$ is the time step size.

We compute the numerical flux with a two-point flux approximation (TPFA) and focus on one-dimensional uniform grids for simplicity. We decompose the numerical flux into a static rock and geometric transmissibility coefficient independent of the primary variables, $T_{ij}$, and a dynamic phase-dependent part that is a function of the primary variables. The static term, $T_{ij}$, is precomputed using a harmonic average of one-sided transmissibilities, denoted by $\widehat{T}_{ij}$, as

[TABLE]

where $k_{i}$ (respectively, $k_{j}$) is the absolute permeability in cell $i$ (respectively, $j$), $d_{i,(ij)}$ (respectively, $d_{j,(ji)}$) is the distance from the center of cell $i$ (respectively, $j$) to the interface, and $\mathcal{A}_{ij}=\mathcal{A}_{ji}$ is the area of the interface. For the simple grids considered in this work, we have $d_{i,(ij)}=d_{j,(ji)}$. The dynamic term at the interface is computed using a first-order upwinding scheme relying on cell-averaged quantities. The two upwinding schemes used in this work are discussed next.

3.1

Phase-Potential Upwinding

In Phase-Potential Upwinding (PPU, Aziz and Settari,, 1979), the numerical flux is constructed as an approximation of the interfacial integral of the phase velocity written with Darcy’s law:

[TABLE]

where the superscript $P$ denotes the PPU scheme. The phase mobility is approximated at interface ($ij$) using the upstream saturation with respect to the discrete phase-potential difference as

[TABLE]

where the discrete phase-potential difference involves a viscous, gravitational, and capillary contribution:

[TABLE]

In fractional flow form, the PPU numerical flux is decomposed into a viscous flux, $V_{\ell,ij}^{n+1}$, a buoyancy flux, $G_{\ell,ij}^{n+1}$, and a capillary flux, $C_{\ell,ij}^{n+1}$. Specifically, for $m\neq\ell$,

[TABLE]

where $\Delta P_{c,w}:=\Delta P_{c}$ and $\Delta P_{c,n}:=0$. This methodology yields a monotone numerical flux that is not differentiable in the presence of counter-current flow due to buoyancy or capillarity (Sammon,, 1988; Brenier and Jaffré,, 1991; Kwok and Tchelepi,, 2008). Previous authors identified this non-differentiability as the cause of nonlinear convergence difficulties for Newton solvers when large time steps are taken (Wang and Tchelepi,, 2013; Li and Tchelepi,, 2015). This limitation has motivated the design of a robust Implicit Hybrid Upwinding detailed in the next section.

In heterogeneous media, the standard PPU scheme presented above is employed for all interfaces without special treatment for the interfaces between two rock regions with different capillary pressure functions. Therefore, this formulation does not strictly enforce the interface conditions at the interfaces, and is referred to as standard PPU in the remainder of this paper. But, in addition to this standard PPU scheme, we also consider the combination of PPU with discrete interface conditions to better resolve the capillary jumps. This extended formulation is discussed in Section 3.3 and is referred to as PPU-C.

3.2 Implicit Hybrid Upwinding

In the Implicit Hybrid Upwinding (IHU) strategy, the viscous, buoyancy, and capillary fluxes are upwinded separately to achieve a differentiable and bounded flux for the transport problem and therefore improve nonlinear convergence. This robust approach was originally developed for two-phase transport with buoyancy (Eymard et al.,, 1989; Lee et al.,, 2015), and was later extended to three-phase flow and transport with capillary pressure (Hamon et al.,, 2016; Lee and Efendiev,, 2016; Hamon et al.,, 2018; Lee and Efendiev,, 2018; Moncorgé et al.,, 2019) and to a multi-point upwinding strategy reducing grid orientation effects (Hamon and Mallison,, 2020). The scheme presented in Hamon et al., (2018) was shown to lead to a more robust Newton convergence than with the standard PPU approach for challenging simulations with strong buoyancy and capillary effects. It is briefly reviewed below for the two-phase case.

The IHU numerical flux is constructed as the approximation of the interfacial integral of the phase velocity written in fractional flow form. The upwinding for the viscous flux is based on the direction of the total flux at interface ($ij$) denoted by $\overline{u}_{T,ij}$ as

[TABLE]

where

[TABLE]

For simplicity, we evaluate the total flux, $\overline{u}_{T,ij}$, using PPU in our 1D examples. However, in multiple dimensions, a smooth weighted averaging for the total flux is advantageous in some cases (Hamon et al.,, 2016, 2018).

Although only mild buoyancy effects are considered in the numerical examples, we describe the discretization of the IHU gravity flux since the improved nonlinear behavior obtained with IHU can offset the additional cost induced by the interface conditions (see Remark 3). The IHU gravity flux is upwinded based on the density difference between phases at interface ($ij$) as

[TABLE]

where

[TABLE]

In the IHU numerical scheme, the discretization of the capillary flux relies on the nonlinear capillary diffusion coefficient defined in (8). This term is evaluated at interface $(ij)$ using the maximum value of the diffusion coefficient between the saturations in cells $i$ and $j$ as follows:

[TABLE]

where

[TABLE]

We emphasize that $D^{(\alpha)}_{ij}$ is not a fixed quantity but depends on $S_{i}$ and $S_{j}$. The capillary flux is differentiable and monotone with respect to saturation. By construction, this flux term is bounded in the entire saturation range because $D^{(\alpha)}_{ij}(S_{i},S_{j})$ is bounded and vanishes at both saturation endpoints. This is in contrast to the PPU capillary flux that is unbounded when the capillary curve has an asymptote.

Remark 1**.**

As explained above, the IHU scheme yields a bounded capillary flux and therefore does not require any modification to the capillary pressure function given in (9). However, the PPU scheme produces an unbounded capillary flux, which creates significant numerical difficulties in the simulations. To overcome this issue, we modify (9) using the methodology detailed in A. Although this modification only aims at improving the behavior of the PPU scheme, we use it for both PPU and IHU to perform a fair comparison between the two schemes.

We have assumed in this section that both cells are in the same rock region denoted by $\Omega_{\alpha}$. When the cells belong to two rock regions with different capillary pressure functions, we employ a specific treatment that enforces an equilibrium at the interface to compute the capillary flux in heterogeneous porous media. The incorporation of the discrete interface condition in the scheme, referred to as IHU-C, is presented in Section 3.3.

3.3 Discrete interface conditions

In this section, we review the discretization of the interface conditions that can be used to enforce a local capillary equilibrium at the interface between the fracture and the matrix. The formulation follows the work of Van Duijn et al., (1995) and the discretization employed in Cancès, (2009); Hamon et al., (2018). We introduce two degrees of freedom located at each interface between two rock regions with different capillary pressure functions. Specifically, $S^{(m)}_{ij}$ and $S^{(f)}_{ij}$ are the saturations located at the interface on the side of the matrix and on the side of the fracture, respectively. The numerical flux is then computed by finding the saturation pair, $(S^{(m)}_{ij},S^{(f)}_{ij})$, that satisfies the following two-by-two nonlinear system for a fixed pair of cell-centered saturations $(S_{i},S_{j})$ and a fixed total flux {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{u}_{T,ij}}:

[TABLE]

The first constraint in (28) guarantees discrete mass flow rate conservation of the wetting phase at the interface. Due to the fixed total flux, equation (28) also yields conservation of mass flow rate for the non-wetting phase since we have F^{n+1}_{\textit{n},ij}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{u}_{T,ij}}-F^{n+1}_{w,ij}. The second constraint enforces the equality of the extended capillary pressures at the interface. In the mass conservation constraint, the fluxes are computed as

[TABLE]

where the third argument of the flux is evaluated with the interface saturation in the same rock region.

We assumed in Section 2 that capillary pressure is a strictly decreasing function of saturation. As a result, the interface conditions yield a unique solution pair whenever the flux approximation used to compute (29) is monotone with respect to saturation. Therefore, as in the homogeneous case, we can consider two different methods to discretize this flux. In the scheme referred to as PPU-C, we use PPU to evaluate the fluxes in the first line of (28). Instead, in the IHU-C scheme, IHU is applied to discretize these fluxes. These fluxes are computed with the methodologies described in Sections 3.1 and 3.2, except that each flux is based on a one-sided transmissibility, $\widehat{T}_{ij}$, and a half depth difference, $\widehat{\Delta z}_{ij}$.

We summarize the numerical procedure employed to solve (28) in Algorithm 1. In our current investigation, we can choose $\Omega_{\alpha}$ to be the matrix since $P_{c,\textit{max}}^{(m)}>p^{(f)}_{\textit{max}}$ and $P_{c,\textit{min}}^{(m)}<p^{(f)}_{\textit{min}}=0$. Then, we use the extended pressure condition stated in the second line of (28) to write

[TABLE]

as a function of $S^{(m)}_{ij}$. We employ Newton’s method to solve the scalar mass conservation equation at the interface using $S^{(m)}_{ij}$ as primary variable. The general case can be handled by switching the saturation variable so that the function $h$ in the local nonlinear system remains well defined whenever one of the phase pressures is discontinuous at the interface (Alali,, 2018).

The derivatives computed in the last two lines of Algorithm 1 are then used when the flux terms are assembled into the global Jacobian matrix. The implicit function theorem is used to compute the derivatives of the interfacial saturation, $S^{(\alpha)}_{ij}$, with respect to the cell-centered variables, $\tau_{i}$ and $\tau_{j}$ (B). Noting that the local nonlinear systems are only solved up to a user-defined tolerance, we have observed that the performance of the global nonlinear solver is better when the left flux (respectively, right flux) is used to compute the derivatives of the variables in the left cell (respectively, right cell).

Remark 2**.**

Our formulation yields a unique solution when capillary pressure is a strictly decreasing saturation function. Extending the formulation to capillary pressure functions that are not strictly decreasing requires modifying (28) as in Brenner et al., (2013, 2017) and will be considered in future work.

Remark 3**.**

Interface conditions are used to improve the accuracy of the scheme at the interfaces between different rock types on coarse grids. As shown in Section 5, achieving the same accuracy without interface conditions would require to significantly refine the mesh in the neighborhood of these interfaces, which, in 3D, becomes extremely expensive. Instead, although the interface conditions require additional computations, they only involve local, single-variable solves at the mesh interfaces where capillary pressure is discontinuous in space. These local systems are eliminated during the assembly and therefore do not need to be introduced in the global Jacobian system. In addition, we propose to combine interface conditions with a hybrid discretization of the flux given in Section 3.2. This hybrid flux discretization yields a robust scheme that is particularly well suited for problems with strong buoyancy and capillary forces. As shown in Brenner et al., (2020) on challenging, three-dimensional problems, the improved global nonlinear behavior achieved with the hybrid scheme can offset the additional cost induced by the local interface conditions.

4 Numerical flux properties

We first compare the PPU and IHU numerical fluxes with viscous and capillary forces in homogeneous one-dimensional media. We consider the fluid properties of the matrix, i.e., quadratic relative permeabilities and the capillary pressure of Fig. 1 defined by $p^{(m)}_{e}=3\,\text{psi}$ and $\theta^{(m)}=4$. In this section, we assume that $\overline{u}_{T}$ is constant in space and time. The behavior of the schemes for heterogeneous test cases will be studied in Section 5.

We first discuss the pure capillary case ($\overline{u}_{T}=0$) shown in Fig. 3. For saturations close to the endpoints, 0 and 1, we observe that the PPU scheme overestimates the absolute magnitude of the capillary force since the numerical flux tends to infinity. This difference stems from the fact that the PPU capillary flux is constructed by approximating directly the gradient of capillary pressure. Conversely, the IHU numerical flux is based on a linearized form involving the capillary diffusion coefficient (8) and remains bounded in the entire saturation range. We show in the spontaneous imbibition test of Section 5.1 that IHU produces more accurate predictions of the imbibition rate on coarse grids than with PPU.

We now consider the viscous-capillary case with $\overline{u}_{T}=1/2$ shown in Figs. 4 and 5. We see in Figs. 4 and 5 that the PPU approach still overestimates the capillary force, yielding a large numerical flux for $S=0$ and $S=1$ while the numerical flux remains bounded with the IHU scheme. But, considering now Fig. 5, we observe that the two schemes produce different numerical fluxes close to the center of the saturation range (for $S_{L}\in[0.3,0.7]$ and $S_{R}\in[0.3,0.7]$). This is because the split treatment of the mobilities in IHU alters the viscous-capillary balance of forces and produces a different cocurrent flow region as shown in Fig. 5. The PPU and IHU cocurrent flow regions contain the saturation pairs $(S_{L},S_{R})$ such that $F^{\textit{PPU}}_{w}(\overline{u}_{T},S_{L},S_{R})F^{\textit{PPU}}_{n}(\overline{u}_{T},S_{L},S_{R})>0$ and $F^{\textit{IHU}}_{w}(\overline{u}_{T},S_{L},S_{R})F^{\textit{IHU}}_{n}(\overline{u}_{T},S_{L},S_{R})>0$, respectively. We quantify this difference with a truncation error analysis and show in Section 5.2 that this leads to different predictions of non-wetting phase trapping with PPU and IHU.

Next, we derive the truncation error for the IHU and PPU schemes. The calculation is done using the problem setup of Section 5.2 since we will apply this analysis to explain the respective behaviors of IHU and PPU for cocurrent flow in a medium undergoing forced imbibition. We assume one-dimensional cocurrent viscous-capillary flow with a (scalar) velocity $u_{T}=1/2$. We further assume that the saturation profile is monotone, since this is the case in the homogeneous region considered in Section 5.2. Let the wetting-phase saturation $x\mapsto S(x)$ be a twice differentiable function of space, with $S_{i}=S(x_{i})$. For cocurrent flow, the viscous truncation error is by construction the same for the two upwindings:

[TABLE]

where $u_{T}$ is the scalar total velocity and $M_{w}$ is defined as

[TABLE]

Next, considering the case $S_{max}<S_{i+1}\leq S_{i}\leq S_{i-1}$, where $S_{max}$ is the saturation for which the capillary diffusion coefficient of the matrix attains its maximum, we use Taylor series expansions to obtain the following expression of the capillary truncation error:

[TABLE]

where $D$ is the capillary diffusion coefficient defined in (8) and $\Delta x:=x_{i+1}-x_{i}:=x_{i}-x_{i-1}$. The same expression holds for the case $S_{i-1}\leq S_{i}\leq S_{i+1}\leq S_{max}$. The IHU capillary truncation error term of (33) is always bounded which is in agreement with Fig. 4(b). With the PPU scheme, considering again the cocurrent flow case, we write the capillary truncation error term as

[TABLE]

We note that, in (34), the term between brackets in the right-hand side contains the second derivative of capillary pressure and is unbounded for small saturations. Therefore, for any $\epsilon_{1}>0$, there exists a $\eta_{1}>0$ sufficiently small such that $S<\eta_{1}$ implies $|\mathcal{E}^{P}_{C}|>|\mathcal{E}^{H}_{C}|+\epsilon_{1}$, leading to large errors in the flux computation with PPU when the wetting-phase saturation is close to zero. This unbounded term in the leading capillary truncation error term is responsible for the large values observed in Figs. 4(a) and 5.

But, since the capillary pressure function considered in this work has an inflection point, there is a saturation, $S^{\text{inflec}}$, such that $\frac{\partial^{2}P_{c}}{\partial S^{2}}(S^{\text{inflec}})=0$. Therefore, considering any small $\epsilon_{2}>0$, one can find a $\Delta x_{2}>0$ and $\eta_{2}>0$ sufficiently small such that $\Delta x<\Delta x_{2}$, $|S-S^{\text{inflec}}|<\epsilon_{2}$ implies $|\mathcal{E}^{P}_{C}-(-\mathcal{E}^{H}_{C})|<\epsilon_{2}$. This means that in a neighborhood of $S^{\text{inflec}}$, the capillary truncation errors produced by PPU and IHU have an opposite sign. Since the two schemes produce the same viscous truncation error, $\mathcal{E}^{P}_{V}=\mathcal{E}^{H}_{V}$, for the cocurrent case, the interaction between the capillary and viscous errors is going to affect the accuracy of the IHU and PPU schemes differently.

Specifically, when we apply this analysis in Section 5.2, we will see that the leading terms in $\mathcal{E}^{\textit{H}}_{V}$ and $\mathcal{E}^{\textit{H}}_{C}$ have the same sign, which means that the capillary truncation error of IHU tends to amplify the viscous truncation error. Instead, the leading terms in $\mathcal{E}^{P}_{V}$ and $\mathcal{E}^{P}_{C}$ have an opposite sign, and therefore the two error terms of the PPU scheme tend to cancel each other, yielding a smaller viscous-capillary truncation error than with IHU.

5 Numerical examples

In this section, we refer to the standard Phase-Potential Upwinding scheme as PPU. In PPU-C, the standard scheme is modified with the introduction of discrete interface conditions at the interface between the matrix and the fracture. The fluxes involved in the discrete interface conditions are also discretized with PPU. Finally, in IHU-C, Implicit Hybrid Upwinding is employed in the homogeneous regions as well as in the discretization of the fluxes involved in interface conditions. The global nonlinear systems are solved with Newton’s method. Global convergence is achieved when

[TABLE]

5.1 Spontaneous imbibition

We first study the accuracy of the numerical schemes in predicting the hydrocarbon recovery rate driven by capillary imbibition in a fractured porous medium to illustrate the importance of the interface conditions. For simplicity, we consider a one-dimensional horizontal matrix-fracture system of length $L=20\,\text{m}$ sketched in Fig. 6. There is no source or sink term such that capillarity is therefore the only driving force of the flow. The left half of the domain is occupied by the matrix whose uniform permeability is set to $k^{(m)}=1\,\text{mD}$. The right half represents the fracture and has a permeability of $k^{(f)}=10^{5}\,\text{mD}$.

The fracture pore volume is artificially amplified by two orders of magnitude compared to that of the matrix to simulate fast wetting-phase injection through the fracture. Our assumption is that in a multidimensional fracture network where counter-current imbibition occurs, the non-wetting phase that enters the fracture will quickly move away due to viscous or buoyancy forces. Hence, the wetting-phase saturation in the fracture should remain close to one. A similar setup was employed in the recent study of Vo et al., (2019).

We consider two sets of saturation-dependent properties in the matrix. In the first set, we use quadratic relative permeabilities, leading to the diffusion coefficient of Fig. 1. In the second set, we use cubic relative permeabilities and we obtain the diffusion coefficient of Fig. 1. The matrix is initially fully saturated with the non-wetting phase, and the fracture is fully saturated with the wetting phase. The saturation profiles obtained with the three numerical schemes until steady state for quadratic relative permeabilities are in Fig. 7.

We compute the recovery percentage at a given dimensionless time as the ratio of the cumulative amount of non-wetting phase that has left the matrix until this time over the cumulative amount at steady state. The recovery percentage is reported as a function of dimensionless diffusion time defined as

[TABLE]

The steady state is reached at $t_{D}\approx 0.8$. In Fig. 8, this numerical experiment is repeated for various levels of spatial refinement in the matrix, from a single cell in the matrix ($N^{(m)}=N^{(f)}=1$) to 128 cells in the matrix ($N^{(m)}=N^{(f)}=128$). In the refinement study, the grid is refined at the same rate in the matrix and in the fracture. The resolution given by $N^{(m)}=128$ is too fine for practical simulations but is used as a reference to compute the error in the recovery percentage as a function of time. We see that the standard PPU scheme significantly over-predicts the recovery percentage for small $N^{(m)}$ – i.e., for coarser levels of spatial refinement. For instance, for $N^{(m)}=1$, PPU predicts that the recovery reaches $80\%$ after $t_{D}=0.063$), compared to $t_{D}=0.171$ with the fine reference. This over-prediction is reduced upon spatial refinement but with a slow convergence to the reference solution as $N^{(m)}$ is increased.

The over-prediction of the recovery percentage is much less severe for the schemes based on interface conditions. For $N^{(m)}=1$, PPU-C predicts the $80\%$ recovery time at $t_{D}=0.126$, as opposed to $t_{D}=0.171$ in the reference solution. In addition, this over-prediction error remains lower than with the standard PPU as the grid is refined. Finally, the IHU-C scheme slightly under-predicts the recovery compared to PPU and PPU-C. For $N^{(m)}=1$, the IHU-C scheme predicts the $80\%$ recovery time at about $t_{D}=0.202$ as compared to the reference solution that predicts about $t_{D}=0.171$. This under-prediction is reduced with spatial refinement at a faster rate than with PPU and PPU-C. Specifically, IHU-C with only two cells in the matrix results in a very good match with the reference solution. With $N^{(m)}=2$ and $N^{(m)}=4$, the estimated times for the $80\%$-recovery are within $5\%$ of those of the reference solution.

To distinguish early-time dynamics from the behavior at steady state, Fig. 9 shows the error in the recovery percentage as a function of dimensionless time for different levels of spatial refinement. The error in the recovery percentage is computed as

[TABLE]

where $R^{n}_{N^{(m)}}$ (respectively, $R^{n}_{\textit{ref}}$) is the recovery percentage at time $n$ for $N^{(m)}$ cells in the matrix (respectively, the recovery percentage at time $n$ for the reference solution). We consider both quadratic and cubic relative permeabilities. Figure 9 illustrates that even though all the schemes predict the steady state solution correctly, the imbibition rate is better captured by the schemes that rely on interface conditions. In particular, we can see by comparing Figs. 9(a), 9(c), and 9(e) that the error obtained with IHU-C converges quickly as the number of cells is increased. For cubic relperms (Figs. 9(b), 9(d), and 9(f)) the error decreases at approximately the same rate for the three schemes but is maller with IHU-C.

This is confirmed by Fig. 10, which compares the recovery error integrated over time obtained with the three schemes for different levels of spatial refinement. We compare the $L_{\infty}$-norm of the error in the recovery percentage for the full simulation until steady state as

[TABLE]

We highlight two key observations on Fig. 10. First, the three schemes produce the same solution upon refinement in space. Second, on coarser grids, large errors are observed with the standard PPU. These errors are significantly reduced with the introduction of interface conditions in PPU-C and IHU-C. For both PPU-C and IHU-C, these results illustrate the importance of interface conditions to capture the flux at the matrix-fracture interface and accurately predict the rate of imbibition in the matrix.

We underline that these results expressed in dimensionless time can be used to extend our analysis to other matrix-fracture configurations with different characteristic variables, such as absolute permeability, porosity, and length of penetration. In simulations not included here for brevity, we considered a modified version of the test case discussed here by doubling the penetration length and quadrupling the absolute permeability in the matrix. The error in recovery for this modified configuration, expressed as a function of dimensionless time until $t_{D}=0.4$, is in very good agreement with those of Figs. 8 to 10.

Remark 4**.**

We have observed that reducing the volume of the fracture to let the non-wetting phase saturation increase in the fracture does not change the conclusions of this section. Therefore, the numerical examples performed with smaller fracture volumes are not included to the present paper for brevity.

5.2 Forced imbibition

We now extend the study of the accuracy and efficiency of the numerical schemes in a fractured porous medium by considering a forced imbibition process. We consider a one-dimensional domain with a length $L=20\,\text{m}$. The domain, sketched in Fig. 11, is tilted by an angle of 15 degrees updip from left to right. The two fracture regions are located in the left quarter and in the right quarter of the domain. The matrix occupies the middle of the domain. Unlike in Section 5.1, the modified problem has a source term and a sink term maintaining the wetting phase saturation close to one in the fractures. Therefore, we do not amplify the fracture pore volume and all the cells in the domain have the same volume.

The values of the absolute permeabilities, relative permeabilities, and capillary pressures in the matrix and the fractures are the same as in the previous case. Initially, the matrix is 50%-saturated with the non-wetting phase, and the fracture is 97%-saturated with the wetting phase. The wetting phase is injected from the left at a fixed rate, and a producer operates at a fixed pressure of $3000\,\text{psi}$ on the right. We are interested in the forced capillary imbibition process, so we make sure that viscous and buoyancy forces are such that the saturation in the matrix explores the negative section of the capillary pressure curve ($S\geq 0.5$). In this configuration, buoyancy displaces the non-wetting phase out of the matrix to the right of the domain, and the wetting-phase injection creates a total flow from left to right. During this process, the local capillary equilibrium at the interface between the matrix and the fracture on the right imposes a non-zero non-wetting phase saturation in the matrix. This capillary end-effect (Richardson et al.,, 1952; Huang and Honarpour,, 1998), visible in the saturation profiles of Fig. 12, leads to non-wetting phase trapping in the matrix at steady state. For the two spatial resolutions of Fig. 12, we observe that the non-wetting phase remains relatively large in the matrix even though it vanishes entirely in the fractures as a result of production. In particular, the non-wetting phase saturation at the right boundary of the matrix does not decrease over time and stays close to 0.5 due to the capillary end-effect.

For this test case, we observe that this interfacial saturation is relatively well captured by standard PPU on the coarse grid. We still see that the schemes based on interface conditions are more accurate than standard PPU at the beginning of the simulation (0.01 PVI), particularly at the boundary between the left fracture and the matrix. But, at later times (0.1 PVI and 1 PVI), the coarse saturation profiles yielded by PPU and PPU-C are in agreement. At steady state, we note a mismatch at $x_{D}\approx 0.6$ in the representation of the sharp gradient between the coarse IHU-C results and the coarse profiles obtained with the PPU-based results. We see in Fig. 13 that the mismatch in the coarse saturation profiles results in an over-estimation of the average non-wetting phase saturation remaining in the matrix at steady state with IHU-C. In Fig. 13, it also leads to an under-estimation of the amount of non-wetting phase that is produced from the fracture on the right.

Next, we use the truncation error analysis of Section 4 to explain the discrepancy between the PPU-based schemes and IHU-C for the forced imbibition process. To simplify the analysis, we omit buoyancy forces and only consider the viscous-capillary equilibrium. The other parameters remain unchanged, resulting in the coarse-grid saturation profiles of Fig. 14. We apply the truncation error analysis of Section 4 to the saturation profiles in the matrix ($x\in[0.55,0.75]$) where the difference between the PPU-C and IHU-C results are relatively large. The flow in this region is cocurrent from left to right, which corresponds to the configuration studied in Section 4. The leading terms of the viscous, capillary, and total truncation errors are shown in Fig. 15.

For the viscous term in Fig. 15(a), we see that the two schemes lead to the same truncation error since they are based on the same discretization in the presence of cocurrent flow. In Fig. 15(b), the PPU and IHU capillary truncation error terms computed in (33)-(34) have a similar absolute magnitude but an opposite sign. This is key to understand the behavior of the total error, obtained by summing the viscous and capillary error terms and shown in Fig. 15(c). We observe that the capillary term tends to reduce the viscous error with PPU, but instead amplifies the viscous error with IHU, yielding a larger total error with IHU than with PPU. This is particularly the case when the saturation gradient is very sharp in the neighborhood of the interface with the fracture on the right ($x\in[0.65,0.75]$). In other words, the discrepancy does not result from the treatment of the interfacial saturations, but instead arises because of the split treatment of the viscous-buoyancy-capillary flux in the homogeneous regions. This methodology introduces splitting errors that affect the accuracy of IHU-C in the presence of sharp saturation gradients. We plan to address this limitation with a high-resolution of IHU-C along the lines of Mykkeltvedt et al., (2017).

6 Conclusion

In this work, we have studied a fully implicit finite-volume scheme for capillary-dominated multiphase flow in matrix-fracture systems. In particular, we have studied the impact of discrete interface conditions enforcing a local capillary equilibrium between the matrix and the fracture on the accuracy of the predicted imbibition rate into the matrix. We have also considered the interaction of the discrete interface conditions with the upwinding schemes – namely, PPU and IHU – used to evaluate the interfacial saturation-dependent coefficients.

We demonstrate that both the PPU and IHU schemes based on discrete interface conditions (and referred to in the paper as PPU-C and IHU-C, respectively) yield an interfacial flux that is more accurate than with the standard scheme. Our results show that using interface conditions to compute the numerical flux between the matrix and the fracture improves the accuracy of the prediction of the imbibition rate compared to the standard scheme without interface conditions. Achieving the same accuracy with the standard scheme would require significant mesh refinement, which would dramatically increase the cost of large-scale, three-dimensional simulations. Future directions include extending these results for the three-phase case to further demonstrate the applicability of the scheme to practical simulation problems.

Appendix A Capillary pressure interpolation

Here, we explain the procedure used to enforce an upper bound, $P_{c,max}$, and a lower bound, $P_{c,min}$, on the capillary pressure function and therefore avoid infinite gradients when capillary pressure is accounted for. This is done by replacing the expression given in (9) with a polynomial function when the saturation is close to one of the endpoints, as follows:

[TABLE]

where $P_{c}$ is the function defined in (9). In a preprocessing step, we compute the triplets $(a,b,S^{-})$ and $(c,d,S^{+})$ such that the capillary pressure function defined in (39) is continuous and twice differentiable. Computing the first triplet involves solving the following system

[TABLE]

The nonlinear system solved to compute $(c,d,S^{+})$ is analogous. All the figures and numerical examples presented in this paper rely on (matrix) capillary pressures computed with (39). We set $P_{c,max}=15\,\text{psi}$ and $P_{c,min}=-15\,\text{psi}$.

Appendix B Derivatives of $S^{(\alpha)}_{ij}$ with respect to cell-centered variables

The last lines of Algorithm 1 require the computation of the derivatives of the saturation at interface $(ij)$, denoted by $S^{(\alpha)}_{ij}$, with respect to the cell-centered primary variables, $p_{i}$, $p_{j}$, $S_{i}$, and $S_{j}$. We remind the reader that the total flux, $\overline{u}_{T,ij}(\Delta p_{ij},S_{i},S_{j})$, is a function of $\Delta p_{ij}$, $S_{i}$, and $S_{j}$. Consider the function

[TABLE]

Given $\Delta p_{ij}$, $S_{i}$, and $S_{j}$, the local nonlinear solver finds the interfacial saturation, $S^{(\alpha)}_{ij}$, such that

[TABLE]

Noting that the monotonicity properties of the flux imply that

[TABLE]

we apply the implicit function theorem to compute the derivatives of $S^{(\alpha)}_{ij}=f_{ij}(\Delta p_{ij},S_{i},S_{j})$ with respect to $\Delta p_{ij}$, $S_{i}$, and $S_{j}$. For $\tau=\{a,b,c\}$, we obtain

[TABLE]