1 Introduction

Two-phase flow in fractured rocks is relevant to various application, ranging from geothermal energy (Babaei and Nick 2019; Erfani et al. 2019), groundwater contamination (Gleeson et al. 2009; Sebben et al. 2015; Aminnaji et al. 2019; Aljuboori et al. 2019; Dejam 2019), carbon capturing and storage (Emami-Meybodi et al. 2015; Babaei and Islam 2018; Erfani Gahrooei and Joonaki 2018; Dejam and Hassanzadeh 2018; Erfani et al. 2020; An et al. 2020b; Ershadnia et al. 2020), and soil salinization (Weisbrod and Dragila 2006) to enhanced hydrocarbon recovery (Farajzadeh et al. 2012; Rokhforouz and Akhlaghi Amiri 2017; Rokhforouz and Amiri 2018; Bakhshi et al. 2018; Rabbani et al. 2020a; Chen et al. 2020; An et al. 2020a). The complexity of the fluid flow in highly heterogeneous rocks is caused by the notion that the fluid flows in two different media, matrix block and fracture, which separates the fluids into flowing and stagnant regions due to property contrast (Erfani et al. 2019; Hasan et al. 2019; Soltanian et al. 2019; An et al. 2020c; Aziz et al. 2020; Rabbani et al. 2020b; Hasan et al. 2020). In recent decades, increase in worldwide oil demand has motivated many studies in this field, since naturally fractured reservoirs (NFR) hold about 60% of the total oil reserves in the world (Joonaki et al. 2016; Hakimov et al. 2019). Unique characteristics and heterogeneous nature of such reservoirs, in addition to different oil production mechanisms, make them relatively complicated to model and understand compared to conventional reservoirs.

Gravity drainage and spontaneous imbibition are important production mechanisms in NFRs (Saidi et al. 1979; Saidi 1983; Austad and Standnes 2003; Ameri et al. 2015; Aljuboori et al. 2019). Gravity drainage is a process in which the gravity acts as the main driving mechanism and the non-wetting phase (i.e., gas phase) pushes the oil out of the rock. It is an effective production mechanism in the gas-invaded zone due to a larger density difference. Spontaneous imbibition is usually effective in the water-invaded zone, where water is imbibed into the rock due to capillary forces and oil is produced to the fractures. Gravity drainage phenomenon can be classified into two main types: free-fall and forced gravity drainage. Free-fall gravity drainage is usually the major oil production mechanism in high permeability and thick reservoirs with a great drawdown pressure (Hagoort 1980; Chatzis and Ayatollahi 1993), while forced gravity drainage is a result of gas injection. Gravity and capillary are significant driving forces in gravity drainage. Along with wettability, the efficiency of the process is controlled by block-to-block interactions such as capillary continuity, which preserves the hydraulic continuity of the oil phase over fractures, and reinfiltration, occurring when the produced oil from a block enters the underneath matrix block due to gravity and/capillary forces (Dullien 2012). When the length of a liquid element is larger than the fracture aperture, the upper face of the lower porous matrix block will touch the liquid element, and consequently, a liquid bridge will be formed between the two blocks (Fig. 4). The existence of liquid bridges and reinfiltration provides the hydraulic continuity between the blocks and decreases the threshold height of the upper block. On the one hand, the bulk flow regime in the matrix and vertical fractures and, on the other hand, film flow in the horizontal fracture strongly influence the two-phase gas–oil flow in fractured oil reservoirs (Rostami et al. 2010).

Numerous experimental studies have been conducted on the gravity drainage mechanisms. The majority of these studies were performed on synthetic porous media such as coarse quartz grains sand packs with high permeability, so that the duration of each experiment, as well as the set-up height, can be adequately decreased (Moon and Bruce 2007; Zendehboudi et al. 2011; Saedi et al. 2015; Teng et al. 2016). Saidi et al. (1979) studied the interaction of matrix blocks and the presence of capillary continuity and reinfiltration between adjacent matrix blocks. Horizontal fracture aperture of 10 \(\mu\)m was shown to be necessary to form stable liquid bridges in the horizontal fractures. Major results reported by researchers in this field could be classified into two general categories: first, increasing the model height and decreasing capillary pressure improve the process efficiency; second, the presence of a driving viscous force leads to earlier breakthrough and can improve the ultimate recovery factor (i.e., forced gravity drainage). Notably, lower recovery factor by increasing the gas injection rate is also reported after a critical gas injection rate (Rostami et al. 2010; Akhlaghi et al. 2012; Saedi et al. 2015; Zobeidi et al. 2016).

Firoozabadi and Markeset (1992) and Firoozabadi and Markeset (1994) used spacers of different thicknesses (up to 1.2 mm) as horizontal fractures in a multi-stack block system. They showed the existence of capillary continuity in the fractures with the aperture of up to 1.2 mm and concluded that liquid bridges play an important role in controlling the gravity drainage mechanism. Effects of horizontal fracture properties on reinfiltration process have been examined by Sajadian et al. (1998), documenting that the reinfiltration rate significantly depends on the oil supply from the upper blocks as well as the spacer specifications. Zobeidi and Fassihi (2018) showed the importance of block-to-block interactions in miscible/immiscible gas–oil gravity drainage as well as solvent injection process.

Pore-scale experimental studies were also conducted by some researchers, which led to a better understanding of the governing parameters, such as horizontal fracture apertures, wettability and tilt angle, which were reported to be crucial for the stability and shape of the formed liquid bridges (Chatzis et al. 1988; Mashayekhizadeh et al. 2011; Bahari Moghaddam and Rasaei 2015; Khorshidian et al. 2018; Harimi et al. 2019).

Gravity drainage modelling can be categorized as analytical (Bech et al. 1991; Di Donato et al. 2006; Abbasi et al. 2018), numerical (Schechter and Guo 1996; Chow and Butler 1996; Saedi et al. 2015; Udoh et al. 2015; Dejam 2018; Rahmati and Rasaei 2019) and empirical–experimental (Li and Horne 2003; Zendehboudi et al. 2011) models. Some researchers investigated the gravity drainage mechanisms in tight rocks by introducing film flow phenomena (Nenniger and Storrow 1958). Pavone et al. (1989) assumed linear function for relative permeability and the logarithmic relation for capillary pressure and proposed an analytical model. Dejam et al. (2011) studied the impact of new capillary pressure and relative permeability functions to linearize the governing equations of gravity drainage. They calculated the saturation distribution and the production rate using inverse Laplace and Fourier transformation of the linearized equations. Dejam (2018) used a fine-grid numerical model to study the influence of different fracture capillary pressure models on oil recovery and therefore the block-to-block interaction process in fractured porous media. They showed that the existence of fracture capillary pressure results in higher ultimate oil recovery factor. Aghabarari and Ghaedi (2019) introduced a new boundary condition for matrix blocks to include reinfiltration phenomena in commercial simulators.

In this research, free-fall gravity drainage in a three-block system is investigated experimentally and numerically. The main objective is scrutinizing the interactions of the adjacent blocks, particularly capillary continuity and reinfiltration phenomena, which control the efficiency of oil recovery by gravity drainage. For this purpose, a three-block experimental set-up was designed, including vertical and horizontal fractures. Then, an extensive set of experiments were conducted to examine the effects of various parameters, such as, matrix tilt angle, spacers’ permeability, wettability and effective contact surface area on the oil production. To provide a better understanding of the reinfiltration process and block-to-block interactions through visualization of the saturation profiles, numerical modelling of free-fall gravity drainage, including block-to-block interaction processes, was performed by developing a CFD model and verified by experimental data obtained in this study.

This research has features, which make it distinct from other studies of this field. First of all, a multi-block experimental set-up is used to quantify the effect of different characteristics of horizontal fracture on the block-to-block interactions. Second, the effect of horizontal fracture wettability was also investigated on matrix block integrity through gas wetting alteration of the utilized spacers. Moreover, COMSOL (Multiphysics 2012) software is used to model the gravity drainage, which was validated against the experimental data. It is shown by visualizing saturation profiles at different times that the reinfiltration process can also be modelled by COMSOL. The presented experimental results can be used in commercial simulators and for validation of numerical models.

In the remainder of the paper first, the numerical simulation backgrounds and governing equations are presented in Sect. 2. Then, experimental materials and method are presented in Sect. 3. The experimental and modelling results and discussions are presented in Sect. 4, followed by the summary and concluding remarks.

2 Numerical Simulation Backgrounds

Among commercial software packages the multiphysics flexibility of COMSOL (Multiphysics 2012) makes it an appropriate option to study different fluid flow in porous media problems. Bogdanov et al. (2007) investigated steam-assisted gravity drainage mechanism for heavy oil reservoirs using COMSOL. Following that, the process was modelled by STARS package from CMG (Stars and Guide 2008) and the results were compared. They aimed to investigate the adequacy of COMSOL in the modelling of unconventional reservoirs and showed that many important areas of EOR methods can be investigated using this package. Finally, they concluded that, although the computational cost and performance of the commercial packages, like CMG, is still better than COMSOL, due to the high variability and flexibility, COMSOL can be used for phenomenological investigation of new EOR methods in pore to core scale (Bogdanov et al. 2007; Rokhforouz and Amiri 2019). Bjørnarå and Aker (2008) implemented two-phase flow equations of porous media in COMSOL, which can be very complex depending on the introduced assumptions. These equations usually consist two partial differential equations with some auxiliary algebraic equations. The main objective was to determine the appropriate equations for complex modelling (Bjørnarå and Aker 2008). They compared six equations, including partial pressure and fluid flow formulation based on pressure and four other equations were phase pressure–saturation formulation (Binning and Celia 1999), fractional flow (Binning and Celia 1999; Chen et al. 2006), weighted pressure formulation (Chen and Huan 2003; Chen et al. 2006), and Buckley and Leverett (1942) equations based on the pressure and saturation. General form of all mentioned equations is as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} \sum {S_\alpha }=1,\,P_c=P_{nw}-P_n, \, S_{e\alpha }=f(P_c) \\ \frac{\partial }{\partial t}(\phi \rho _\alpha S_\alpha )-\nabla \cdot \left[ \rho _\alpha \, {{\mathbb {K}}} \frac{k_{r\alpha }}{\mu _\alpha }(\nabla P_\alpha +\gamma _\alpha \nabla z )\right] =\rho _\alpha q_\alpha \\ S_{w\alpha }=\frac{S_{w}-S_{wr}}{1-S_{wr}}, \end{array}\right. } \end{aligned}$$
(1)

where \(S_\alpha\) and \(S_{e \alpha }\) denote the saturation and the effective saturation of the phase \(\alpha\), respectively. \(P_c\) is the capillary pressure, nw, w and r subscripts denote non-wetting, wetting phases and residual saturation, respectively. \(\phi\) is porosity, t is time, \(\rho\) is density, \({{\mathbb {K}}}\) is the medium intrinsic permeability tensor (absolute permeability in case of homogeneous and isotropic assumption), \(k_r\) is the phase relative permeability, \(\mu\) is phase viscosity, \(\gamma\) is fluid specific weight, z is the depth and q is the volumetric flow rate. Zendehboudi et al. (2012) modelled free-fall and forced gravity drainage mechanisms using COMSOL, and the simulation results were compared with experimental data. The main objective of this study was to develop a numerical model for gravity drainage in fractured porous media, in addition to describing the effects of various parameters on the process.

2.1 Governing Equations

In this study, the matrix was considered as a homogeneous (constant porosity and isotropic permeability) medium. As the pressure gradient was negligible, fluid density was assumed to be constant over time and space (incompressible flow). All similar mathematical models were developed based on continuity equation, momentum balance (Darcy law in porous media) and capillary pressure functions. The general form of the governing equations for incompressible multiphase flow (without sink and source terms) for non-wetting and wetting phases, in the SI system can be formulated as:

$$\begin{aligned} -\phi \frac{\partial S_w}{\partial t}+\nabla \cdot \left( \frac{{K}k_{rw}}{\mu _w}(\nabla P_w+\rho _w)\right)= 0 \end{aligned}$$
(2)
$$\begin{aligned} -\phi \frac{\partial S_{nw}}{\partial t}+\nabla \cdot \left( \frac{{K}k_{rnw}}{\mu _{nw}}(\nabla P_{nw}+\rho _{nw})\right)= 0 \end{aligned}$$
(3)

where K is the absolute permeability and \(k_r\) is relative permeability. Other used equations for two-phase flow simulation, are capillary pressure and saturation equations, which can be expressed as follows:

$$\begin{aligned} P_c= P_{nw}-P_w \end{aligned}$$
(4)
$$\begin{aligned} S_{nw}+S_w= 1 \end{aligned}$$
(5)

For implementation of the above governing equations into COMSOL software, there is a need for rearrangement of Eqs. 2–5 (Multiphysics 2012). So, the obtained equations can be expressed as follows:

$$\begin{aligned} \phi \frac{\partial S_w}{\partial S_e}\frac{\partial S_e}{\partial P_c}\frac{\partial (P_{nw}-P_w)}{\partial t}+\nabla \cdot \left( \frac{{K}k_{rw}}{\mu _w}\nabla P_w+\frac{{K}k_{rw}}{\mu _w}\rho _w g \right)= 0 \end{aligned}$$
(6)
$$\begin{aligned} \phi \frac{\partial S_{w}}{\partial S_e}\frac{\partial S_e}{\partial P_c}\frac{\partial (P_{nw}-P_w)}{\partial t}+\nabla \cdot \left( \frac{{K}k_{rnw}}{\mu _{nw}}\nabla P_{nw}+\frac{{K}k_{rnw}}{\mu _{nw}}\rho _{nw} g \right)= 0 \end{aligned}$$
(7)

2.2 Capillary Pressure Function

To eliminate time derivatives of saturation and rearranging the equation in pressure form, \(\frac{\partial S_e}{\partial P_c}\) was used (Eqs. 6 and 7 ). Capillary pressure is a function of saturation and is used for saturation profile calculation. There are several sets of capillary pressure and relative permeability functions in the literature such as Brooks and Corey (1964), Brutsaert (1966), Mualem (1976) and Van Genuchten (1980). Brooks and Corey (1964) function have been chosen for this study, which can be expressed as follows:

$$\begin{aligned}&{\left\{ \begin{array}{ll} S_e=\left( \frac{P_{ct}}{P_c}\right) ^\lambda , &{} {P_c>P_{ct}} \\ =1, &{} P_c\le P_{ct} \end{array}\right. } \end{aligned}$$
(8)
$$\begin{aligned}&{\left\{ \begin{array}{ll} K_{rw}=S_e^{\frac{2+3\lambda }{\lambda }}, &{} {P_c>P_{ct}} \\ =1, &{} P_c\le P_{ct} \end{array}\right. } \end{aligned}$$
(9)
$$\begin{aligned}&{\left\{ \begin{array}{ll} K_{rnw}=(1-S_e)^2 \left( 1-S_e^{\frac{2+\lambda }{\lambda }}\right) , &{} {P_c>P_{ct}} \\ =1, &{} P_c\le P_{ct} \end{array}\right. } \end{aligned}$$
(10)

where \(P_{ct}\) (threshold pressure) and \(\lambda\) (pore size distribution index) are characteristic constants of the medium. By defining two separate media in the model (matrix and fracture), with different porosity, permeability and flow function the model geometrical properties were developed. The continuity at the boundary was considered. Capillary pressure value in fractures was assumed nonzero and relative permeability of the fracture was set to be a linear function (Firoozabadi and Ishimoto 1994). Threshold pressure was obtained from Leverett J-function (Leverett 1941), so due to physical property differences, the threshold pressure of the matrix and fracture were different (Sahimi 2011; Zendehboudi et al. 2012; Dejam and Hassanzadeh 2018; Dejam 2018). The numerical model boundary conditions, mesh configuration and validation are further discussed in Sect. 4.2.

3 Experimental Materials and Methods

To understand how reinfiltration and capillary continuity mechanisms affect the gravity drainage process, a three-block apparatus was designed and constructed. It was used to speculate the contribution of the block-to-block interactions, which led to higher ultimate recovery factor. The schematic diagram of the fractured packed multi-block model is illustrated in Fig. 1. The set-up has three inner cylinders of 40 cm height each and the outer cylinder is 120 cm in height. 0.4 mm diameter spots on the inner blocks walls were created using a laser cutter, which connects the matrix blocks to the vertical fracture. The number and diameter of the spots were controlled by the existence of hydraulic continuity between the two media. According to the provided explanation, 400 spots per 40 cm of the model was considered and uniformly performed. Two horizontal fractures were designed between the blocks. The use of spacers represents these horizontal fractures, and as the properties of the spacers changed, the significance of the block-to-block interaction mechanism could be modified accordingly. The utilized spacers and their characteristics are presented in Table 1. The annulus between inner cylinders and the outer cylinder represented the vertical fracture. The flow was controlled by two inlet and outlet valves and the outlet stream accumulated in a glass container, which was closed with a cap to avoid fluid evaporation. Weight of the produced fluid was measured using an electronic balance and the production data were logged by a data acquisition system.

Table 1 Spacers used for the three-block experiments and their specifications
Full size table
Fig. 1
figure 1

Schematic of the stack of three-block experimental set-up

Full size image

Free-fall gravity drainage experiments were conducted at ambient condition (20 °C, 1 atm). During the experiment, the top of the model was open to the atmosphere by opening the inlet valve. The details of the experimental procedure, including sand packing, fluid saturation and data acquisition are as follows:

  1. 1.

    Specific amount of sand (with a size range of 0.2–0.4 mm) was weighed and placed inside of each inner cylinder as the matrix block. The cylinder was placed on a shaker to obtain a uniform and homogeneous sand pack. The process continued until no more sand could be added.

  2. 2.

    The top and the bottom of each inner cylinder were blocked by 400 mesh metal screens. Then, they were placed in the outer cylinder. Afterwards, the upper and lower caps were tightened.

  3. 3.

    The bottom valve was closed and the upper one was opened to the vacuum pump.

  4. 4.

    A glass container of a specific amount of test fluid was placed on balance and the vacuumed system of matrix and fracture begins to saturate using a line, which connected the bottom valve and test fluid in the container. Subtracting the weight of the beaker after and before the saturation process yielded the amount of imbibed fluid throughout the saturation process. having the bulk volume of the system and the void volume using the imbibed fluid, we were able to calculate the porosity of the sand pack (porous media).

  5. 5.

    The experiment commenced when the upper valve was opened to the atmosphere and the bottom one to the glass container. An electronic balance was used to measure the weight of the produced fluid, which was connected to the data acquisition system, used for logging the data. Moreover, the glass container was closed using parafilm to avoid fluid evaporation during the experiment.

The weight of the container was measured continuously, and the recovery factor was calculated from the indirect volume measurements. It should be noted that using crude oil as the test fluid was restricted due to Plexiglas properties. The inner cylinders were representative of the matrix blocks, and specifications of the multi-block model and the test fluids are presented in Tables 2 and 3, respectively. To investigate the role of horizontal fracture surface wetness, four spacers (thin section type I) were coated by MariSeal 800 chemical as a gas wetting agent (Erfani Gahrooei and Ghazanfari 2017, 2018; Azadi Tabar and Ghazanfari 2019).

Table 2 Physical specifications of the three-block experimental set-up
Full size table
Table 3 Test fluids specifications
Full size table

4 Results and Discussion

4.1 Experimental Results

4.1.1 Effect of Tilt Angle

To examine the effect of tilt angle on the gravity drainage process in the presence of block-to-block interactions, three angles of 15°, 6°, as well as 0° from the vertical direction, were selected for the experiments. Moreover, to activate the reinfiltration mechanism, direct contact spacers were placed between the blocks. Figure 2 presents the obtained experimental data for different tilt angles, which shows that increasing the tilt angle of the system leads to a higher ultimate recovery factor. Increasing the tilt angle decreases the effective height of the system so weakens the gravity force. Moreover, in higher tilt angles the travelling liquid blocks will be transferred easier to the vertical fracture rather than forming liquid bridges in the horizontal fracture. So, increasing the tilt angle can decrease the capillary continuity. Noteworthy, the obtained ultimate recovery factors were in the same range of the obtained experimental data by other researchers (Firoozabadi and Markeset 1992, 1994), while quantitative comparison was avoided because of different set-up configuration (such as block height and annulus thickness), rock and fluid properties. The capillary gravity ratio (CGR) for the gravity drainage process can be defined as

$$\begin{aligned} CGR=\frac{\Delta \rho g H}{2\sigma cos(\theta )/r} \end{aligned}$$
(11)

where H is the block height, \(\sigma\) is kerosene-air surface tension (30 dyne/cm), \(\theta\) is kerosene-rock-air contact angle and r is the average pore size. The CGR was calculated as \(\approx\)10 for the experimental condition, which shows the dominance of gravity force over capillary in our experiments. In higher tilt angles the effective block height (H) and the CGR decrease, which results in lower recovery factors.

Mentioned mechanisms are the most significant fluid transfer mechanisms from one block to another. Moreover, the effective height of the system is decreased by the increase of the tilt angle from the vertical plane. In higher tilt angles, the velocity of the travelling droplets decreases and, as a result, the reinfiltration process will be weakened. Thus, in higher tilt angles it is easier for droplets to enter the vertical fracture rather than being imbibed into the lower block, so the liquid bridges become weaker.

Fig. 2
figure 2

Effect of tilt angle on recovery factor in the three-block system

Full size image

4.1.2 Effect of Spacers’ Permeability

The oil recovery process from a block in a fractured reservoir can be divided into two regimes. The period before gravity–capillary equilibrium (i.e., early time, transient) and the period after these forces balance (i.e., late time, steady state). In early time, block-to-block interactions result in a lower gravity drainage rate, due to reinfiltration and the formation of liquid bridges (i.e., capillary continuity) mechanisms. Interaction between blocks increases the effective height and the ultimate recovery factor. The experiments showed that the upper block did not have a significant role in cumulative recovery at the first steps of the process, and it can even postpone the drainage of the lower block. However, by draining the oil from the upper block to the lower, the cumulative recovery of the system is increased. In other words, the drainage of a multi-block system can take more time rather than a single block. Hence, one of the crucial factors in any gravity drainage process is determining the effective height of the blocks in the system.

Hence, it can be concluded that the formation of liquid bridges leads to higher gravity drainage efficiency in late time. It is founded that this mechanism will enhance the drainage rate until the gravity–capillary equilibrium takes place. Overproduction of a multi-blocks system is related to the liquid film phenomenon, which is obtained by direct contact of the blocks, which is controlled by the capillary number and the fracture permeability. Dejam and Hassanzadeh (2011) investigated the critical length of a droplet, which affected the formation of liquid bridges. They showed that the critical length \((L_c)\) depends on a threshold bond number \((B_0^*=\frac{\rho g r^2}{\sigma })\) and fracture capillary pressure \((P_{cf})\). The relation between these parameters is expressed as follows:

$$\begin{aligned} L_c=\frac{\sigma }{P_{cf}}\left( \frac{1}{B_0^*}-2 \right) \end{aligned}$$
(12)

According to Eq. 12, increasing the fracture capillary pressure leads to a decrease in the critical length. To evaluate the effect of capillary pressure which can be related to the pore size distribution and fracture permeability, spacers with different permeabilities were used in the experiments. Figure 3 depicts the effect of spacers’ permeability on the oil recovery factor, in the early time stage the gravity force was dominating the fluid flow. Reducing the fracture permeability, as expressed in Eq. 12 leads to greater critical length and will increase the number of liquid bridges (Dejam and Hassanzadeh 2011, 2018). This phenomenon temporarily decreases the oil production rate. The oil droplets expelled from the upper block will be imbibed into the lower block through liquid bridges. However, in the late time, recovery would increase because of the reduction in threshold height for spacers with lower permeability.

Fig. 3
figure 3

Effect of spacer permeability on the recovery factor in the three-block system

Full size image

4.1.3 Effect of Spacers’ Wettability

Spacers’ wettability can be considered as a parameter, controlling the horizontal fracture capillary pressure. Equation 13 determines the fracture capillary pressure (Kralchevsky and Nagayama 2001; Ghezzehei and Or 2005) as:

$$\begin{aligned} P_{cf}\approx \sigma \left( \frac{1}{r_0}-\frac{2\cos \theta }{b_f} \right) \end{aligned}$$
(13)

where \(r_0\) is droplet radius and \(\theta\) is the contact angle, which represents the wettability of the fracture, \(b_f\) is the fracture aperture and \(\sigma\) is the interfacial tension. Figure 4 shows the schematic representation of a liquid droplet in a horizontal fracture. According to Eq. 13, which applies to an open fracture, the contact angle of the wetting phase, would be between 0° and 90°. Therefore, the value of the capillary pressure is less than the condition in which \(\theta\) is between 90° and 180°. As a result, an oil-wet matrix may result in a stronger capillary continuity and a stable liquid bridge formation, which ultimately results in higher oil recovery. Yet, the equation for the fracture capillary pressure is not applicable in this research, because the utilized spacers were porous (thin sections) and capillary pressure model for this condition can be expressed as:

$$\begin{aligned} P_{cf}\approx \sigma _f \left( \frac{2\cos \theta }{b_f} \right) \end{aligned}$$
(14)
Fig. 4
figure 4

Schematic representation of a liquid droplet in a horizontal fracture (Dejam and Hassanzadeh 2011)

Full size image

Figure 5 illustrates the experimental results of the recovery factor versus time for oil-wet and gas-wet spacers. In each experiment, four thin sections of type (I) were placed on the upper and lower faces of the horizontal fractures (direct contact spacer was used as the horizontal fracture). As depicted in Fig. 5 in early time stage lower recovery was obtained in the presence of oil-wet spacer that may occur due to lower fracture capillary pressure and the corresponding block-to-block interactions, but the ultimate recovery factor was higher in this case; which could be explained as a result of better communication between blocks through liquid bridges between the matrix blocks, and subsequent decrease of the threshold height of the system. Additionally, in the gas wetting condition of the horizontal fractures, there is a higher chance of communication with the vertical fracture for the produced oil from a matrix block as the liquid droplets are not absorbed by the fracture walls, so the chance of channelling to the vertical fracture is increased.

Fig. 5
figure 5

Effect of spacers’ wettability on the recovery factor in the three-block system

Full size image

4.1.4 Effect of Spacers’ Effective Surface Area

Figure 6 presents the obtained experimental results obtained by sparse and direct contact spacers. As shown, direct contact spacers obtained higher ultimate recovery for the system compared to sparse spacers (for the spacers specifications refer to 1). One of the important factors, controlling the reinfiltration phenomena, is the spacers’ effective surface area, providing contact between the upper and the lower matrix blocks. As this factor have a direct influence on capillary pressure, reducing the fracture effective area will hypothetically weaken the capillary continuity. Because in such condition, there is less room for the formed oil droplet to increase in size. Besides, heterogeneity and roughness of the fracture surface hinder the film flow and the travelling droplet mechanism. So, greater contact between the blocks results in higher production.

Fig. 6
figure 6

Effect of spacers’ effective surface area on the recovery factor in the three-block system

Full size image

4.2 Numerical Simulation Results

4.2.1 Model Geometry and Boundary Conditions

Figure 7 illustrates the schematic view of the boundaries for multi-stack blocks numerical model. At initial condition, the model was saturated with oil and boundary#1 was opened into the air and atmospheric pressure and boundary#3 had atmospheric pressure as well. Boundary#2 was non-wet for both phases and no-flow boundary condition was considered at this face. Boundary and initial conditions for the model are summarized in Table 4.

Fig. 7
figure 7

Schematic view of the boundaries for multi-stack blocks numerical model

Full size image
Table 4 Numerical model boundary and initial conditions
Full size table

4.2.2 Mesh Configuration

For a model with a small grid size, the software will provide more accurate results, but the runtime increases for each simulation. Since fluid mostly flows in the vertical direction, a rectangular grid leads to acceptable results compared to triangular mesh. Using the finite element method, which can simply handle the front displacement, fine-grid mapping may not seem to be necessary due to low gradient fluid flow (i.e., free-fall gravity drainage). To ensure the accuracy of the analysis, the effect of grid size on the obtained results was investigated by varying the mesh size to find the optimal size for further calculations.

4.2.3 Model Validation

In this section, model validation against the experimental results, and estimation of the Corey model saturation function parameters are presented. As the model parameters were estimated based on the experimental results and the characteristics of the conducted experiments, some model properties especially the fracture and matrix permeability are different from typical values in the field. Table 5 presents the obtained model parameters.

Figure 8 shows an acceptable match between the experimental data and the numerical model results. The experiment was conducted with sea packed sand and kerosene using direct contact spacers. The horizontal fracture aperture in the numerical model was 0.01 mm and the matrix permeability was 600 Darcy. Table 6 compares the recovery factor results for the experiment and the numerical model for three different times. The presented data illustrate the fact that the numerical modelling results are in good agreement with the experimental data. The average relative error was less than 5% throughout the simulation compared to the experiment, which is acceptable.

Table 5 Multi-block COMSOL model specifications
Full size table
Fig. 8
figure 8

Validation of the three-block numerical model with the experimental data

Full size image
Table 6 Results of experiment and numerical modelling at different times
Full size table

4.2.4 Effect of Mesh Size on the Numerical Solution

Figure 9 presents the numerical modelling results for three different element numbers; the presented results prove that the obtained results are independent of both mesh size and mesh number, in the selected range. So, the coarse mesh (150 elements) was selected for further calculations as it had lower runtime. The mesh was refined at horizontal fractures and we made sure to use at least 5 mesh elements in the fracture aperture.

Fig. 9
figure 9

Effect of mesh size and number on the numerical results

Full size image

4.2.5 Application of COMSOL for Simulating the Reinfiltration Phenomena

Figures 10 and 11 illustrate the oil saturation profiles along the three-block system used in the experiments with uniform block permeability and horizontal fracture aperture of 0.1 mm obtained from the numerical model, at early and late times, respectively. Before any discussion, it may be helpful to briefly review the concepts of reinfiltration phenomenon. Equation 15 presents one dimensional drainage rate at the depth of \(z=0\) (lower surface) and reinfiltration at \(z=H\) (upper surface):

$$\begin{aligned} q=\frac{Kk_{ro}}{\mu _o}\left[ \Delta \rho g-\frac{\partial P_c}{\partial S_o}\frac{\partial S_o}{\partial z}\right] \end{aligned}$$
(15)

The term \(\frac{\partial S_o}{\partial z}\) can be positive, negative or zero. It can be negative when the threshold pressure is zero. At \(z=0\), where the block is fully saturated with oil, \(\frac{\partial S_o}{\partial z}\) is zero. If the threshold pressure is zero, the maximum rate of production at \(z=0\) which is the initial rate reads as:

$$\begin{aligned} q=\frac{Kk_{ro}}{\mu _o}\left[ \Delta \rho g\right] \end{aligned}$$
(16)

When the production from the block at \(z=0\) begins, \(\frac{\partial S_o}{\partial z}\) is negative and the rate of drainage decreases, so the rate of drainage at \(z=0\) follows:

$$\begin{aligned} q\leqslant \frac{Kk_{ro}}{\mu _o}\left[ \Delta \rho g\right] \end{aligned}$$
(17)

If reinfiltration exists at the top of the block, a sufficient amount of fluid can be fed from the upper boundary and \(S_o\) at \(z=0\) equals to 1, as the saturation gradient at \(z=H\) is zero, while if an insufficient amount of fluid was fed from above, \(S_o\) at \(z=0\) equals to 1 and the saturation gradient at \(z=H\) will be positive, so:

$$\begin{aligned} q\geqslant \frac{Kk_{ro}}{\mu _o}\left[ \Delta \rho g\right] \end{aligned}$$
(18)

As it is shown in Fig. 10 the saturation profiles indicate a change or return at 120 s time at the location of the horizontal fracture. The reason for this sudden increase in oil saturation at the horizontal fracture is the emergence of inter-block phenomena such as the formation of liquid bridges and reinfiltration of the produced oil from the upper block and it can increase the oil saturation in the lower block because of the imbibed oil into this block. Since pressure continuity boundary condition was applied at the boundary of the horizontal fractures and the adjacent blocks, the simulator was able to model inter-block phenomena, causing such sudden change in properties (e.g,. permeability, porosity and saturation functions) of two media (matrix block and fracture). In other words, with mesh refinement at the boundary of the fracture and the use of the finite element method, we could study the effect of the block-to-block interactions on the gravity drainage process efficiency. At times of 180 and 210 s feeding the lower block from the upper one through inter-blocks phenomena is obvious from the saturation profile. For instance, oil saturation at 210 s for approximate heights of 0.73 to 0.8 is equal despite oil production. This represents compensation of the oil production with the oil fed from the upper block, being imbibed to the lower one. At times 240 and 270 s the upper block was completely drained and is at residual oil saturation, so the saturation profile of the middle block followed a normal distribution.

It can be concluded that the oil saturation of the upper block should be less than the lower block in both gravity drainage and reinfiltration process; however, these phenomena are more obvious in reinfiltration process because of the imbibed oil into the underneath block.

Fig. 10
figure 10

Early time simulation results of oil saturation profiles for the three-block system. (Note: the horizontal axis shows the model elevation and the oil saturation is shown in the vertical axis)

Full size image
Fig. 11
figure 11

Late time simulation results of oil saturation profiles for the three-block system. (Note: the horizontal axis shows the model elevation and the oil saturation is shown in the vertical axis)

Full size image

4.2.6 Effect of Horizontal Fracture Aperture

In this section, we study the effect of horizontal fracture aperture on the gravity drainage performance, which was experimentally difficult to study. As Fig. 12 shows, decrease in the horizontal fracture aperture stimulates the block-to-block interactions such as reinfiltration, film flow, travelling liquid droplets and capillary continuity through the formation of the liquid bridge; thus, the rate of oil production in early times decreases due to presence of such mechanisms. However, due to capillary continuity, the threshold height decreases and the ultimate recovery is increased. Moreover, when the fracture aperture increases, the capillary pressure in the fracture is decreased which leads to unstable liquid bridges (refer to Eqs. 12–14). Such phenomena decrease the oil drainage rate in early times. Another reason behind increasing the ultimate recovery in late time is lower liquid saturation in the horizontal fractures. In conclusion, stable liquid bridges, which occur in smaller fracture apertures, can maintain the capillary continuity and increase the ultimate recovery factor of the gravity drainage process.

Fig. 12
figure 12

Effect of horizontal fracture aperture on the recovery factor of the multi-stack block system. (All fracture aperture units are in mm)

Full size image

5 Summary and Conclusions

In this work, an experimental set-up, as well as a numerical model of a multi-blocks system, was successfully designed and used to analyse the gravity drainage mechanism in fractured systems. The effects of capillary continuity and reinfiltration phenomena, which mainly control block-to-block interactions, were investigated. The impacts of important parameters such as tilt angle, spacers’ wettability, permeability, effective surface were investigated experimentally using an intensive experimental set. COMSOL multiphysics package was used for numerical modelling of the process to understand the block-to-block interactions through visualization of the saturation profiles. The numerical model was validated by the experimental data. Based on the experimental data and modelling results presented in this work, the following conclusions can be drawn:

  • Decreasing the horizontal fractures aperture decreases oil drainage rate in early time stages but increases the ultimate recovery due to the greater extent of capillary continuity between the adjacent blocks.

  • Tilt angle from the vertical axis decreases the recovery factor, due to higher velocity of travelling liquid element and liquid film flow as well as a higher chance for oil droplets to enter the vertical fracture rather than being imbibed by the underneath matrix block. Moreover, increasing the tilt angle of the system decreases the effective height of the system and lowers the gravity compared to capillary force and hence decreases the ultimate recovery factor.

  • Gas wetting alteration of the horizontal fractures weakens the capillary continuity and decreases the ultimate recovery factor, compared to original liquid wetting condition. Conversely, the model with gas-wet horizontal fractures showed a higher oil production rate in early time stages, due to better communication of the matrix blocks with the vertical fracture.

  • Reinfiltration contribution in the block-to-block interaction process is mainly controlled by the horizontal fracture specifications. Higher spacers’ permeability decreases the oil production rate at early times, while it increases the ultimate recovery factor due to higher capillary continuity and lower system threshold height.

  • Good agreement was observed between the numerical model results and the experimental data of oil recovery, showing that COMSOL multiphysics can successfully be used for the modelling of the gravity drainage process in a multi-block system. The validated numerical model was used to visualize the saturation profiles and study the effect of horizontal fracture aperture on the gravity drainage process in a multi-block system.

  • The behaviour of the saturation profiles along with the matrix blocks, obtained from the numerical modelling, showed that block-to-block interaction phenomenon was successfully modelled by COMSOL CFD software.

The provided experimental data and modelling results shed light into the effect of different fracture properties on the block-to-block interactions and their interplay on the ultimate recovery factor of the system. Yet there are many other important factors which can be further studied by designing specific experimental set-ups or numerical models. Some of these factors include the effect of production rate on block-to-block interactions, the heterogeneity of the matrix blocks and the effect of heterogeneity sequence and layering, fracture roughness and the vertical fracture properties.