Modeling the effects of capillary pressure with the presence of full tensor permeability and discrete fracture models using the mimetic finite difference method

Abstract

Capillary dominated flow or imbibition—whether spontaneous or forced—is an important physical phenomena in understanding the behavior of naturally fractured water-driven reservoirs (NFR’s). When the water flows through the fractures, it imbibes into the matrix and pushes the oil out of the pores due to the difference in the capillary pressure. In this paper, we focus on modeling and quantifying the oil recovered from NFR’s through the imbibition processes using a novel fully implicit mimetic finite difference (MFD) approach coupled with discrete fracture/discrete matrix (DFDM) technique. The investigation is carried out in the light of different wetting states of the porous media (i.e., varying capillary pressure curves) and a full tensor representation of the permeability. The produced results proved the MFD to be robust in preserving the physics of the problem, and accurately mapping the flow path in the investigated domains. The wetting state of the rock affects greatly the oil recovery factors along with the orientation of the fractures and the principal direction of the permeability tensor. We can conclude that our novel MFD method can handle the fluid flow problems in discrete-fractured reservoirs. Future works will be focused on the extension of MFD method to more complex multi-physics simulations.

1 Introduction

Naturally fractured reservoirs (NFR’s) have been an important research topic in the field of multiphase flow in porous media due to the complex nature of the multi-phase flow in the fracture–matrix system. Such reservoirs form the vast majority of oil and gas reserves in the world (Saidi 1983; Gong and Rossen 2018) as they are a target for EOR applications Massarweh and Abushaikha (2020), and are characterized by their heterogeneous distribution of porosity and permeability Thomas et al. (1983). Thus, it is critical to understand the underlying physical phenomena that control the amount of recoverable hydrocarbons from NFR’s. A fractured reservoir is composed of two main domains: a matrix block and a fracture block. Typically, the fracture network is highly permeable compared to the matrix block that is characterized with high porosity, but other classification with different porosity/permeability ratios can exist Allan and Sun (2003). In such systems, fractures should be closely examined and their properties accurately evaluated to ensure a realistic representation of the reservoir during numerical simulations.

Fluid communication between the matrix and the fracture blocks in NFR’s is established through various phenomena: diffusion, advection, gravity drainage, and imbibition whether it's forced or spontaneous. The physical phenomena of spontaneous imbibition (SI) is generally of great interest since no external force is applied to induce the fluid exchange, and the oil production is determined by the mass transfer rate (Qasem et al. 2008; Hatiboglu and Babadagli 2007). The oil is carried through the high permeability conduit thus preventing the buildup of pressure differentials across the reservoir, and limiting the influence of viscous displacement. We call the production process in this case as capillary dominated flow. In simple words, spontaneous imbibition is defined as “the process by which a wetting fluid is drawn into a porous medium by capillary action” Morrow and Mason (2001). The main driving force in this process is the differential capillary pressure between the matrix and the fracture, leading to efficient imbibition of water into the matrix which in turn displaces oil into the fracture. The capillary pressure drives the wetting phase into the matrix, and thus the degree of wettability of the porous media greatly influences the ultimate recovery of the reservoir (Austad and Milter 1997; Alyafei 2019).

Spontaneous imbibition usually occurs in two main forms: counter-current and co-current flow. As the name imply, counter-current is when the wetting phase and the non-wetting phase flow in opposite directions resulting in a total net flow of zero. The boundaries are well sealed and closed except for the inlet boundary and no external pressure is applied. Is some cases, the non-wetting phase escapes from the inlet through which the water imbibes in which can be limited by applying a force equivalent to the capillary back-pressure. On the contrary, fluids flow in one direction from the inlet through the outlet with the other boundaries sealed in co-current imbibition. In water-wet reservoirs, counter-current imbibition is predominant and requires sensitive care when modeled. The work of Abd et al. (2019) highlights the importance of counter-current imbibition in hydrocarbon recovery and describes extensively all the major breakthroughs in the experimental, analytical and numerical aspects of this topic. Early efforts to understand the physical significance of the imbibition process focused on experimental work (Mattax and Kyte 1962; Iffly et al. 1972; Du Prey 1978; Hamon and Vidal 1986; Zhang et al. 1996), while subsequent efforts highlighted the key differences between the two modes of imbibition assisted oil recovery (co-current and counter current SI) and their unique characteristics (Bourbiaux and Kalaydjian 1990; Pooladi-Darvish and Firoozabadi 2000; Morrow and Mason 2001). These experimental efforts were later tested and verified using numerical techniques and semi-analytical solutions (Fischer and Morrow 2006; Schmid et al. 2011, 2016; Nooruddin and Blunt 2016; Khan et al. 2018; Abd and Alyafei 2018; Abd and Alyafie 2019) allowing for better interpretation of spontaneous imbibition processes (Fig. 1).

Fig. 1

Reproduced from Abd et al. (2019)

Different imbibition scenarios. a and b display geometries where there is a combination of counter- and co-current spontaneous imbibition. c and d show counter-current spontaneous imbibition.

Numerically, NFR’s can be represented by two distinctive models: dual porosity and discrete fracture/discrete matrix (DFDM) approach. In dual porosity models, the petrophysical properties of the matrix and the fracture are averaged, and a transfer term is utilized to establish fluid flow between the domains (Barenblatt et al. 1960; Warren and Root 1963; Abushaikha and Gosselin 2008). The main assumption in the development of this model is that a unified transfer rate is used implying that the fracture is filled with the wetting phase instantly Di Donato and Blunt (2004). However, recent modifications have been implemented to account for multiple transfer rates, and create a more realistic representation of the physical problem Geiger et al. (2013). On the other hand, the DFDM approach relies on discretizing the fracture and the matrix domain through numerical approximation methods, thus allowing to generate a simpler problem in terms of numerical complexity. Although this approach provides highly accurate results, it cannot be used for systems with complex network of interconnected fractures.

In the DFDM approach, the system can be discretized using various numerical schemes including finite element method (FEM) and finite volume method (FVM). In a matrix–fracture framework, FEM is generally used to discretize the domain, and is categorized into different families including classical FEM and Mixed-FEM (MHFE) methods. The main difference between FEM and MHFE is that the latter is locally conservative Abushaikha et al. (2017). This property allows for the extension of the method to account for anisotropic properties of the reservoir lithology. However, the resultant algebraic system is numerically challenging to solve, and thus mixed-hybrid finite element (MHFE) method was introduced. In MHFE method, Lagrange multipliers are used to account for the physical properties at the interface. This method has been tested for the effect of fractured media on the imbibition of water into the matrix in counter-current flow and proved to be accurate. However, MHFE is limited to some types of elements because of the evaluation of the MFE basis functions. Among these three methods, unlike MHFE and FVM, FEM is globally mass conservative and not locally, the property that rises from the intrinsic property of the method. For FVM, the local conservation is satisfied, however, especially for large-scale cases, the application of FVM may suffer from two mesh systems and constrained by memory consumption Bathe and Zhang (2002). MHFE is another locally mass conservative method Brezzi and Fortin (2012). It has the advantages of FEM and also naturally satisfies the locally mass conservation, which is suitable for the complex flow. However, in MHFE, the construction of the basis functions can be only applied to simple meshes due to the analytical evaluation of basis functions. If a simple and good grid is necessary for FVM and MHFE, this is not necessarily required for mimetic finite difference (MFD) method whose formulation resembles the mixed-hybrid formulation (Brezzi et al. 2005, 2005; Brezzi and Fortin 2012; Brezzi et al. 2006). To overcome these restrictions, mimetic finite difference (MFD) method was introduced to model highly unstructured polygons (Da Veiga et al. 2009; Lipnikov et al. 2014) based on numerical computation of the basis thus broadening its applicability. This gives the mimetic scheme an advantage for reservoir simulation applications (Lipnikov et al. 2008; Campbell and Shashkov 2001) that can handle complex unstructured grids. In fractured media where the mesh size varies considerably, and unstructured grids are required to capture the changes in the geology of the domain, MFD is accurate in predicting hydrocarbon recovery (Abushaikha and Terekhov 2018, 2020; Hjeij and Abushaikha 2019) even when chemical reactions are considered Abd and Abushaikha (2021) contrary to the other methods (Abushaikha et al. 2015; Abd and Abushaikha 2020b, a). Moreover, fractured porous media consists of matrix and fractures. The fracture distribution is complex and variable in these reservoirs with strong characters of anisotropy. This kind of anisotropic porous medium is usually both heterogeneous and anisotropic. Full permeability tensor is therefore needed in modeling flow in anisotropic medium. FV schemes and MHFE schemes are currently widely applied to solve the flow model with a full permeability tensor. As we mentioned above, there are some limitations in the application of these two methods. MFD method allows for the utilization of full tensor permeability in both the matrix and the fracture which maintaining the flexibly of polygonal grids, thus mimicking closely the fluid flow in the actual reservoir given that it is implemented properly.

In this work, we aim to study the flow kinetics of imbibition processes for different rock wetting states, using full tensor permeability and descritized fractures using the MFD method. The investigation of the impact of capillarity on multiphase flow for highly heterogeneous and anisotropic reservoirs is a novel approach, specifically with the application of unstructured grids using the MFD. The paper will discuss the governing equations and the physical problem formulation of the domain using MFD and DFDM approaches, then multiple cases will be examined and analyzed for different systems with different fracture structures and rock properties. The results will show that this newly developed and novel approach yields highly accurate simulations in modeling spontaneous imbibition processes for fractured reservoirs while preserving the integrity of the complex geometry and the heterogeneous petrophysical properties using MFD.

2 Governing Equations

The governing equations for the two phase flow of a fluid in porous media are given by the conservation of mass and Darcy’s law. It is assumed that the diffusion and dispersion forces are very small so they can be neglected. The continuity equation is expressed as follows:

$$\begin{aligned} \frac{\partial (\phi \rho _{\alpha } S_{\alpha })}{\partial t} + \nabla \cdot (\rho _{\alpha } {\mathbf {u}}_{\alpha }) =\rho _\alpha q_{\alpha },\qquad \alpha =w,o \end{aligned}$$

(1)

where $\phi$ is the porosity of the medium, $S_{\alpha }$ denotes the saturation of the phase, $\rho _{\alpha }$ is the density of the phase, ${\mathbf {u}}_{\alpha }$ is Darcy’s velocity, and $q_{\alpha }$ is the mass flow rate. Darcy’s velocity can be written for each phase as:

$$\begin{aligned} {\mathbf {u}}_{\alpha }= - \frac{k_{r \alpha }}{\mu _{\alpha }}{\mathbf {k}}(\nabla p_{\alpha } -\rho _{\alpha }g\nabla {\mathbf {D}}) \end{aligned}$$

(2)

where $k_{r \alpha }$ is the phase relative permeability, $\mu _{\alpha }$ is the phase viscosity, ${\mathbf {k}}$ is the intrinsic permeability tensor, $\rho _{\alpha }g$ is the vertical pressure gradient and $\nabla {\mathbf {D}}$ is the vertical depth.

We then introduce the local constraints to fully close and define the system, represented by the saturation and capillary pressure relationships.

$$\begin{aligned} S_{\rm w}+S_{\rm o}= & 1 \end{aligned}$$

(3)

$$\begin{aligned} p_{\rm c}(S_{\rm w})= & p_{\rm o}-p_{\rm w} \end{aligned}$$

(4)

The phase mobility is expressed by:

$$\begin{aligned} \lambda _{\alpha }(S_\alpha )=\frac{k_{r\alpha }(S_\alpha )}{\mu _\alpha }, \qquad \alpha =w,o \end{aligned}$$

(5)

2.1 Numerical Discretization Using MFD

The starting point for the MFD methods are the governing equations of mass balance and Darcy velocity described in Eqs. 1 and 2 earlier. In this method, we utilize a vectorial basis function $\mathbf {V_i}$ to descritize the momentum balance. These basis functions have the following properties:

• The vectorial functions have a flux of {1} at interface i and zero elsewhere.

• The divergence of the basis function in constant over an element.

• The velocity field ${\mathbf {u}}$ can be estimated using ${\mathbf {u}}=\sum ^{N_{\rm f}}_{i=1} \mathbf {V_i}Q_i$ where $Q_i$ is the flux at the interface.

The key element of the MFD method is the local inner product $W_i$. According to the MFD method, Darcy’s law (Eq. 2) over a mesh element $\varOmega _i$, the flux variable $Q_i$ is written in the terms of element-average and face-average pressures:

$$\begin{aligned} {\varvec{Q}}_i=\varvec{\varLambda }{\varvec{W}}_i({\varvec{e}}_i X_i-{\varvec{\pi }}_{i}) \, \end{aligned}$$

(6)

where ${\varvec{Q}}_i=[Q_1,Q_2,\ldots ,Q_m]^{\rm T}$ is the vector of fluxes on the interface, m is the number of faces of element $\varOmega _i$, $\varvec{\varLambda }$ is the matrix of relative permeabilities, $e_i=(1,\ldots ,1)^{\rm T}_{1\times m}$, $X_i$ is the scalar quantity at the centroid of element $\varOmega _i$, ${\varvec{\pi }}_i$ is the vector of scalar quantities defined at the centroid of faces of element $\varOmega _i$. ${\varvec{W}}_i$ is a positive definite matrix, which is a key part of MFD method. The quantities for MFD are shown in Fig. 2. A linear pressure can be written in the form $p=X_i a+b$, where a and b are constant vector and scalar constant scalar quantity. Here, we consider the capillary pressure and gravity, thus their quantities at the centroid and face are satisfied by the linear relationships below:

$$\begin{aligned} p_i-\pi _{k,i}= & ({\varvec{x}}_i -{\varvec{x}}_k)\cdot {\varvec{a}}_p, \end{aligned}$$

(7)

$$\begin{aligned} p_{c,i}-\pi _{c,k,i}= & ({\varvec{x}}_i -{\varvec{x}}_k)\cdot {\varvec{a}}_{\rm c}, \end{aligned}$$

(8)

$$\begin{aligned} z_i-\pi _{z,k,i}= & ({\varvec{x}}_i -{\varvec{x}}_k)\cdot {\varvec{a}}_z, \end{aligned}$$

(9)

where $p_i$ is the pressure at the centroid of element $\varOmega _i$, $\pi _{k,i}$ is the interface pressure at interface $A_k$; $p_{c,i}$ is the capillary pressure at the centroid of element $\varOmega _i$, $\pi _{c,k,i}$ is the interface capillary pressure at interface $A_k$; $z_i$ is the depth at the centroid of element $\varOmega _i$, $\pi _{z,k,i}$ is the interface depth at interface $A_k$; ${\varvec{x}}_i$, ${\varvec{x}}_k$ are the coordinate vector of the centroid of element $\varOmega _i$ and the centroid of interface $A_k$, respectively.

According to Darcy’s law (Eq. 2), the flux through the interface $A_k$ of element $\varOmega _i$ is

$$\begin{aligned} Q_{i,k}=\int _{A_k}{\varvec{v}} {\rm d}s=-|A_k|{\varvec{v}}_k \cdot \hat{{\varvec{n}}}_k, \end{aligned}$$

(10)

where ${\varvec{v}}$ is the velocity, ${\varvec{v}}_k$ is the Darcy velocity on the interface $A_k$, $|A_k|$ is the area of interface $A_k$, $\hat{{\varvec{n}}}_k$ is the outward unit normal vector on interface $A_k$.

Based on the linear relationships (Eqs. 7–9) and Eq. 2, the above equation can be rewritten as:

$$\begin{aligned} Q_{\alpha ,i,k}=-|A_k|\lambda k \cdot ({\varvec{a}}_{\alpha ,p}-{\varvec{a}}_{\alpha ,c}-\rho _\alpha g {\varvec{a}}_z) \end{aligned}$$

(11)

where $\alpha$ is oil or water phase.

Substituting Eq. 11 and Eqs. 7–9 into Eq.6, we can see that the matrix $W_i$ satisfies the following conditions:

$$\begin{aligned} {\varvec{W}}_i {\varvec{X}}={\varvec{N}} {\varvec{K}}, \end{aligned}$$

(12)

where ${\varvec{X}}=[{\varvec{x}}_1-{\varvec{x}}_i,\ldots ,{\varvec{x}}_k-{\varvec{x}}_i,\ldots ,{\varvec{x}}_m-{\varvec{x}}_i ]^{\rm T}$, ${\varvec{N}}=[ |A_1|\hat{{\varvec{n}}}_1,\ldots , |A_k|\hat{{\varvec{n}}}_k,\ldots , |A_m|\hat{{\varvec{n}}}_m ]^{\rm T}$. m is the number of interfaces within element $\varOmega _i$.

The way to obtain a symmetric positive-definite matrix $W_i$ has been discussed by Lie et al. (2012). Here, we use the following inner product for $W_i$:

$$\begin{aligned} {\varvec{W}}_i=\frac{1}{|\varOmega _i|}{\varvec{N}} {\varvec{K}} {\varvec{N}}^{\rm T}+ \frac{6}{d}tr({\varvec{K}}) {\varvec{A}}_i ({\varvec{I}}_i-\bar{{\varvec{Q}}_i}\bar{ {\varvec{Q}}_i}^{\rm T}) {\varvec{A}}_i. \end{aligned}$$

(13)

where d denotes the spatial dimension. ${\varvec{A}}_i$ is a diagonal matrix composed by interface areas. $\bar{{\varvec{Q}}_i}=orth({\varvec{A}}_i {\varvec{X}}_i)$ is orthonormal basis, where ${\varvec{X}}_i$ is a diagonal matrix composed by the distance from centroids to interface centroids. $tr({\varvec{K}})$ is the trace of ${\varvec{K}}$. ${\varvec{I}}_i$ is an identity matrix.

After calculating the matrix ${\varvec{W}}_i$, the interface flux $Q_{\alpha ,i}$ of phase $\alpha$ can be rewritten as

$$\begin{aligned} {\varvec{Q}}_{\alpha ,i}=\varvec{\varLambda }_{\alpha }{\varvec{W}}_i[{\varvec{e}}_i ({\varvec{p}}_{\alpha }-{\varvec{p}}_{\alpha ,c}-\rho _\alpha g \nabla {\varvec{z}})- ({\varvec{\pi }}_{\alpha }-\pi _{\alpha ,c}-\rho _\alpha g \nabla \varvec{\zeta }) ], \end{aligned}$$

(14)

where ${\varvec{p}}_\alpha$, ${\varvec{p}}_{\alpha ,c}$ and ${\varvec{z}}$ are element pressure vector, capillary pressure vector and depth vector for element $\varOmega _i$, respectively. ${\varvec{\pi }}_\alpha$, ${\varvec{\pi }}_{\alpha ,c}$ and $\varvec{\zeta }$ are interface pressure vector, capillary pressure vector and depth vector for interfaces of element $\varOmega _i$, respectively.

We re-arrange and simplify the interface flux $Q_{\alpha ,i}$ of phase $\alpha$, and can get the numerical expression of eq. 14:

$$\begin{aligned} Q_{\alpha ,i}=\frac{k_{r,\alpha }}{\mu _\alpha }\left[ (L_{W_i} p_{\rm o} - p_{c,\alpha }-\rho _\alpha gD) -\sum ^{N_{\rm f}}_{j=1}(W_{ij}\pi _{o,j}-\pi _{c,\alpha ,j}-\rho _{\alpha }g\pi _{D,j})\right] \qquad \alpha = o,w \end{aligned}$$

(15)

where $L_{W_i}=\sum _{j=1}^{N_{\rm f}}W_{i,j}$, $p_{\rm o}$ is the oil pressure, $p_{c,\alpha }$ is the capillary pressure of phase $\alpha$, D is the depth, $\pi _{o,j},\ \pi _{c,\alpha ,j},\ \pi _{D,j}$ are the Lagrange multipliers on interface j for oil pressure, capillary pressure of phase $\alpha$ and depth, respectively.

The standard approach in modeling the flow within Darcy’s equation focuses on treating the permeability as a scalar value. Although such a treatment results in simplified equations that be solved with less computation times compared to full tensor permeability implementation, the latter is needed naturally fractured-reservoir modeling to correctly solve fluid flow problems in a variety of realistic settings. The full tensor permeability ${\mathbf {K}}$ is then written as,

$$\begin{aligned} {\mathbf {K}}= \left[\begin{array}{lll} k_{xx} & k_{xy} & k_{xz} \\ k_{yx} & k_{yy} & k_{yz} \\ k_{zx} & k_{zy} & k_{zz} \\ \end{array}\right] \end{aligned}$$

(16)

The off-diagonal permeability elements $k_{xy}, k_{xz}, k_{yx}, k_{yz}, k_{zx}, k_{zy}$ account for the dependence of flow rate in one direction on pressure differences in orthogonal directions while the diagonal permeability elements $k_{xx}, k_{yy}, k_{zz}$ represent the dependence of flow rate in one direction on pressure differences in the same direction. Note that the full tensor permeability ${\mathbf {K}}$ is embedded in the formulation of the ${\mathbf {W}}$ matrix Zhang and Abushaikha (2019).

Combining with Eq. (15), the velocity in Eq. 1 can be expressed in the terms of flux, . We can then represent the governing equations (Eq. 1) for oil and water in space and time as:

$$\begin{aligned}&V_e \frac{\partial \phi (1-S_{\rm w}) \rho _{\rm o}}{\partial t} -\sum _{i=1}^{N_{\rm f}}Q_{o,i}=\rho _{\rm o} q_{\rm o} \end{aligned}$$

(17)

$$\begin{aligned}&V_e \frac{\partial \phi S_{w} \rho _{\rm w}}{\partial t}- \sum _{f=1}^{N_{\rm f}}Q_{w,i}=\rho _{\rm w} q_{\rm w} \end{aligned}$$

(18)

$V_e$ is the volume of the element $\varOmega _i$. $Q_{\alpha ,i}$ is the face flux of the element $\varOmega _i$. The rock is assumed to be slightly compressible, and thus approximated by:

$$\begin{aligned} \phi \approx \phi _{\rm o}(1+c_r(p-p_{sc})) \end{aligned}$$

(19)

Assuming that the fluxes are continuous at a givens shared interface by two elements, the governing equations are then discretized in time using Euler approximation, and written as a residual term:

$$\begin{aligned} R_{o}= & \frac{ V_e }{\Delta t } \left[ \phi (p_{\rm o})^{n+1}(1-S_{w,e}^{n+1})-\phi (p_{\rm o})^{n}(1-S_{w,e}^{n})\right] -\lambda ^{*}_{o}(S_{\rm w})^{n+1}\sum _{i=1}^{N_{\rm f}}Q_{o,i} -V_e\rho _{\rm o} q_{\rm o} \end{aligned}$$

(20)

$$\begin{aligned} R_{w}= & \frac{ V_e }{\Delta t } \left[ \phi (p_{\rm o})^{n+1}(S_{w,e}^{n+1})-\phi (p_{\rm o})^{n}(S_{w,e}^{n})\right] -\lambda ^{*}_{w}(S_{\rm w})^{n+1} \sum _{i=1}^{N_{\rm f}}Q_{w,i} -V_e\rho _{\rm w} q_{\rm w} \end{aligned}$$

(21)

where e represents the element number, $\Delta t$ is the time step, $\lambda _{\alpha }^*$ is the upstream mobility of the phase and indices n and n + 1 refers to the previous and current newton iteration, respectively.

In the MFD method, we need to ensure flux continuity at the interface by adding a constraint represented by the momentum balance equations. If the interface of a certain element is not shared by another element (i.e., interface is at boundary), then the flux at that interface is assumed to be zero:

$$\begin{aligned} Q_{\alpha ,i}=\lambda _{\rm o}^{*}(S_{\rm w})^{n+1}\left[ (L_{W_i} p_{\rm o} - p_{c,\alpha }-\rho _\alpha gD) -\sum ^{N_{\rm f}}_{j=1}(W_{ij}\pi _{o,j}-\pi _{c,\alpha ,j}-\rho _{\alpha }g\pi _{D,j})\right] =0 \end{aligned}$$

(22)

However, if an interface is shared by two adjacent element, then $Q_{\alpha ,i}=-Q_{\alpha ,i}^{'}$. This implies that the Lagrange multipliers enforced at the interface are equal as well. The expanded total flux $Q_{t_{\alpha ,i}}$ at the shared interface for each phase $\alpha$ can be written as:

$$\begin{aligned} Q_{t_{\alpha ,i}}=Q_{\alpha ,i}+Q_{\alpha ,i}^{'} \end{aligned}$$

(23)

The details of the derivation and the final equations for approximating solutions of the original system have been clarified in Zhang and Abushaikha (2019).

Fig. 2

Schematic of matrix and fracture elements. a Grid analysis used to define the mimetic inner product. b Discretization of the matrix and fractures

2.2 Fracture Discretization Using DFDM

For the discrete fracture treatment, the above flux equations are only applicable for fractures in a lower space-dimension. Thus, the interface flux $Q^{\rm F}_{\alpha ,i}$ of phase $\alpha$ for fracture system can be rewritten as

$$\begin{aligned} Q^{\rm F}_{\alpha ,i}= & \lambda ^{\rm F}_{\rm w} \left[ (L^{\rm F}_W (p^{\rm F}_{o,i}-\rho _\alpha g \nabla D^{\rm F}_i)-\sum ^{N^{\rm F}_{\rm f}}_{j=1}W^{\rm F}_{ij}(\pi ^{\rm F}_{i,j}-\rho _\alpha g \nabla \pi ^{\rm F}_{D,i,j})\right] \\&+\lambda ^{\rm F}_{\rm w} \left( L^{\rm F}_W p^{\rm F}_{c,i}-\sum ^{N^{\rm F}_{\rm f}}_{j=1}W^{\rm F}_{ij}\pi ^{\rm F}_{c,i,j}\right) \end{aligned}$$

(24)

where $L^{\rm F}_W=\sum _{j=1}^{N^{\rm F}_{\rm f}} W^{\rm F}_{i,j}$ and $N^{\rm F}_{\rm f}$ is the number of edges of fractures. $p^{\rm F}_{c,i}$ is the capillary pressure of phase $\alpha$ for fracture system. $D^{\rm F}_i$ is the depth, $\pi ^{\rm F}_{o,j},\ \pi ^{\rm F}_{c,i,j},\ \pi ^{\rm F}_{D,i,j}$ are the Lagrange multipliers on interface j for oil pressure, capillary pressure of phase $\alpha$ and depth for fracture system, respectively. The superscript F represents fracture system.

For the fracture–matrix connections, a fracture F is consider as a lower-dimensional object as interior face, and the fracture element usually represented as a face of a matrix face, as shown in Fig. 3. Obviously, the fracture element pressure $p^{\rm F}_{\rm o}$ are part of the matrix interface pressure $\pi _i$ ($\pi _j$ ). When calculating the final linear system, the fracture interface pressures are eliminated and only the matrix element pressures are kept.

Fig. 3

a Relative permeabilities for the strongly water-wet SWW, weakly water-wet WWW and mixed-wet MW cases, as well as b the respective capillary pressures based on Blunt (2017). The spontaneous imbibition process stops at the respective $Sw^*$

Simultaneously, the flux in Eq. 15 is written locally for all faces within each matrix element. In order to couple the matrix and fracture elements together, the following continuity conditions for the flux are imposed at each interface of two neighboring elements (taking the interface F of $\varOmega _i$ and $\varOmega _j$ as an example in Fig. 3:

• (1) If F is a matrix element, the flux satisfies:

  $$\begin{aligned} Q_{i}+Q_{j}=0, \end{aligned}$$

  (25)

• (2) If F is a fracture, the fracture flux satisfies:

  $$\begin{aligned} \sum _i Q^{\rm F}_i= Q_{i,k}+ Q_{j,k}+q^{\rm F}, \end{aligned}$$

  (26)

where $q^{\rm F}$ is the sink/source for fracture. The above equation is the mass conservation relation in the fracture.

Thus, the equations for matrix and fracture systems can be coupled through the matrix–fracture transfer function (Eq. 25). The details of the derivation for this part can be referred to Zhang and Abushaikha (2019).

3 Building the Numerical Model

The rock properties for the matrix were developed using the parameters presented in Schmid et al. (2016) study. The relative permeability curves were constructed based on the power law model:

$$\begin{aligned} k_{\rm rw} &= k_{\rm rw max}\left[ \frac{S_{\rm w} - S_{\rm wi}}{1-S_{\rm wi} - S_{\rm or}}\right] ^n \end{aligned}$$

(27)

$$\begin{aligned} k_{\rm ro} &= k_{\rm rwo max}\left[ \frac{1-S_{w} - S_{\rm or}}{1-S_{\rm wi}-S_{\rm or}}\right] ^m \end{aligned}$$

(28)

where $k_{{\rm rwo max}}$ is the maximum relative permeability of oil, $k_{\rm rw max}$ is the maximum relative permeability of water, n and m are the relative permeability exponents, $S_{\rm w}$ is the water saturation, $S_{\rm wi}$ is the initial water saturation, $S_{\rm or}$ is the residual oil saturation, $k_{\rm rw}$ is water relative permeability and $k_{\rm ro}$ is the oil relative permeability.

On the other hand, the capillary pressure prediction model for generalized mixed-wet systems Skjaeveland et al. (1998) is written as:

$$\begin{aligned} P_{\rm c}=\frac{c_{\rm w}}{\left( \frac{S_{\rm w}-S_{\rm wr}}{1-S_{\rm wr}}\right) ^{a_{\rm w}}} + \frac{c_{\rm o}}{\left( \frac{S_{\rm o}-S_{\rm or}}{1-S_{\rm or}}\right) ^{a_{\rm o}}} \end{aligned}$$

(29)

where $P_{\rm c}$ is the capillary pressure, a and c are constants representing either drainage or imbibition processes. All constants are positive except for $c_{\rm o}$, and the capillary curve consists of two asymptotic branches: positive water branch and negative oil branch. A summary of the test parameters used for different wetting cases is found in Table 1, while the relative permeability and capillary curves are plotted in Fig. 3

Table 1 Rock and fluid properties for rocks with different wetting states

4 Numerical Convergence: Model Verification

We aim to demonstrate the accuracy and the converge of the fully implicit MFD method with discrete fractures for different wetting states in an oil-water system when spontaneous and forced imbibition are considered. The ability of the method to accurately predict [rgb]0,0,1thermodynamic equilibrium between phases and map the changes in the water and oil saturation profiles are presented next.

We consider a 3D domain with a homogeneous matrix of size 20 m $\times$ 20 m $\times$ 100 m. Initially, the domain is saturated with equal volumes of oil and water; the water phase lies on top of the oil phase. The domain is fully sealed and the fluids do not interact with the surrounding thus eliminating any external influence on the flow mode. Since gravity forces dictate the flow in the domain along with the capillary forces of the oil and water phases, the water tend to flow vertically downwards due to the imbalance between the capillary pressure and gravity until the equilibrium is reached. This behavior is governed by the equation below that predicts the height of the water column as it moves:

$$\begin{aligned} H=\frac{P_{\rm c}(S_{\rm w})}{(\rho _{\rm w}-\rho _{\rm o})g}, \end{aligned}$$

(30)

In order to study the behavior of spontaneous imbibition flow in the domain with capillary and gravity forces for the MFD method, three different cases have been considered to account for varying fracture locations as shown in Fig. 4. Each model is tested for three different wetting states described in Table 1 while Table 2 shows the fluids and matrix properties.

Fig. 4

3D reservoir model for the validation test: a model with no fracture. (The mesh is composed by 1043 tetrahedral elements.) b model with one fracture (aperture = 0.1 cm) (The mesh is composed by 1191 tetrahedral elements, where 40 triangles are for the discrete fracture model); c model with three fractures (aperture = 0.1 cm). (The mesh is composed by 1643 tetrahedral elements, where 122 triangles are for the discrete fracture model.)

Table 2 Test 1 fluid and other parameters

The comparison between the numerical and the analytical solutions is demonstrated in Fig. 5 while Fig. 6 shows the water saturation profiles at equilibrium for the Models I, II and III generated through the simulation of the MFD method in a 3D domain. Based on Eq. 30, the computed height of the water column from the numerical solutions agrees well with the analytical solution of the capillary curves. The water imbibes downwards because of the imbalance between capillary forces and gravity forces. After those forces reach equilibrium, the water saturation remains unchanged. The good agreements show that this approach models the capillary pressure and gravity cases correctly for different wetting states of the rock and different fracture placements.

Fig. 5

The analytical and numerical solution both plotted for Models I, II and III. The solutions coincide for different wetting states regardless of the existence and the location of fractures in the domain, proving an effective estimation of the water column height using the MFD method. The solution for all models are the same since only gravity and capillary forces are allowed to influence the flow of fluids. The wetting state will only change the rate at which a solution will be achieved, but does not affect the final equilibrium conditions

Fig. 6

The oil and water saturations at the end of the simulation are practically similar for the three different models regardless of the changes in the wetting state and fracture existence. Initially, water is positioned on top of the oil; as the simulation starts, fluid exchange takes place and the water displaces the oil upwards due to gravity segregation and capillary forces. The location of the fractures in the system, and the rock wetting states dictates the time taken to reach equilibrium. At the end of the simulation, the oil is fully mobilized to the upper part of the domain and sits on top of the water as the visualizations show

5 Applied Computational Cases

The main objective of the tests is to model the behavior of the flow in naturally fractured reservoirs for different wetting states. The capillary forces have a significant impact on the amounts of oil recovered from a domain and the pressure distribution profiles as will be shown. We designed the tests to compare and highlight the effect of various capillary curves on the recovery factor, and we further added multiple full permeability tensor cases with different rotation angles and compared it to the scalar approach to show the ability of the MFD scheme used in incorporating the effects of the full tensor. The recovery factors were discussed in the light of the wetting state and the direction of the permeability tensor.

In this section, we present two tests to demonstrate the capability of the MFD method in robustly and accurately predicting the flow of oil and water in complex domains with the existence of multiple fractures, and for different capillary pressure profiles. First, we test the method for a fractured domain where the fractures are aligned along the injector and the producer, as we observe closely the behavior of the flow while varying the capillary pressure and the permeability tensor. The second test comprises of two full reservoirs connected with a communicating fracture, with a producer and an injector placed in each reservoir.

We analyze the simulation results in terms of the pore volume injected (PVI) written as:

$$\begin{aligned} {\rm PVI}=\frac{1}{V_{\rm p}}\int _0^t q_{\rm t} {\text {d}}t \end{aligned}$$

(31)

where $V_{\rm p}$ is the total pore volume of the model, $q_{\rm t}$ is the injected rate. We also measure Courant–Friedrichs–Lewy (CFL) number as a condition for convergence while solving for the governing equations (Fig. 7).

Fig. 7

Test I: the reservoir model in a has the dimensions of 10.0 $\times$ 10.0 $\times$ 1.0 m) with three fractures present (aperture = 1.0 cm). The mesh in b is an unstructured grid composed of 1583 tetrahedral elements with 3636 interfaces, where 98 triangles are for the discrete fractures

5.1 Test I: Two Phase Flow with 3 Fractures

In the first test, we examine the effect of permeability tensor on the spontaneous imbibition process while varying the wetness of the matrix in the domain. The matrix permeability tensor is represented as:

$$\begin{aligned} K_m=R_{\rm o}(-\theta )\text {diag}(30.0,1,1)R^{\rm T}_{\rm o}(-\theta ), \end{aligned}$$

(32)

where $R_{\rm o}(-\theta )$ is the rotation matrix. The permeability tensor of fracture is $K_{\rm f}$ = diag(500,500,500) and the aperture is 0.1 cm. We place two wells on the opposite sides of the domain, while aligning the flow from the injector to the producer along the existing fractures. The movement of the fluids in the well is controlled by a constant flow rate of 0.02 m$^3$/day). The other parameters are listed in Table 3.

Table 3 Fluid and other parameters for Test 2

Here, we consider three cases; a case where the permeability is a scalar and fixed at value of 30 md (C1), and two other cases represent by the following rotation angles:

$$\begin{aligned} \theta = {\left\{ \begin{array}{ll} 0^{\circ } &\text{ for } o=x \quad \text{ Case } \text{2 } \text{(C2) } \\ 45^{\circ } & \text{ for } o=z \quad \text{ Case } \text{3 } \text{(C3) } \end{array}\right. } \end{aligned}$$

(33)

The domain is initially fully saturated with oil and is swept through slow continuous water injection. The water saturation profiles are illustrated in Fig. 8 for three different permeability scenarios. We notice that when capillary forces are ignored, the flow of water from the injector to the producer is mainly controlled by the positioning of the faults and their relative placement to the major principal direction of the permeability tensor. The water tends to flow in the least resistant path (i.e., faults with higher permeability than the matrix). When the angle of the tensor is aligned parallel to the faults and the direction of flow (45$^{\circ }$), the water flows into the major fault sweeping the surrounding oil. We notice that the water breakthrough at the producer is fast compared to the other permeability cases where the flow occurs in the major x-direction of the domain. The pressure and velocity profiles in Figs. 9 and  10, respectively, confirm the earlier observation as the rate of flow along the major fracture when a rotation angle of 45$^{\circ }$ is used is almost double compared to the scalar and 0$^{\circ }$ cases.

Fig. 8

Test I: water saturation after 0.25 PVI for different permeability cases. First row: C1, middle row: C2, last row: C3. Each column represents no Pc, SWW, WWW and MW cases, respectively

The general observations on the behavior of the flow as the permeability tensor changes holds true even when capillarity is introduced and the rock wetness is varied. If we look at the second column of visualizations in Fig. 8, we notice that the sweep efficiency of water is much higher when capillarity is considered. The water front moves uniformly across the domain and mobilizes the oil. The direction of the water front is dictated by the angle of the permeability tensor, and the water does not have a concentrated velocity in the fracture conduits as compared to the case when capillary effects are absent. This behavior is attributed to the strongly water wet property of the rock, where water tend to be imbibed by the throat pores of the domain, and thus expelling the oil out towards the producer. The strong capillary forces associated with the wetness of the rock towards water overcomes the high permeability of the fracture channel, thus causing the water to move uniformly in the matrix and recover more oil. This is confirmed by the corresponding pressure and velocity profiles in Figs. 9 and 10, respectively, as the velocity of the flow in mainly uniform across the domain.

Fig. 9

Test I: Pressure profiles after 0.25 PVI for different permeability cases. First row: C1, middle row: C2, last row: C3. Each column represents no Pc, SWW, WWW and MW cases, respectively

Fig. 10

Test I: Velocity profiles after 0.25 PVI for different permeability cases. First row: C1, middle row: C2, last row: C3. Each column represents no Pc, SWW, WWW and MW cases, respectively

Subsequently, the same analysis can follow for when WWW and MW wettability cases are encountered. The WWW behavior is closely similar to that of the SWW case, but we tend to see some overshooting in water saturation through the fracture conduits and at the producer. This is explained by the decrease in the capillary forces due to the wetting nature of the rock, thus allowing the water to escape through the fractures. This observation is more evident in the MW case where the rock have regions favoring both oil and water equally. This causes the flow of water in the matrix to be much slower compared to the SWW and WWW cases, but still faster than when no capillarity is introduced. In the MW case, the water flows rapidly in the fractures channels and an early breakthrough happens at the prouder well. This is most evident when the rotation angle for the permeability tensor is 45$^{\circ }$. The pressure and velocity profiles dictate that the oil is drained slowly while the water travels fast across the domain.

The curves plotted in Fig. 11a, b show the water cut and the recovery factors, respectively. The hydrocarbon recovery factor is computed using:

$$\begin{aligned} RF=\frac{\text {Produced Hydrocarbons}}{\text {Initial Hydrocarbons in Place}} \end{aligned}$$

(34)

The water cut profiles agrees with the water saturation profiles analyzed earlier. When capillarity is ignored and a tensor permeability with 45$^\circ$ rotation angle is used, we notice an early breakthrough of water at 18 days which increases slowly and steadily till the end of the simulation. However, when the permeability is treated as a 0$^\circ$ tensor or scalar, water breakthrough is delayed till 90 days but happens instantly as the water at the producer jumps from 0 to 80% in 1 day and then keeps increasing gradually. This indicates, that once the water approaches the fracture channels, it flows all the way to the producer as no resistance from the rock matrix is encountered. The early breakthrough means that oil is not efficiently displaced yielding only a 50% recovery factor for all permeability cases.

Fig. 11

a Water cut at the producer for different cases versus the simulation time in days b Recovery factor versus the simulation time for different stability cases and varying permeability tensor. The water breakthrough is the fastest when capillarity is ignored, but a significant amount of hydrocarbon is left unswept. Recovered oil is higher with capillarity especially when the rock is strongly water wet

On the other hand, introducing capillarity does not affect the final water content at the producer but rather controls the rate at which the water breakthrough is achieved. We notice from Fig. 11a that water arrives at the producer first when MW conditions are used, followed by WWW and SWW subsequently. In SWW and WWW cases, and once the water front approaches the fracture channels, the water cut increases rapidly in the producer from 0 to 40%, followed by a gradual increase till 100%. The SWW case yields the highest recovery factor at around 78% as the water tends to imbibe the pores of the rock and expel the oil out. In brief, the final values of the water cut are not greatly affected by the changes in the wettability of the rock, but rather the path of the curve and the time at water breakthrough happens. On the other hand, the wetting state affects remarkably the amount of the oil produced from a certain model.

All the permeability scenarios had the same CFL number and the same number of time steps and newton iterations regardless whether capillarity was considered or not as shown in Fig. 12a–d. This indicates that the variations in the permeability tensor does not affect the computational performance of the simulator but only changes the physics of the problem under study. We also notice that the newton iterations and the CPU time needed to arrive to a plausible solution almost doubled when capillarity is considered. The same observation holds true regardless of the wetting sate of the rock. The capillary forces poses physical complexities on the problem especially with the existence of fractures, and thus more computational effort is needed to converge to a physical solution.

Fig. 12

Test I computational results. All the scenarios have the same average CFL when capillarity is considered, but the CPU time triples compared to wen Pc = 0

Capillarity forces greatly affect the kinetics of the recovery in the hydrocarbons reservoirs as we explained here. In this test, we explored the differences in the solution generated by different capillary pressure curves, and demonstrated the robustness of the MFD methods in predicting the flow and recovery factors for varying permeability tensor in the presence of faults. In the next test, we impose a more challenging domain with two reservoirs connected by a fault to further test the accuracy of the MFD method with capillarity present.

5.2 Test II: Two Reservoirs Connected with a Geological Fault

In this test, we model two reservoirs connected by one fracture, as shown in Fig. 13a. The fault has the dimensions of 933.0 $\times$ 500.0 $\times$ 1.7 m and is treated as a 2D surface. The resultant mesh is shown in Fig. 13b. The model consists of 8960 hexahedral elements for the matrix and 256 quadrilateral elements for the fracture. Initially, both upper and lower compartments of the reservoir are saturated with oil, and the other physical properties of the fluids are listed in Table 4. An injector well is placed in the bottom reservoir operating at a constant rate of 1000 m$^3$/day. Also, a producer flowing at the same rate is positioned in the upper compartment of the reservoir. The direction of the flow is mapped by the relative permeability and the capillary pressure curves for which a weakly water rock is chosen and follows the models presented in Table 1. The capillary pressure is assumed to be zero in the fracture while the relative permeability profiles are linear.

Fig. 13

Test II: The domain a is composed of two reservoirs and one connecting fractures while the mesh b is composed of 8960 hexahedral elements representing the matrix and 256 quadrilateral elements for the fracture

Table 4 Test 3 fluid and other parameters

The matrix permeability tensor is fixed $K_m=\text {diag}(60,60,60)$ while the permeability of the fracture is written as:

$$\begin{aligned} K_{\rm f}=R_{\rm o}(-\theta )\text {diag}(60,60,10000)R^{\rm T}_{\rm o}(-\theta ), \end{aligned}$$

(35)

where $R_{\rm o}(-\theta )$ is the rotation matrix. Here, we consider three cases; a case where the permeability is a scalar and fixed at value of 1000 md (C1), and two other cases represent by the following rotation angles:

$$\begin{aligned} \theta = {\left\{ \begin{array}{ll} 0^{\circ } &\text{ for } o=x \quad \text{ Case } \text{2 } \text{(C2) } \\ 30.96^{\circ } & \text{ for } o=y \quad \text{ Case } \text{3 } \text{(C3) } \end{array}\right. } \end{aligned}$$

(36)

In Figs. 14, 15 and 16, we show the water saturation, pressure and velocity profiles at 700 and 1200 days of simulation time. Let us first examine C1 and C2 after 700 days of continuous water injection. In these cases, we notice that the oil starts to mobilize slowly in the bottom reservoir without reaching the connecting fault when capillarity is not considered. The vertical direction of the flow is perpendicular to the alignment of the fault, causing resistance to the flow. This means that the velocity of the flow will be low as water travels though the matrix domain. On the contrary, accounting for the WWW capillary curves will permit the water to interact with the matrix and thus get imbibed into the pores of the rock. This capillary force dictates a faster movement of the water, thus reaching the fracture earlier. However, if we look at C3, we see clearly that the water has not reached the fault after 700 days regardless whether capillarity is considered or not. This is explained by the alignment of the fault permeability tensor with the direction of the flow at an angle of 30.96$^\circ$. This means that water will travel slower due to fluid-matrix interactions and the permeability tensor alignment; an observation that is evident when the velocity profiles are observed in Fig. 16a. After 1200 days of continuous injection, water reaches to the top of the reservoir for C1 and C2 when both $P_{\rm c}$ is zero and capillarity is considered. The main difference is that the WWW nature of the rock implies higher water saturation at the top reservoir, and faster flow of fluids in the fracture due to the direction of the permeability tensor.

Fig. 14

Test II: Water saturation profiles after a 700 days and b 1200 days. The top row refers to the cases when capillarity is ignored whole the bottom row represents the WWW condition. Also, each column refers to cases C1, C2 and C3 for varying fault permeability

Fig. 15

Test II: Pressure profiles after a 700 days and b 1200 days. The top row refers to the cases when capillarity is ignored whole the bottom row represents the WWW condition. Also, each column refers to cases C1, C2 and C3 for varying fault permeability

Fig. 16

Test II: Velocity profiles after a 700 days and b 1200 days. The top row refers to the cases when capillarity is ignored whole the bottom row represents the WWW condition. Also, each column refers to cases C1, C2 and C3 for varying fault permeability

In Fig. 17, the water cut and the recovery factors are plotted. The observed trends of the curves are quite similar to those examined in Test I, as early breakthrough is noticed for all permeability variations when capillarity is not considered. However, the distance between the water cut profiles is not large, because the vertical upwards direction of the flow assisted by the imbibition process into the matrix for capillarity cases counteracts the downward gravity forces in the domain. The capillary forces assist the flow in traveling through the fracture and all the way to the producer. The recovery profiles for the oil produced show the average recovery factor is almost 5% higher when capillarity is introduced compared to when Pc = 0. This indicates that the fluid-matrix interactions support a better drainage of the hydrocarbons in place. One interesting observation that is different from Test I is that the water breakthrough is earlier when the permeability is treated as a scalar. This is attributed to the angles and the rotation axis considered when the permeability is represented as a tensor in the fault, spanning both the x and y planes. This means that the direction of the full permeability tensor will not be perfectly aligned with the direction of the flow in the fracture and hence will create a resistance to the flow. This necessitates that the water will mover slower in the fault plane and will take longer time to reach to the producer, and hence a delayed breakthrough is expected. This opposes the results of the previous test as the major principal direction of the permeability tensor in the matrix was parallel to the fracture placement, creating a highly conductive channel for water to reach the producer faster when permeability tensors are modeled.

Fig. 17

a Water cut at the producer for different cases versus the simulation time in days and b recovery factor versus the simulation time for different stability cases and varying permeability tensor. The water breakthrough is the fastest when capillarity is considered, but a significant amount of hydrocarbon is recovered. The velocity profiles at the wells and the fault are shown in a slice section of the domain in c

In order to examine the movement of fluid in the fracture plane more closely, we take a sliced cross section in the y-direction of the reservoir-fault system. The velocity profiles in the fracture and the two wells are plotted in Fig. 17c for the different capillarity and permeability scenarios. The velocity of the flow is computed from the information on the fluxes passing through the system and the resolution of the grids. The velocity is observed to be much higher in the fracture compared to that in the well, because a finer grid is used to represent the fracture elements. We notice that for all permeability cases, the velocity is higher when capillary pressure is used; an observation consistent with the conclusion drawn from the saturation and water cut profiles earlier. Furthermore, the angle of the permeability tensor did not affect the velocity of the flow in the fault as the profiles of cases C2 and C3 overlap. However, for case C1, the velocity is the highest when scalar permeability is used, as the flow occurs naturally without any opposition.

Finally, we show some computational results in Fig. 18 which highlights that the introduction of capillarity into the problem add physical complexities to the modeling process and thus requiring a larger number of steps and nonlinear iterations to converge to a desirable solution. However, the final answer would be more representative of the real nature of fluid flow in the porous media.

Fig. 18

Test II computational results. The nonlinear iterations increase in b because of the high complexity of the problem due to the imbibition phenomena in b

6 Conclusions and Recommendations

The results of the numerical study presented in this work show the robustness and the efficiency of the novel fully implicit numerical approach combining MFD and DFM techniques to model two phase flow in fractured porous media in the presence of capillarity effects. The imbibition phenomenon was investigated during the production of oil from fractured reservoirs, and was successfully analyzed in the light of the wettability of the rock while considering full permeability tensors and unstructured grids.

The following conclusions can be drawn from this work:

• (1) The comparison between the numerical and the analytical solutions for spontaneous imbibition processes in a closed system with different wetting rock states showed a great match as the oil and water phases mobilize to achieve equilibrium in the system.

• (2) The orientation of the fracture and the full permeability tensor affects greatly the speed of water breakthrough in certain producing systems. If the principal direction of the permeability tensor is aligned with the fracture orientation, water tends to flow in the fracture first thus missing oil spots in the matrix (depends greatly on the wetness of the rock.

• (3) The amounts of recovered water at the producer are not greatly affected by the changes in the wettability of the rock, but rather the path of the curve and the time at water breakthrough happens. On the other hand, the wetting state affects remarkably the amount of the oil produced from a certain model.

• (4) Introducing capillary forces between the fractures and the matrix impose physical complexities on the problem, which in turn adds to the computational time and the newton iterations needed by the MFD method to converge to a satisfactory solution.

In summary, this work introduced a new approach to study the effect of the wettability of the rock and the associated capillary forces on the performance of the fully implicit MFD method to predict oil recovery ratios with fractures present in the system. The MFD method proved to be robust in preserving the physics of the problem, and accurately mapping the flow path in the investigated domains. Future work will focus on comparing the performance of the MFD method against other discretization schemes for the same application of modeling spontaneous imbibition, and investigating the effect of heterogeneous permeability on the wettability of the rock and the oil recovery process.

Appendix

According to the no-flow boundary conditions and continuous conditions for matrix and fracture systems, we can get the residual terms as follows:

$$\begin{aligned} R_{u_\alpha 1}= & \biggl [ \lambda _{\alpha ,i}[(L_{T_i} (p_{o,i}-\rho _\alpha g \nabla D_i)-\sum ^{N_{\rm f}}_{l=1}T_{il}(\pi _{i,l}-\rho _\alpha g \nabla \pi _{D,i,l})] \biggr ] ^{n+1} \end{aligned}$$

(37)

$$\begin{aligned} R_{u_\alpha 2}= & \sum _\alpha \left[ \lambda _{\alpha ,i}\left[ (L_{T_i} (p_{o,i}-\rho _\alpha g \nabla D_i)-\sum ^{N_{\rm f}}_{l=1}T_{il}(\pi _{i,l}-\rho _\alpha g \nabla \pi _{D,i,l})\right] \right. \\&\left. + \lambda _{\alpha ,j}\left[ (L_{T_j} (p_{o,j}-\rho _\alpha g \nabla D_j)-\sum ^{N_{\rm f}}_{l=1}T_{jl}(\pi _{j,l}-\rho _\alpha g \nabla \pi _{D,j,l})\right] \right] ^{n+1} \\&+ \sum _\alpha \left[ \lambda ^{\rm F}_{\alpha ,i}\left[ (L^{\rm F}_{T_i} (p^{\rm F}_{o,i}-\rho _\alpha g \nabla D^{\rm F}_i)-\sum ^{N_{\rm f}}_{l=1}T^{\rm F}_{il}(\pi ^{\rm F}_{i,l}-\rho _\alpha g \nabla D^{\rm F}_{i,l})\right] \right. \\&\left. +\lambda ^{\rm F}_{\alpha ,j}\left[ \left( L^{\rm F}_{T_j} (p_{o,j}-\rho _\alpha g \nabla D^{\rm F}_j)-\sum ^{N_{\rm f}}_{l=1}T_{jl}(\pi ^{\rm F}_{j,l}-\rho _\alpha g \nabla \pi ^{\rm F}_{D,j,l})\right] \right] ^{n+1}\right. \end{aligned}$$

(38)

$$\begin{aligned} R_{u_{\rm c}}= & \left[ \left( L_{T_i} p_{c,i}-\sum ^{N_{\rm f}}_{l=1}T_{ij}\pi _{c,i,l}\right) +\left( L_{T_j} p_{c,j}-\sum ^{N_{\rm f}}_{l=1}T_{ij}\pi _{c,j,l}\right) \right] ^{n+1} \end{aligned}$$

(39)

$$\begin{aligned} R^{\rm F}_{u_\alpha 1}= & \left[ \lambda ^{\rm F}_{\alpha ,i}\left[ \left( L^{\rm F}_{T_i} (p^{\rm F}_{o,i}-\rho _\alpha g \nabla D^{\rm F}_i)-\sum ^{N_{\rm f}}_{l=1}T^{\rm F}_{il}(\pi ^{\rm F}_{i,l}-\rho _\alpha g \nabla \pi ^{\rm F}_{D,i,l})\right] \right] ^{n+1}\right. \end{aligned}$$

(40)

$$\begin{aligned} R^{\rm F}_{u_{\rm c}}= & \left[ \left( L^{\rm F}_{T_i} p_{c,i}-\sum ^{N_{\rm f}}_{l=1}T^{\rm F}_{ij}\pi ^{\rm F}_{c,i,l}\right) +\left( L^{\rm F}_{T_j} p_{c,j}-\sum ^{N_{\rm f}}_{l=1}T^{\rm F}_{ij}\pi ^{\rm F}_{c,j,l}\right) \right] ^{n+1} \end{aligned}$$

(41)

$$\begin{aligned} R_{S_\alpha }= & |\varOmega _i |\frac{( \phi S_\alpha )^{n+1}-(\phi S_\alpha )^n}{\bigtriangleup t} \\&+ \sum ^{N_{\rm f}}_{i=1}\left[ \lambda _\alpha \left[ (L_T (p_{o,i}-\rho _\alpha g \nabla D_i)-\sum ^{N_{\rm f}}_{j=1}T_{ij}(\pi _{i,j}-\rho _\alpha g \nabla \phi _{D,i,j})\right] \right] ^{n+1}-|\varOmega _i |(q_\alpha )^{n+1} \end{aligned}$$

(42)

$$\begin{aligned} R_{S_{\alpha }^{\rm F}}= & a_F|A_F |\frac{(\phi ^{\rm F} S^{\rm F}_\alpha )^{n+1}-(\phi ^{\rm F} S^{\rm F}_\alpha )^{n} }{\bigtriangleup t} \\&+ \sum ^{N_F}_{i=1}\left[ \lambda _\alpha \left[ (L_T^{\rm F} (p^{\rm F}_{o,i}-\rho _\alpha g \nabla z^{\rm F}_i)-\sum ^{N_F}_{j=1}T^{\rm F}_{ij}(\pi ^{\rm F}_{i,j}-\rho _\alpha g \nabla \phi ^{\rm F}_{D,i,j})\right] \right] ^{n+1}-a_F|A_F|(q^{\rm F}_\alpha )^{n+1} \end{aligned}$$

(43)

Consider the set of primary variables $y=\{p_{\rm o},\pi _{\rm o},\pi _{\rm c},\pi ^{\rm F}_{\rm o},\pi ^{\rm F}_{\rm c},S_{\rm w},S^{\rm F}_{\rm w}\}$, the linear system on each nonlinear iteration has the following form:

$$\begin{aligned} \left[\begin{array}{lllllll} \frac{\partial R_{u_\alpha 1} }{\partial p_{\rm o}} & \frac{\partial R_{u_\alpha 1} }{\partial \pi _{\rm o}} & \frac{\partial R_{u_\alpha 1} }{\partial \pi _{\rm c}} & \frac{\partial R_{u_\alpha 1} }{\partial \pi ^F_{\rm o}} & \frac{\partial R_{u_\alpha 1} }{\partial \pi ^F_{\rm c}} & \frac{\partial R_{u_\alpha 1} }{\partial S_{\rm w}} & \frac{\partial R_{u_\alpha 1} }{\partial S^F_{\rm w}} \\ \frac{\partial R_{u_\alpha 2} }{\partial p_{\rm o}} & \frac{\partial R_{u_\alpha 2} }{\partial \pi _{\rm o}} & \frac{\partial R_{u_\alpha 2} }{\partial \pi _{\rm c}} & \frac{\partial R_{u_\alpha 2} }{\partial \pi ^F_{\rm o}} & \frac{\partial R_{u_\alpha 2} }{\partial \pi ^F_{\rm c}} & \frac{\partial R_{u_\alpha 2} }{\partial S_{\rm w}} & \frac{\partial R_{u_\alpha 2} }{\partial S^F_{\rm w}} \\ \frac{\partial R_{u_{\rm c}} }{\partial p_{\rm o}} & \frac{\partial R_{u_{\rm c}} }{\partial \pi _{\rm o}} & \frac{\partial R_{u_{\rm c}} }{\partial \pi _{\rm c}} & \frac{\partial R_{u_{\rm c}} }{\partial \pi ^F_{\rm o}} & \frac{\partial R_{u_{\rm c}} }{\partial \pi ^F_{\rm c}} & \frac{\partial R_{u_{\rm c}} }{\partial S_{\rm w}} & \frac{\partial R_{u_{\rm c}} }{\partial S^F_{\rm w}} \\ \frac{\partial R^F_{u_\alpha 1} }{\partial p_{\rm o}} & \frac{\partial R^F_{u_\alpha 1} }{\partial \pi _{\rm o}} & \frac{\partial R^F_{u_\alpha 1} }{\partial \pi _{\rm c}} & \frac{\partial R^F_{u_\alpha 1} }{\partial \pi ^F_{\rm o}} & \frac{\partial R^F_{u_\alpha 1} }{\partial \pi ^F_{\rm c}} & \frac{\partial R^F_{u_\alpha 1} }{\partial S_{\rm w}} & \frac{\partial R^F_{u_\alpha 1} }{\partial S^F_{\rm w}} \\ \frac{\partial R^F_{u_{\rm c}} }{\partial p_{\rm o}} & \frac{\partial R^F_{u_{\rm c}}}{\partial \pi _{\rm o}} & \frac{\partial R^F_{u_{\rm c}}}{\partial \pi _{\rm c}} & \frac{\partial R^F_{u_{\rm c}}}{\partial \pi ^F_{\rm o}} & \frac{\partial R^F_{u_{\rm c}}}{\partial \pi ^F_{\rm c}} & \frac{\partial R^F_{u_{\rm c}}}{\partial S_{\rm w}} & \frac{\partial R^F_{u_{\rm c}}}{\partial S^F_{\rm w}} \\ \frac{\partial R_{S_\alpha }}{\partial p_{\rm o}} & \frac{\partial R_{S_\alpha }}{\partial \pi _{\rm o}} & \frac{\partial R_{S_\alpha }}{\partial \pi _{\rm c}} & \frac{\partial R_{S_\alpha }}{\partial \pi ^F_{\rm o}} & \frac{\partial R_{S_\alpha }}{\partial \pi ^F_{\rm c}} & \frac{\partial R_{S_\alpha }}{\partial S_{\rm w}} & \frac{\partial R_{S_\alpha }}{\partial S^F_{\rm w}} \\ \frac{\partial R_{S^F_\alpha }}{\partial p_{\rm o}} & \frac{\partial R_{S^F_\alpha }}{\partial \pi _{\rm o}} & \frac{\partial R_{S^F_\alpha }}{\partial \pi _{\rm c}} & \frac{\partial R_{S^F_\alpha }}{\partial \pi ^F_{\rm o}} & \frac{\partial R_{S^F_\alpha }}{\partial \pi ^F_{\rm c}} & \frac{\partial R_{S^F_\alpha }}{\partial S_{\rm w}} & \frac{\partial R_{S^F_\alpha }}{\partial S^F_{\rm w}} \\ \end{array}\right] \left[\begin{array}{l} \delta p_{\rm o} \\ \delta \pi _{\rm o} \\ \delta \pi _{\rm c} \\ \delta \pi ^F_{\rm o}\\ \delta \pi ^F_{\rm c}\\ \delta S_{\rm w}\\ \delta S^F_{\rm w} \\ \end{array}\right] = \left[\begin{array}{l} R_{ p_{\rm o} } \\ R_{ \pi _{\rm o} }\\ R_{ \pi _{\rm c} }\\ R_{ \pi ^F_{\rm o}}\\ R_{ \pi ^F_{\rm c}}\\ R_{ S_{\rm w}}\\ R_{ S^F_{\rm w} } \\ \end{array}\right] \end{aligned}$$

(44)

where $\frac{\partial R_x}{\partial y}$ is Jacobian matrix of unknowns y, $R_y$ is the residual for the mass balance and saturation equations. As a result, we simultaneously solve all the governing equations fully implicit in time.