Gravity drainage mechanism in naturally fractured carbonate reservoirs; review and application

Gravity drainage is one of the essential recovery mechanisms in naturally fractured reservoirs. Several mathematical formulas have been proposed to simulate the drainage process using the dual-porosity model. Nevertheless, they were varied in their abilities to capture the real saturation profiles and recovery speed in the reservoir. Therefore, understanding each mathematical model can help in deciding the best gravity model that suits each reservoir case. Real field data from a naturally fractured carbonate reservoir from the Middle East have used to examine the performance of various gravity equations. The reservoir represents a gas–oil system and has four decades of production history, which provided the required mean to evaluate the performance of each gravity model. The simulation outcomes demonstrated remarkable differences in the oil and gas saturation profile and in the oil recovery speed from the matrix blocks, which attributed to a different definition of the flow potential in the vertical direction. Moreover, a sensitivity study showed that some matrix parameters such as block height and vertical permeability exhibited a different behavior and effectiveness in each gravity model, which highlighted the associated uncertainty to the possible range that often used in the simulation. These parameters should be modelled accurately to avoid overestimation of the oil recovery from the matrix blocks, recovery speed, and to capture the advanced gas front in the oil zone.


Introduction
Fossil fuels constitute a significant portion of world energy consumption, and the world demand on this type of energy is expected to grow in the coming decades [1,2]. Many major economic sectors are significantly influenced by oil prices. The global oil market has attributed prices as an essential factor to all world economics [3]. However, the currently producing fields are not able to meet the future increment to the fuel energy demand, which requires an increase in production from low-recovery reservoirs or to develop new prospects. Carbonate reservoirs host about 70% of the conventional oil reserves in the Middle East, and most of them are naturally fractured [4]. The Naturally Fractured Carbonate Reservoirs (NFCR) are well known for their low recoveries [5,6] compared to the clastic reservoirs due to the multiscale of geological heterogeneities, [5,7,8]. Therefore, accurate modelling of reservoir heterogeneities is a critical and an essential stage for further improvements in the modelling and simulation process, which enhance historical matching or future prediction. Furthermore, sustaining oil and gas production is key to a successful secondary or tertiary recovery project.
NFCRs are frequently simulated using the dual medium. The matrix, which is often known as a stagnant medium, provides most of the fluid reserves. Meanwhile, fractures provide faster corridors for the flowing fluids. The interaction between the two mediums has a growing interest as it represents the primary pillar of the dual-porosity model. It is represented by different physical phenomena such as diffusion, gravity drainage, imbibition, and capillarity, which are the main recovery mechanisms that control the flowing fluid between the matrix blocks and the surrounding fractures. In a gas-oil system, gravity drainage represents the major recovery mechanism in the reservoir during natural depletion or in gravity-assisted processes [9][10][11]. Therefore, it is necessary to understand the transfer function evolution before evaluating the effectiveness and performance of the various mathematical expressions of gravity drainage using the dual-porosity model. The dual-porosity model has the advantage of requiring significantly less CPU time for the simulation run compared to high-resolution models. However, high-resolution models are superior in representing the reservoir heterogeneity and accurately depicting the gravity drainage process. Therefore, high-resolution models are often used as a reference scenario to evaluate the improvements in the developed formulas for the dual-porosity model (e.g., [12][13][14]).
Saturation profiles and recovery speed are among the challenges that face gravity drainage models. Inappropriate selection of gravity model or misuse of matrix characteristics results in overestimation of oil recovery from the matrix blocks, hence inaccurate representation or matching actual reservoir behavior. Furthermore, several authors (e.g., [15][16][17]) have proposed transfer function formulas to calculate the fluid exchange rate between the matrix and the fractures. Most of these formulas include a gravity term to account for the contribution of gravity drainage in the overall oil recovery from the matrix blocks. The suggested principle by Barenblatt [18] and Warren and Root [19] established the foundation of most of the mathematical expressions of the transfer functions that have Darcy-equation forms. Despite the suggested improvements to the transfer functions, they are still responsible for the accurate mimicry of the recovery speed from the matrix blocks.
Barenblatt [18] described the concept of liquid transfer between the matrix and the fissures, in addition to highlighting the pressure discontinuity (i.e., fluid pressure differences between the fractures and the adjacent matrix cells) as shown in Equation (1).
where τ is the liquid rate that flowed from the matrix cells to the fractures, and α * is a dimensionless characteristic of the fractures that depends on the fluid viscosity µ in addition to the pressure difference between P m , and P f , which represents the liquid pressure in the matrix and fracture respectively. The dimensionless α * in Equation (1) was redefined by Warren and Root [19] to be (k m [L 2 ].σ[L −2 ]), where they defined σ as a shape factor, and k m is the matrix permeability as illustrated in Equations (2) and (3): where n represents the number of the natural fracture sets and L related to the matrix dimensions. It could be highlighted that Equation (2) has the form of the Darcy equation. This equation form represents one of the suggested solutions to simulate the fluid exchanges between the matrix and fractures. However, other approaches (e.g., empirical models) have also been suggested for the same purpose, but they have not been investigated in the current work. The Darcy-form transfer functions have explained in detail as they are available in the commercial simulators (e.g., ECLIPSE) which can apply at both small and field scales.
Kazemi [20] extended the Warren and Root model to the multiphase flow to account for capillary and gravity forces, where he illustrated the multiphase flow in both fractures and matrix mediums as he proposed in his simulator. Moreover, Kazemi explained that the shape factor σ in the extension of Equation (2) for the three-dimensional case of cubic cells could be calculated using the formulas [15]; Further improvement in representing the gravity effect was suggested by Gilman and Kazemi [16] by using the subdomain approach, where they demonstrated that using matrix grid refinement could rigorously capture the saturation profile and more accurate representation of fluid exchange between the fracture and matrix subdomain; see Figure 1. The subdomain approach suggested by Gilman and Kazemi to improve the gravity modelling which shown the matrix grid refinement in addition to the complete fluid segregation in the surrounding fractures when the capillary pressure is zero, after Gilman [16].
They suggested that the potential differences between the matrix subdomain and fracture can be calculated using the equation: where the s refers to the subdomain system and α refers to the fluid phase. The pressure in the fractures P f s can be calculated in an oil-water system according to the following formulas: The term P f I in Equation (6) refers to the pressure at the interface between the oil and the water phase while the (∆z/2) represent the gridblock center. Gilman and Kazemi validated their improved gravity model by comparing its results with fine-grid models.
Quandalle and Sabathier [17] referred to the non-equal influence of the three active forces of recovery (viscosity, gravity, and capillarity) that should be adjusted, where three coefficients were defined (K v , K g , and, K c ) to represent the former mentioned active forces which should be multiplied by the matrix-fracture exchange to adjust the flow for the corresponding force. Abushaikha and Gosselin [21] derived the transfer function proposed by Quandalle and Sabathier [17] as each matrix block face has its own potential, mobility, and permeability. The equation has the following form, where the vertical transfer function is separated from the horizontal transfer function: Moreover, Abushaikha and Gosselin [21] explained that the horizontal shape factor (σ Q ) and vertical shape factor (σ GD ) may be defined as follows: Gilman and Kazemi [16,22] proposed a modification to the Sonier [23] formula to improve the representation of the gravity effect by replacing the gridblock thickness (∆z) with matrix block thickness (L z ), in addition to adding the depth for both fracture and matrix to approximate the gravity segregation in the fractures as illustrated in the following equation: where h f α and h m α are the phase α saturation height in the fracture and matrix respectively, which could be calculated using the formula for both fracture and matrix: Uleberg and Kleppe [24] referred to the importance of capillary pressure continuity in the matrix blocks that significantly impacted oil recovery, especially in the gas-oil gravity drainage process. Furthermore, they highlighted the importance of other parameters such as the shape of the matrix block, re-infiltration process, and diffusivity, which required accurate representation to depict the physical flow phenomena occurring between the fracture and matrix.
Alkandari [25] worked on scaling up the centrifugal experiments into the reservoir scale in a gas-oil system to evaluate recovery under gravity drainage. He concluded a hyperbolic decline of the oil recovery due to gravity drainage and he recommended construction of the transfer function using the oil recovery curves of the experimental results. Balogun [26] used several scenarios to demonstrate the effect of different recovery mechanisms and to improve the transfer function by matching the results of a fine-grid model. Abushaikha and Gosselin [21] studied the effect of the shape factor using different matrix shapes and they examined several sensitivity scenarios, in addition to comparing and evaluating different forms of shape factors and transfer functions including the gravity term.
Su [27] suggested using a dynamic shape factor derived from a comparative study using fine-grid model results. The author proposed a simple correlation to fit the shape factor by using fitting parameters as illustrated below: where α * * and β in Equation (12) are fitting parameters that adjust the dual-porosity model results to match the fine-model outcome. Furthermore, as the formula illustrates, these results are limited for specific vertical matrix block dimension L z of 6.34 m (20.8 ft) and it is required to be adjusted for different matrix block heights. Despite the good results obtained by the new suggested shape factor, it was only able to match the recovery profile of the uppermost cell of the vertical stacked cells 1D model, while a block-to-block connection is used for the rest of the cells. Moreover, the error evaluation of the model validation cannot be supported as it was based on the same fine model that the improved shape factor was fitted to in order to match its results. This was in addition to the difficulties of determining the uppermost cells in the real reservoir model, which were made using their suggestions that were impractical for real case studies. Other transfer functions with non-Darcy-form equations (e.g., [28,29]) have not been tested in the current work due to their unstable performance.
The transfer function has been intensively studied to represent the physical parameters accompanying fluid exchange (e.g., gravity, capillary pressure, diffusion) that should be accurately accounted for in order to improve the calculation of the dual-porosity model. Despite the abovementioned progress in modelling the exchange rate between the matrix and fractures for the dual-porosity model, it is still impossible to capture the transient flow effect as a single value only of fluid saturation assigned to the matrix. Moreover, Gilman and Kazemi's suggestion [16] of matrix refinement resulted in an overwhelming number of cells that can be equal to the number of cells in the fine model that required significantly higher computation time compared to the original dual-porosity model and nullified the advantage of the model. However, Quandalle's model [17] is preferable due to separation of the vertical flow potential from the horizontal flow, as the gravity drainage works mainly in the vertical direction, which reduces oil recovery overestimation and improves the model outcomes due to elimination of the gravity effect in the horizontal flow direction.
The current work has two objectives. The first is to illustrate the importance of the gravity drainage mechanism in the gas-oil system using fine-scale modelling as well as investigating the consequences of using various mathematical models on the field performance and saturation profiles across the reservoir. This objective was achieved using various modelling scales-fine-, medium-, and reservoir-scale models-to demonstrate these differences. The second is to evaluate the impact of gravity drainage parameters, such as matrix block dimensions and matrix permeability, on the speed and amount of the oil recovery. The effect of the abovementioned parameters has been shown through several sensitivity scenarios to compare the oil recovery for an acceptable range of each parameter. Furthermore, this work helps to gain an understanding to improve the modelling of the gravity drainage mechanism, and carefully select the matrix characteristics that control the oil recovery and recovery speed from the matrix blocks.

The Recovery Mechanism of the Gravity Drainage-Model Comparisons and Parameters Sensitivities
Despite the presented improvements in some mathematical formulas of the gravity drainage process or the transfer function in general. It would be impossible to apply these equations at the field scale. Furthermore, their performance and stability have not been verified at this scale. Therefore, the available options in the ECLIPSE simulator have been used to compare the behavior of different gravity models and their consequences on the traditional history-matching parameters such as Bottom Hole Pressure (BHP), Gas-Oil Ratio (GOR), and oil production rate. Furthermore, sensitivity cases have been implemented to evaluate the role of both matrix permeability in the vertical direction (k z ) and the height of the matrix blocks (L z ) on the ultimate oil recovery.
The current work was divided into three evaluation steps. The evaluation started from the medium scale to illustrate the differences between various gravity drainage models. Then, a fine grid was constructed to examine the fluid exchanges between the matrix blocks and fractures driven by gravity effect due to density differences between the gas and oil until the system reaches to an equilibrium state. Finally, there was an investigation into the gravity drainage effect at a full-field scale using a real case of a fractured reservoir.

The Effect of the Gravity Drainage Medium-Scale Modelling
The second 2D model of the 6th SPE comparative solution project [30] has been used to illustrate the effect of the different gravity drainage expressions that were mentioned previously on the ultimate recovery from the matrix gridblocks and to estimate the remaining oil saturation in the matrix cells. The 6th SPE comparative solution project consists of a simple cross-sectional dual-porosity model which comprises five layers with uniform thickness as shown in Figure 2. Further model details, dimensions, and rock and fluid properties can be found in [30]. A natural depletion scenario was simulated in the model using a single well producing from the fifth layer. The following comparison and sensitivity cases have implemented to evaluate the gravity formulas effect on the model performance:

Comparison of the Gravity Drainage Models
There are several gravity drainage formulations suggested for evaluating the fluid exchange between the matrix and fracture domain due to gravity force. Light is shed on some of the gravity expressions in Equations (8) and (10) where, in these equations, the oil is allowed to flow from the matrix to the fractures but the reverse flow of oil is not allowed (e.g., [15,17]). Moreover, Quandalle and Sabathier formulations can be set to allow the flow of oil from the fractures to the beneath matrix block (i.e., re-infiltration as explained in Figure 3) to represent the fluid flow between the two domains due to gravity. In the 6th comparative solution project, the Kazemi gravity model was used in the project, and the results were compared to the other results of various gravity drainage models, as listed below: The simulation results of the abovementioned gravity drainage options (A, B, C, and D) for the 6th comparative solution project are illustrated in Figure 4. The simulation outcome exhibited remarkable differences in model behavior, which highlighted the role of the gravity mechanism. Oil recovery could be a clear example to demonstrate that oil recovery increased by approximately 10% when the gravity option switched from No Gravity to Alternative Gravity Drainage. The oil allowed to flow from the matrix to the fractures as illustrated by black arrows (e.g., [15,17] models), while the allowed re-infiltration option is demonstrated by red arrows. Comparison of the model behavior using different gravity drainage options, remarkable differences were observed in the results which highlight the significant roles of the gravity mechanisms. The results of the gravity drainage scenario (blue line legend) represent the published results of [30] work, where the data file is available in the ECLIPSE package with the name SPE6_FRAC.
Furthermore, the performance of each gravity drainage formula can be demonstrated by comparing the oil and gas saturation distribution at the end of the simulation time as shown in Figure 5. The three scenarios (No Gravity, Gravity Drainage and Alternative Gravity Drainage-allow re-infiltration) have demonstrated a lower performance by lower oil recovery (Figure 4), high gas saturation in the fractures, and poor matrix block recovery compared to the Alternative Gravity Drainage scenario.

Figure 5.
Oil and gas saturation profile at the end of simulation time, which illustrates how different gravity drainage options change the saturation profile. Furthermore, the saturation profile also highlighted the underestimation of the exchange rate between the matrix and the fracture domain when the gravity effect was not accounted for, or when inappropriate representation of the gravity mechanisms was used.

Vertical Permeability of the Matrix
As highlighted in Equations (8) and (9), Quandalle and Sabathier [17] have suggested division of the shape factor into two forms: • Vertical shape factor (σ GD ): used in the fluid exchange between the matrix and fracture in the vertical direction that significantly affected by the gravity force. • Horizontal shape factor (σ Q ): used in the fluid exchange in the horizontal direction.
Therefore, the vertical permeability can be used in this formula to accurately represent the exchange rate in the top and bottom faces. However, other suggested formulas use the horizontal permeability only for all flow directions. The following results illustrate the effect of vertical matrix permeability on the ultimate oil recovery, where a range of vertical permeability of the matrix (0-10 mD) was tested to investigate its effect by selecting the option of Quandalle and Sabathier expression (Alternative Gravity Drainage- Figure 6A) and Gilman and Kazemi Expression (Gravity Drainage- Figure 6B). The obvious conclusion of the simulation outcomes is that the vertical permeability of the matrix has no effect on the oil recovery in the Gilman and Kazemi formula, which means an inappropriate representation of the fluid exchange in the vertical direction by using the horizontal matrix permeability instead of the vertical permeability for all gridblock faces. Furthermore, Abushaikha and Gosselin [21] have compared several transfer function formulas in their work using single-cell models, and they have concluded that the Quandalle [17] formulas are the best representations of gravity drainage compared to the other available formulas (e.g., [15,20,28]).
2.1.3. Matrix Block Dimensions (L x , L y , and, L z ) Actual matrix block dimensions (which are different from the grid dimensions as illustrated in Figure 7) have an impact on the values of both shape factor and gravity drainage term, hence on the transfer function. The matrix block dimensions are related to the fracture frequency or fracture intensity (i.e., number of fractures per line length, scanline measurements). The scanline measurements of the Fracture Intensity (FI) could be converted into average fracture spacing (S mean ) by using the formula S mean = 1/FI [31]. The average fracture spacing represents the average matrix block size in the measured direction. However, the matrix block height could be related to the mechanical layer thickness in strata-bound fractures, where the fractures develop independently of other geological units [32]. Natural fractures may occur at various spaces which create various sizes of matrix blocks, [6]. These sizes, which represent the actual dimensions of the matrix block, are preferably considered in dual-porosity and dual-permeability models. However, in the modelling workflow, the matrix block dimensions that resulted from the upscaled fracture network are typically used, which represents the local average of the fracture spacing. Moreover, the matrix block height (L z ) is one of the essential parameters in gravity drainage terms in the previously mentioned formulas (e.g., Equations (8)-(11)). Therefore, various heights of the matrix blocks have investigated using Quandalle and Sabathier formula [17] and by using multipliers to maintain the variation of the originally distributed heights in the grid model as illustrated in Table 1 , while the simulation results are shown in Figure 8:

The Effect of the Gravity Drainage Fine-Scale Modelling
A small-scale model was constructed to simulate the gravity drainage process using ideal matrix blocks with the dimensions 3.048 × 3.048 × 3.048 m 3 (10 × 10 × 10 f t 3 ) surrounded by fractures with 0.0762 m (0.25 ft) width, as illustrated in Figure 9. Each matrix block was refined into ten cells in the I, J, and K directions. The total number of matrix cells is 8000 cells. The simulation of the fine-grid model aims to illustrate the fluid exchanges between the matrix and fractures due to the gravity differences of the contained fluids that are solely attributed to the gravity forces. Furthermore, the effect of matrix block size on the ultimate oil recovery from the matrix by drainage process was investigated in addition to the effect of matrix permeability in (x, y) directions and (z) direction.
To simulate the oil recovery from matrix blocks under the gravity effect, the fractures have assumed to be fully saturated with gas while the matrix blocks contain oil saturation in addition to irreducible water saturation (S wirr = 0.22). These conditions were provided to mimic the gas advances in the fractures, leaving the oil trapped in the matrix cells. A unity slope of the relative permeability curves was used for the fractures besides zero capillary pressure. However, a relative permeability set of Qamchuqa formation (which is discussed later) was used for the matrix medium. Further rock and fluid properties are illustrated in Table 2 below:

Full-Field Application and Sensitivities
One of the Middle East's naturally fractured carbonate reservoirs was used to evaluate the effect of different gravity drainage formulas on the transfer function at the field scale. The reservoir is known as the Cretaceous Qamchuqa reservoir and is one of the producing reservoirs in the Jambur field. Jambur field is in the foothill zone to the northeast of Iraq; see Figure 10A. It is represented by a simple structure of asymmetrical anticline with northwest-southeast axis direction. The field extends over 30 km length and 4.5 km width. The Cretaceous reservoir is the reservoir of interest in the current investigation, and it consists of two producing formations known as Upper Qamchuqa and Lower Qamchuqa, which are considered to be highly prolific formations in Iraq. The Qamchuqa terminology has often used in the north of Iraq. However, their equivalent formations in the south of Iraq and the Arabian Peninsula are known as Mauddud and Shauiba formations, respectively, [33][34][35]. Qamchuqa formations consist of thick-bedded limestone [34], deposited in a shallow marine setting. The two Qamchuqa formations are separated by Upper Sarmord formation, which is a non-reservoir shaly unit; see Figure 10B. The reservoir thickness is estimated to be approximately 600 m in the main reservoir region, which represents the Neritic deposition environment, as illustrated in Figure 11. The thickness dramatically decreases toward the southern region to less than 160 m. However, the Cretaceous reservoir does not exist in the Basinal deposition environment where Qamchuqa facies (i.e., reservoir rocks) passes into the basinal non-reservoir facies of Balambo, [34], Figures 10B and 11. A history-matching scenario was implemented for the Cretaceous reservoir using the gravity model of Gilman and Kazemi [15]. Very good matching results were reported for several producers in the reservoir for oil production rate, bottom hole pressure (BHP), gas-oil ratio (GOR), and water cut. However, no evaluation was conducted regarding the oil recovery from the matrix cell that has been surrounded by the advanced gas in the fractures.
The current investigation aimed to compare the different gravity drainage formulas and to evaluate their performance in addition to highlighting the saturation profiles of both oil and gas in the matrix blocks. Furthermore, sensitivity scenarios were implemented to investigate the effect of matrix vertical permeability and matrix block height. This evaluation highlighted the significant uncertainty in estimating the remaining oil reserve behind the advanced gas-oil contact (GOC) level that could be a target for the secondary or tertiary recovery processes, and it helps to make a sound decision for such costly development options.

Fine-Scale Modelling
A five-year scenario was applied by allowing the gravity segregation to occur between the gas and the oil due to the density differences between the gas and the oil. The gravity-driven flow allowed the fluid exchange between the matrix and the fractures. The gradual reduction in the oil saturation in the matrix represents the oil flow out of the matrix blocks toward the fractures and substituted by the counter flow of the gas as illustrated in Figure 12, which shows the time-lapse of oil saturation changes in the matrix blocks for five years. The oil recovery from the matrix blocks depends on several factors such as matrix block dimensions, vertical permeability, and relative permeability curves that control the speed and recovery factor. Four different sizes of matrix blocks were investigated (1.524, 3.048, 6.096, and 15.24 m (5, 10, 20, and 50 ft)) where the matrix dimensions are assumed to be equal (lx = ly = lz). Significant differences in recovery speed were observed as well as ultimate recovery that results from use of various block dimensions, (e.g., 40% of oil recovery from the matrix required 167 days for matrix dimensions (1.524 × 1.524 × 1.524 m 3 ), while it required 1795 days for the matrix dimensions (15.24 × 15.24 × 15.24 m 3 ). Therefore, inappropriate assumption of the matrix blocks during the fracture modelling workflow can accelerate or degrade the recovery speed and value as shown in Figure 13.
The matrix permeability, on the other hand, has demonstrated contrasting responses. Horizontal permeability in both directions (x and y) has tested for the values (0.01, 0.1, 1, 5, 10, and 20 mD), as illustrated in Figure 14. The observed results have shown negligible differences for the recovery speed where all the scenarios overlapped each other, and they have achieved the same ultimate recovery Figure 14A. However, the results of vertical permeability scenarios have shown remarkable differences in both oil recovery and recovery speed (e.g., the oil recovery from the matrix increased by 31% when the vertical permeability increased from 0.1 mD to 1 mD).

Comparison of the Gravity Formulas
The previously investigated formulas (gravity drainage, alternative gravity drainage, and alternative gravity drainage-allow re-infiltration) were applied to the Qamchuqa reservoir to evaluate their performance at the field scale. The traditional comparison of field performance was examined, as illustrated in Figure 15 for the BHP, where insignificant differences were observed for the three formula outcomes (the pressure curves overlapped above each other) except the scenario where no gravity effect was considered and a lower pressure response (the blue line in Figure 15) was obtained compared to the observed data. As highlighted in the 6th SPE comparative solution project, the gravity drainage term in the Quandalle and Sabathier formula uses the vertical permeability in the vertical fluid exchange calculation for the upper face and lower face of the matrix gridblock coupled with vertical shape factor (σ GD ), where the matrix recovery due to the gravity drainage is a function of the matrix vertical permeability and matrix block height (L z ). A cross-section from the field model was used to compare the saturation profiles for both oil and gas using different gravity drainage models as its location is shown in Figure 16A, while the vertical permeability of the matrix is shown in Figure 16B to compare with the saturation profiles in the reservoir using Quandalle and Sabathier models. Under the same conditions, the simulation results referred to notable differences in the oil recovery from the matrix blocks and in the gas saturation increases, as shown in Figures 17 and 18. This dramatic variation in matrix recovery is related to the gravity model performance and how each model can accurately address the gravity drainage effect. The previous conclusion of [21], that the Gilman and Kazemi model [15] (i.e., gravity drainage option) is overestimating the speed of recovery, can be seen in the current evaluation at the field scale; see Figure 19.
In Gilman and Kazemi model, the oil saturation has depleted uniformly along the cross-section, while the gas saturation increases firmly; see Figures 17B and 18B. However, the alternative gravity drainage models (i.e., Quandalle and Sabathier model [17]) have shown almost similar performance and the drainage variation of the fluid saturation in the matrix blocks. Figures 17C,D and 18C,D have a strong relationship with the vertical permeability of the matrix; see Figure 16B. Meanwhile, the no-gravity-effect results were provided for comparison purposes; see Figures 17A and 18A to illustrate the effect of the gravity drainage contribution in general.
Further illustration was made by quantitative evaluation of the oil saturation in the matrix blocks by comparing the formerly mentioned gravity drainage options. A single block (15 130 11) was selected from the cross-section of Qamchuqa reservoir to demonstrate the variations in the oil recovery from the matrix block as illustrated in Figure 19. Where the reduction in the oil saturation is related to both oil production and pressure depletion, that causes the liberation of the dissolved gas, and changes the relative fluid volumes within a gridblock.

The Vertical Matrix Permeability (K z )
Vertical permeability is one of the most important properties of the reservoir that should be characterized accurately especially when the reservoir heterogeneity cannot be neglected [36]. Furthermore, the diagenesis processes in the carbonate reservoirs alters the depositional properties of the rocks significantly and creates a faster path for the fluid flow, such as solution channels or flow baffles due to cement precipitations. Therefore, vertical permeability has a significant role in the fluid exchange between the matrix and fracture system in the naturally fractured carbonate reservoirs.
Sensitivity scenarios for the vertical permeability effect on the oil recovery from the matrix blocks using the alternative gravity drainage model [17] have investigated. A multiplier factor was used to examine the effect of different values of vertical permeability, and to maintain the variation in the distributed property, which is populated based on the geological characterization of the reservoir. The simulation results illustrate that the higher the vertical permeability, the higher the oil recovery from the matrix blocks, as shown in Figures 20 and 21.  The alternative gravity drainage model of [17] was used to evaluate the effect of the matrix block height (L z ) on the oil saturation in the matrix blocks hence the oil recovery. Modifying the matrix block height will change the vertical shape factor (σ GD ) and affect the gravity mechanism. The simulation results indicate that for a smaller matrix block height, the oil saturation significantly decreases while the taller matrix blocks exhibited a slower oil recovery as shown in Figure 22. Furthermore, the matrix block height shows negligible differences when reaching a certain height where if (L z ) is multiplied by 5 or by 10 in this particular reservoir. The results almost overlapped and no considerable increment in the recovery is observed (cf. [37]). Therefore, overestimating the matrix block height after this point will have no effect on the results due to the capillary pressure threshold between the matrix and fracture becoming negligible [37].

Discussion
The effect of the gravity drainage mechanism has been illustrated through different scenarios using different modelling scales-fine grid, medium grid, and in a real naturally fractured carbonate reservoir model. Various formulas exhibited different abilities to represent the effectiveness of the gravity drainage process accurately. The effect of the gravity term in the Gilman and Kazemi equations tend to overestimate the matrix recovery in the early time (see Figure 19), because they use the horizontal permeability for gravity drainage calculation compared to Quandalle and Sabathier equations, which use the vertical permeability. Furthermore, the Quandalle and Sabathier suggestion of using two different shape factors have improved the representation of the fluid exchanges across each face of the matrix blocks and better representation of the gravity drainage in the vertical direction.
Although the selected gravity drainage models using the dual-porosity model hugely affect the saturation profile in the matrix block, other parameters have also exhibited a comparable role. These parameters are vertical permeability of the matrix, when Quandalle and Sabathier formulas are activated, and matrix block heights which are among the matrix characteristics that can change the saturation profile tremendously. The sensitivity scenario results implemented by altering (k z ) and (L z ) demonstrated a wide range of changes of the fluid saturation distribution in the matrix blocks across the reservoir as shown in Figures 17 and 18. Further confirmation was illustrated through the fine-scale modelling results, where similar results were concluded using the conventional single-porosity model. Therefore, it is essential to address and evaluate the impact of such parameters. The geological characterization of the reservoir can be used as a guidance to reduce the tested range of for (k z ) and (L z ) values and to narrow the uncertainty in the simulation outcome.
Tuning the history-matching results in the naturally fractured reservoirs can be quickly achieved by switching between the gravity drainage options or by altering the vertical permeability or the height of the matrix blocks. Nevertheless, the selected options or values should be justified and tested to prove their validity and to sustain the model reliability for future prediction, where the recovery value and recovery speed can be highly biased and unreliable as illustrated in the fine-model outcomes when inappropriate parameters used. Moreover, the success of any secondary or tertiary recovery methods requires simulation models with reliable prediction ability, where recovery speed using the dual-porosity model is very sensitive to the reservoir characterization that can substantially impact the ultimate recovery and nullify the efforts of increasing the recovery factor from such fields [38].
Although one example of a naturally fractured reservoir was used (i.e., Qamchuqa reservoir) the discussion applies to all the naturally fractured reservoirs where the contrast between fracture and matrix flow capacity make the gas advances in the fractures always ahead compared to the matrix. Gas advancement in the fractures leads to activating the gravity role and accelerating the fluid exchange between the fracture and matrix and to contribute remarkably to the recovery mechanism. Therefore, similar results are expected if other data sets of a fractured reservoir are used in the analysis.
In vuggy fractured reservoirs, the reservoir performance can be varied and depends on the characterization of the vugs. Separate vugs increase the reservoir porosity, but it does not enhance the matrix permeability [39,40]. The touching vugs in the matrix enhance the low matrix permeability up to ten times of that expected from the matrix [39], while improving the exchange rate between the matrix and the fractures; hence the gravity drainage mechanism. Moreover, the fractures could connect the vugs in the system and increase fracture storage and flow capacity [40][41][42]. Furthermore, higher touching vug density in the reservoir may behave equivalently to the fracture performance, which improves the matrix connectivity to the fracture system significantly.

Conclusions
The gravity drainage process is one of the significant drive mechanisms in naturally fractured reservoirs, especially in the gas-oil system as a natural depletion or as a gravity-assisted process [9][10][11]. The contribution of the process has been evaluated qualitatively through the saturation profiles and the recovery speed. However, a quantitative evaluation of the gravity drainage process can be achieved