Document Type: Research Paper
Authors
^{1} Assistant Professor, Department of Chemical Engineering, Jask Branch, Islamic Azad University, Jask, Iran
^{2} Professor, Department of Chemical and Petroleum Engineering, Sharif University of Technology, Tehran, Iran
^{3} Associate Professor, Faculty of Chemical, Petroleum, and Gas Engineering, Semnan University, Semnan, Iran
Abstract
Keywords
Main Subjects
1. Introduction
The hydrocarbon recovery from underground reservoirs highly depends on the porous media structure. Some books are available which discuss the variety of features of the porous medium (Pop and Ingham, 2001; Ingham and Pop, 2005; Vafai, 2005; Nield, Bejan et al., 2006; Vadász, 2008; Vafai, 2010). The significant portion of petroleum reserves in the world is hosted by fractured reservoirs that are considered to have a composition of two distinctive media: the first medium of fractures has high permeability and a low volume of storage, and the second medium is less permeable and has a high volume storage of rock matrix (Warren and Root, 1963). Imbibition in such reservoirs, which entail the displacement of the nonwetting phase (NWP) with the wetting phase (WP), is a significant mechanism for recovery of gas/oil. Gas/oil phases are considered as NWP during water flooding, aquifer invasion, and in any general enhanced methods of recovery; on the other hand, the WP is a fluid which is aqueous. When the saturated matrix blocks of NWP are contacted, the aqueous fluid can imbibe into the matrix blocks from the enclosed fractures (Cai, Yu et al., 2010; Cai and Yu, 2011; Cai, Hu et al., 2012).
Spontaneous imbibition can either be cocurrent or countercurrent varying on the conditions of the boundary. Countercurrent spontaneous imbibition (COUCSI) occurs when the permeable faces of all matrix blocks and a wetting phase come into contact. Cocurrent spontaneous imbibition (COCSI) only occurs when a permeable surface portion contacts with the wetting phase while the nonwetting phase covers the rest of permeable surfaces. The process type is mainly controlled by some factors, including fracturing degree, gravitational forces strength, the relative position of matrix block with respect to the wetting phase source (i.e. aquifer or injection wells), and the wetting phase injection rate into the network of fracture. COCSI is considered where there are dominant capillary forces. This indicates that there is no change in the matrix surface enclosed by the wetting and nonwetting phases in the experiment setup during the process. Also, there are other studies on the recovery performance of a matrix blocks stack or matrix block under the proceeding level of the wetting phase in the fracture network (Schechter, Zhou et al., 1994; PooladiDarvish and Firoozabadi, 1998; PooladiDarvish and Firoozabadi, 2000; Qasem, Nashawi et al., 2008; Standnes, 2010).
Significant tools to study physical processes are analytical models. Handy obtained a onedimensional watergas COCSI model based on the pistonlike displacement assumption which showed the linear variance of gas recovery with the times square root (Handy, 1960). Li and Horne employed an approximate model in onedimensional COCSI analysis considering gravity (Li and Horne, 2006). The basis of their model was on the restricting assumption of the pressure profile of linear capillary. A unique type of their estimated model, when gravity has no importance, was later offered by Li (Li, 2007). Bourbiaux (Bourbiaux, 2009) suggested an easy approximate model of the transfer of matrixfracture wateroil. Schmid et al. (Schmid, Geiger et al., 2011) realized that the McWhorter and Sunada (McWhorter and Sunada, 1990) proposed onedimensional analytical solution to the twophase immiscible flow applied to the problem of onedimensional COCSI without simplifying assumptions. SaboorianJooybari and Khademi proposed a travelling wave solution meant for COCSI issue; The suggested solution was, on the other hand, limited to specific capillary pressure and relative permeabilities functional forms (SaboorianJooybari and Khademi, 2014). The finding of Schmid et al. (Schmid, Geiger et al., 2011) was used by MirzaeiPaiaman and Mashi (Mirzaei‐Paiaman and Masihi, 2014) to apply universal scaling equations to onedimensional COCSI. A simple scaling equation that is free of relative permeabilities and capillary pressure characteristics was extracted using universal equations, applied to some routine applications where there is unavailability of these data. A simple theory was presented by Haugen et al. (Haugen, Fernø et al., 2014) based on the assumption of pistonlike displacement. Simple capillary tube models have been used to model COCSI in porous media (Cai, Yu et al., 2010; Cai, Perfect et al., 2014; MirzaeiPaiaman, 2015). These models generally do not deem the full characteristics of real porous media. The literature review shows no precise analytical solution to COCSI general case where displacement might happen in over onedimension. Even the exact analytical solution to the COCSI problem of onedimension only applies only to pure COCSI without considering the nonwetting phase backflow production at the face of inlet (McWhorter and Sunada, 1990; Schmid and Geiger, 2012). Analytical and semianalytical expressions are extremely significant in transferring laboratory tests performance in field behavior prediction. The major challenge of previous studies is their lack of considering gravity effect (Kashchiev and Firoozabadi, 2003; Mason, Fischer et al., 2010; Standnes, 2010; Standnes, 2010; Li, 2011; Ruth and Arthur, 2011).
Similar to COCSI in fractured reservoirs, a lot of attention has been paid to COUCSI (PooladiDarvish and Firoozabadi, 1998; Firoozabadi, 2000; PooladiDarvish and Firoozabadi, 2000; Morrow and Mason, 2001; Karimaie, Torsaeter et al., 2006; Qasem, Nashawi et al., 2008; Karpyn, Halleck et al., 2009; Schmid, Geiger et al., 2011; Mason and Morrow, 2013; MirzaeiPaiaman, Masihi et al., 2013; Schmid and Geiger, 2013). Iffly et al. (Iffly, Rousselet et al., 1972) reported wateroil COCSI and COUCSI experiments and numerical simulations results. The COCSI tests, as shown by their experiments, were faster and gave a higher recovery of oil than their equivalent COUCSI pairs. As shown by their numerical simulations, the relative permeabilities required by COCSI tests to match experimental recovery were higher compared to relative permeabilities required by COUCSI tests. Bourbiaux and Kalaydjian (Bourbiaux and Kalaydjian, 1990) conducted COCSI and COUCSI studies on waterwet sandstones in a linear system. There was a faster recovery of oil and slightly higher final recovery of oil in the COCSI experiments compared to COUCSI experiments. The COCSI displacement relative permeabilities as revealed by numerical simulations were around 30% higher than the COUCSI relative permeabilities. Such behaviors were attributed to the phenomenon of viscous coupling (Bentsen and Manai, 1993). In their COCSI experiment, 1.2% of the original oil was countercurrently produced in water contact from the bottom face. There was a cocurrent production of 37.2% of the original oil that was in place from the top face. Kantzas et al. (Kantzas, Pow et al., 1997) performed COCSI and COUCSI experiments on carbonate formation and Berea sandstone in Alberta. In the COCSI test, one plug’s face was in contact with water as air freely came out from other surfaces. Kantzas et al. compared the displacement rate and ultimate recoveries of COCSI and COUSI. The analysis of numerical simultaneous performance on one and twodimensional COCSI and COUCSI were also done by PooladiDarvish and Firoozabadi (PooladiDarvish and Firoozabadi, 2000). As to their observations, COCSI oil recovery was more efficient compared to COUCSI as per the rate of displacement and final recovery. The small contribution of the face in water contact to the oil recovery was also noticed in the COCSI. Standness (Standnes, 2004) used COCSI and COUCSI to study the effect of different boundary conditions on the displacements in chalk samples. As per the observations, COCSI oil recovery was faster and efficient as per final recovery compared to COUCSI. Also, it was qualitatively observed that almost the entire oil was cocurrently produced in the COCSI tests. The surface production of countercurrent backflow in water contact constituted an approximate estimation of 510% of the entire production. The oil was quickly produced after the imbibition test initiation. Tests of onedimensional waterair COCSI and COUCSI were performed by Hatiboglu and Babadagli (Hatiboglu and Babadagli, 2010) at two temperatures on Berea sandstone samples. The displacement rate at a low temperature and COCSI final recovery was observed to be faster and higher compared to COUCSI. Nevertheless, even though COCSI ultimate recovery was higher at a higher temperature, the displacement of COUCSI was faster compared to the COCSI test. They claimed to have conducted wateroil COCSI tests too. However, both of the used boundary conditions were not according to the boundary conditions of common pure COCSI. In a single case, saturated porous medium of oil bottom face was in contact with water and was open to atmosphere via the top face. In a different case, both the top face and bottom surface of the onedimensional medium were in contact with water. As discussed previously, dominant COCSI happens when barely a permeable surface and wetting phase are in contacts, and the rest of the permeable surfaces are enclosed by a nonwetting phase. Bentsen and Trivedi (Haugen, Fernø et al., 2014) stated that the observed reduced flux in the countercurrent flow in comparison to cocurrent flow, i.e. a reduction in driving force for every unit volume, maybe because of the capillary coupling. The conducted literature review reveals that while there are similarities between COCSI and COUCSI, the two processes dynamically differ in the residual saturation and displacement rate of the nonwetting phase. In the majority of research during the past century which has addressed spontaneous imbibition, the local equilibrium has been assumed to be instantaneously reached. Conversely, several studies have implicated the existence of nonequilibrium effects in the course of spontaneous imbibition indicating that this phenomenon is basically an unsteadystate procedure (Barentblatt, 1971; Barenblatt and Vinnichenko, 1980; Barenblatt and Gilman, 1987; Barenblatt, Entov et al., 1989; Barenblatt, Patzek et al., 2003; Silin and Patzek, 2004; Le Guen and Kovscek, 2006; Juanes, 2008; Jafari, Masihi et al., 2017; Jafari, Masihi et al., 2018; Jafari, Masihi et al., 2017). The main feature of an unsteadystate process is that capillary pressure and relative permeability vary as a function of time.
In this study, considering the imbibition process as well as the water and oil movement, as the wet phase and nonwet phase, in a block saturated by oil and surrounded by two vertical fractures full of water, a differential equation having partial and nonlinear derivatives is introduced using Darcy and mass balance equations. On the other hand, as there is no analytical solution for this equation, a new equation is introduced by considering a wide range of viscosity values for the wet and nonwet phases and by selecting the best suitable functions for relative permeability and capillary pressure. Considering the boundary conditions governing the countercurrent imbibition, an analytical solution (equation) is developed. Finally, the developed equation is validated. The results of this research can be very important for a better understanding of the imbibition process and the water and oil movement in the fractured environments.
2. Imbibition in a single block surrounded by two vertical fractures filled with water
A schematic view of the considered single block surrounded by two vertical fractures filled with water is presented in Figure 1; a homogeneous matrix saturated with oil (NW), the onedimensional displacement from the matrix block is understood as the left and right sides open to flow, while the lateral faces are impermeable. It should be mentioned that the imbibition process is done through the block.
Figure 1
A schematic view of imbibition process in a block surrounded by vertical fractures.
The velocities of water and oil phases in each specific length can be calculated by considering the Darcy law and mass balance equations as follows:
(1) 

(2) 
In this process, the total velocity can be considered as the summation of water and oil velocities:
(3) 
The combination of the abovementioned equations will lead to the following equation:
(4) 
As we know, the relationships between water and oil potentials and capillary pressure can be defined as follows:
(5) 

(6) 

(7) 

(8) 
If the block is considered in a horizontal position, the gravitational force can be ignored and Equation 8 is summarized as below:
(9) 

(10) 
Therefore, by rearranging the velocity equation for the nonwet phase, the following equation can be obtained:
(11) 
On the other hand, the mass balance equation for the nonwet phase in porous media can be written as follows:
(12) 
By assuming incompressibility for the rock and fluid, Equation 12 can be summarized into the following equation:
(13) 
The combination of Equations 11 and 13 result the following equation:
(14) 
Considering that both oil and gas phases are assumed as incompressible fluid, the flow rate in each section of the rock is equal. Therefore, the derivative of the total velocity term with respect to the length of core is equal to zero and Equation 14 can be summarized as follows:
(15) 
Consequently, the differential equation of the fluid flow in porous media for imbibition process in a single block surrounded by two vertical fractures will be given by:
(16) 

(17) 
Equation 17 is a nonlinear partial differential equation and cannot be solved through analytical methods. In the following, by considering different intervals for the viscosity of the wet and nonwet phases as a new assumption, we can summarize the equation to obtain a new linear upscaling equation.
3. Considering the near values of the water and oil viscosities ( )
By considering the near values of water and oil viscosities, one can use an average value for the viscosity of the system. There are several methods to estimate the average value of the viscosity of a system consisting of different fluids. Many of these methods recommend the volume averaging technique. Gambill presents the following relation to calculate the average value for the viscosity of a system containing oil and water (Gambill, 1959; Maples, 2000).
(18) 
Therefore, the average viscosity can be given by:
(19) 
Hence, the differential equation can be expressed by:
(20) 
By defining the following dimensionless parameters:
(21) 
The form of Equation 20 changes as follows:
(22) 
Thus, the dimensionless time for imbibition process in a core surrounded by two vertical fractures can be defined as follows:
(23) 
Consequently, the final dimensionless form of the differential equation can be obtained by:
(24) 
4. Determining saturation functions for the near values oil and water viscosities
In order to solve Equation 24, determining proper functions for saturation, relative permeability parameters, and capillary pressure is needed. Relative permeability can be presented as an exponential function of saturation as follows:
(25) 
Saturation exponent (n) is calculated empirically through laboratory data matching. Many analytical studies consider a constant value for exponent n to simplify the differential equation (Equation 24). For example, in many studies, linear relative permeability functions reasonably calculate the saturation distribution (Pavone, 1989; Firoozabadi and Ishimoto, 1994; Firoozabadi, Ishimoto et al., 1994). Generally, capillary pressure shows a complex behavior, and there are several functions to describe it (Brooks and Corey, 1964; Bentsen and Anli, 1976)
(26) 

(27) 
Introducing the capillary pressure and relative permeability functions (Equations 26 and 27) to the governing differential equation for imbibition process (Equation 24) results in the following equation:
(28) 
The simplifications of Equation 28 leads to Equation 29:
(29) 
5. Assuming a water viscosity higher than the oil viscosity ( )
Assuming a water viscosity higher than the oil viscosity means that:
(30) 
Therefore, governing differential equation can be written as:
(31) 
Applying the abovementioned assumption (Equation 30) changes Equation 31 to the following form:

(32) 
The remarkable point is that all the expressions related to the viscosity and relative permeability of the nonwet phase are eliminated from the equation. Because the wet phase viscosity is much larger than the notwet phase viscosity, the injection flow of the wet phase is a limiting factor and is the main flow.
6. Determining saturation functions a water viscosity higher than the oil viscosity
By comparing Equations 32 and 24, one can find that Equation 32 has a simpler solution than Equation 24. Therefore, assuming more realistic functions of relative permeability and capillary pressure is easier for Equation 32 than Equation 24. As mentioned above, the goal is to transform the differential equation into a linear and solvable differential equation. To this end, it is needed to express capillary pressure function for different saturation exponents (n) in such a way that the differential equation is linearized and becomes compatible with the physics of the actual capillary pressure function. The result of plotting capillary pressure function at different values of n is tabulated in Table 1.
Table 1
Water relative permeability and capillary pressure values for different saturation exponents (n)
P_{c} 
K_{rw} 
n 
1 

2 

4 
According to Figure 2, despite the difference in the form of capillary pressure formulas presented in Table 1, the predictions show logical trends.
Figure 2
Calculated capillary pressures at different powers (n) of relative permeability equation.
If the equation is solved for the general form of capillary pressure ( ), the below result is acquired:
(33) 
The general form of capillary pressure is defined as follow:
(34) 
The final result of Equation 34 can be presented as follows (WC, 1967):
(35) 
The general hyper geometric function is also defined by [61]:
(36) 
where, a and b are two vectors with a length of p; (a_{1})_{k} and (b_{1})_{k} are Puchhamer signs define as follows:
(37) 
where, n is an integer; therefore:
(38) 
Based on the abovementioned equations, the total form of capillary pressure equation can be presented by:
(39) 
Simplifying Equation 39 changes the form of this equation to:
(40) 
By applying the derivative of capillary pressure and relative permeability to the differential equation, the simplified equation will be given by:
(41) 
The dimensionless form of Equation 41 is as follows:
(42) 
Comparing Equations 29 and 42 showed that by considering two different assumptions for viscosities, the same differential equations were derived. By considering different boundary conditions, this differential equation must be solved analytically for the countercurrent imbibition process, and the new scaling equation derived should be validated.
7. Analytical solution for countercurrent imbibition process in a single block surrounded by two vertical fractures filled with a wet phase
Figure 3 displays a block saturated with oil and surrounded by two vertical fractures filled with a wet phase. At each boundary, the water phase enters the block, and the oil phase exits from the block countercurrently.
Figure 3
A schematic view of the countercurrent imbibition process in a core plug.
In order to investigate the countercurrent imbibition process, the following boundary conditions are considered:
(43) 
The saturation function is also defined by:
(44) 
By introducing these variables to the abovementioned equation, the following equation is obtained:
(45) 
By solving this equation through variable separation method, the following equation is achieved:
(46) 
The new boundary conditions for this equation are given by:
(47) 
Applying the first boundary condition will results in the following equation:
(48) 
We can now apply the second boundary condition as follows:
(49) 
Also, we have:
(50) 

(51) 
The Fourier function corresponding to the saturation function is expressed by:
(52) 
Consequently, the general solution can be written as follows:
(53) 
We can find the coefficients by using these two initial conditions:
(54) 

(55) 
Therefore, C_{n} coefficient must be found in a way that the following Taylor series equals 1.
(56) 
Consequently, the final form of the saturation equation for this model will be expressed by:
(57) 
8. Results and discussion
As shown in the abovementioned sections, using Darcy law and mass balance equations for imbibition process in a block saturated with oil as the nonwet phase and surrounded by two vertical fractures filled with water as the wet phase leads to a nonlinear partial differential equation. This equation cannot be solved by analytical methods; hence, by considering different intervals for the viscosities of the wet and nonwet phases and by finding the best functions for relative permeability and capillary pressure, a new equation was introduced and solved by applying the governing boundary conditions of countercurrent imbibition process. As can be seen in Figure 3, the boundary conditions of the countercurrent imbibition process in a block saturated with oil and surrounded by two vertical fractures filled with water are in such a way that block is initially saturated with oil; then, over the time, water enters the core, and oil is produced around the block. Therefore, oil saturation decreases in the block. Figure 4 illustrates the results of the analytical model for the prediction of saturation. As shown in Figure 4, initially, the oil saturation is equal to one, and then the oil saturation decreases until the dimensionless oil saturation approaches irreducible oil saturation at infinite time. The remarkable point that can be seen in Figure 4 is that the oil production from the block occurs symmetrically. Because similar properties are assumed for the vertical fracture, at any time, the maximum saturation occurs in the center of the block.
Figure 4
Saturation prediction results by the analytical method.
Figure 5 shows the normal saturation prediction results of the countercurrent imbibition process in a block surrounded by two vertical fractures. The result presented in Figure 5 is compatible to the governing conditions, which shows the reliability of the obtained scalingup equation.
Figure 5
The normal saturation prediction result of the countercurrent imbibition process in a block surrounded by two vertical fractures.
The summary of the modeling results of the countercurrent imbibition process in a block surrounded by two vertical fractures filled with water are also listed in Table 2.
Table 2
Modeling results of the countercurrent imbibition process in a block surrounded by two vertical fractures filled with water.


K_{r} 

P_{cD} 

t_{D} 

PDE 

New scalingup equation 
9. Conclusions
Based on the main subject of this study, a new partial differential equation, which is the governing equation of fluid flow in porous media, was achieved by employing the Darcy’s law and material balance. As this equation is highly nonlinear, presenting an analytical solution is not possible. Thus, we considered a few assumptions for the gravity values of the fluids and saturation functions so as to be able to solve it. Afterwards, the partial differential equation was solved in the case of countercurrent imbibition for a matrix block which is surrounded by two vertical fractures filled with the wetting phase. The validation was performed using MATLAB software showing that the proposed equation is accurate enough.
Nomenclature
COCSI 
Cocurrent spontaneous imbibition 
COUCSI 
Countercurrent spontaneous imbibition 
NWP 
Nonwetting phase 
WP 
Wetting phase 