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

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.


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 1 3 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 cocurrent flow. As the name imply, counter-current is when the wetting phase and the nonwetting 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 backpressure. 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, countercurrent imbibition is predominant and requires sensitive care when modeled. The work of  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. 2011Schmid et al. , 2016Nooruddin 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).
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 Fig. 1 Different imbibition scenarios. a and b display geometries where there is a combination of counterand co-current spontaneous imbibition. c and d show counter-current spontaneous imbibition. Reproduced from  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 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 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.

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: where is the porosity of the medium, S denotes the saturation of the phase, is the density of the phase, is Darcy's velocity, and q is the mass flow rate. Darcy's velocity can be written for each phase as: where k r is the phase relative permeability, is the phase viscosity, is the intrinsic permeability tensor, g is the vertical pressure gradient and ∇ 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.
The phase mobility is expressed by:

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 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 can be estimated using = 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 i , the flux variable Q i is written in the terms of element-average and face-average pressures: where Q i = [Q 1 , Q 2 , … , Q m ] T is the vector of fluxes on the interface, m is the number of faces of element i , is the matrix of relative permeabilities, e i = (1, … , 1) T 1×m , X i is the scalar quantity at the centroid of element i , i is the vector of scalar quantities defined at the centroid of faces of element i . 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: where p i is the pressure at the centroid of element i , k,i is the interface pressure at interface A k ; p c,i is the capillary pressure at the centroid of element i , c,k,i is the interface capillary pressure at interface A k ; z i is the depth at the centroid of element i , z,k,i is the interface depth at interface A k ; x i , x k are the coordinate vector of the centroid of element i and the centroid of interface A k , respectively.
According to Darcy's law (Eq. 2), the flux through the interface Based on the linear relationships (Eqs. 7-9) and Eq. 2, the above equation can be rewritten as: where 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: 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 : is orthonormal basis, where X i is a diagonal matrix composed by the distance from centroids to interface centroids. tr(K) is the trace of K . I i is an identity matrix.
After calculating the matrix W i , the interface flux Q ,i of phase can be rewritten as where p , p ,c and z are element pressure vector, capillary pressure vector and depth vector for element i , respectively. , ,c and are interface pressure vector, capillary pressure vector and depth vector for interfaces of element i , respectively.
We re-arrange and simplify the interface flux Q ,i of phase , and can get the numerical expression of eq. 14: where is the capillary pressure of phase , D is the depth, o,j , c, ,j , D,j are the Lagrange multipliers on interface j for oil pressure, capillary pressure of phase 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 is then written as, 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 is embedded in the formulation of the matrix .
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: V e is the volume of the element i . Q ,i is the face flux of the element i . The rock is assumed to be slightly compressible, and thus approximated by: 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: where e represents the element number, Δt is the time step, * 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: However, if an interface is shared by two adjacent element, then Q ,i = −Q � ,i . This implies that the Lagrange multipliers enforced at the interface are equal as well. The expanded total flux Q t ,i at the shared interface for each phase can be written as: The details of the derivation and the final equations for approximating solutions of the original system have been clarified in Zhang and Abushaikha (2019).

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 F ,i of phase for fracture system can be rewritten as i,j are the Lagrange multipliers on interface j for oil pressure, capillary pressure of phase 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 F o are part of the matrix interface . When calculating the final linear system, the fracture interface pressures are eliminated and only the matrix element pressures are kept. 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 i and j as an example in Fig. 3: (1) If F is a matrix element, the flux satisfies: (2) If F is a fracture, the fracture flux satisfies: where q 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 .

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: . 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 * where k rwo max is the maximum relative permeability of oil, k rw max is the maximum relative permeability of water, n and m are the relative permeability exponents, S w is the water saturation, S wi is the initial water saturation, S or is the residual oil saturation, k rw is water relative permeability and k 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: where P c is the capillary pressure, a and c are constants representing either drainage or imbibition processes. All constants are positive except for c 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

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 × 20 m × 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: 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.
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.

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 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 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: 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 where V p is the total pore volume of the model, q 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).

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: where R o (− ) is the rotation matrix. The permeability tensor of fracture is K 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.
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:  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 • ), 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 • is used is almost double compared to the scalar and 0 • cases.
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.
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, 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 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 • . 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: Produced Hydrocarbons Initial Hydrocarbons in Place (a) (b) 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 The water cut profiles agrees with the water saturation profiles analyzed earlier. When capillarity is ignored and a tensor permeability with 45 • 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 • 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.
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.
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.

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 × 500.0 × 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.
The matrix permeability tensor is fixed K m = diag(60, 60, 60) while the permeability of the fracture is written as: where R o (− ) 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: 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 • . 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 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.
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 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 conductive channel for water to reach the producer faster when permeability tensors are modeled.
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.

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 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.

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 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 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.

(a) (b)
(c) 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 (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.
Consider the set of primary variables y = {p o , o , c , F o , F c , S w , S F w } , the linear system on each nonlinear iteration has the following form: where R x 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. (44)