Modelling the transient effect in naturally fractured reservoirs

Modelling of Naturally Fractured Reservoirs (NFR’s) is a very challenging task. NFR’s are often simulated using the dual-porosity (DP) model, which requires significantly lower computational time and less simulation cost compared with the fine-grid modelling. However, the DP model is unable to capture the transient effect, the saturation front, and its tendency to overestimate the predicted oil recovery, as the invaded fluid is immediately reaching the gridblock centre once its saturation exceeds the critical saturation. The matrix block discretisation and modifying the transfer function are among the widely investigated areas to represent the transient effect and to improve the simulation accuracy. Adjusting the transfer function often results in complicated and unstable solutions, which make this approach limited in use. However, the matrix discretisation technique, such as Multiple Interacting Continua (MINC), is one of the utilised approaches to improve the simulation of NFR’s. The model’s layers provide the required mean to capture the transient effect and to include the matrix heterogeneities. In this work, we present an improvement to the original MINC model to enhance its accuracy and stability using a particular case of two subdomains model. We suggested using an equal volume of the matrix sub-blocks, besides performing an adjustment to the calculation of fluid saturation. The new adjustments have provided a stable solution and improved simulation results compared with the original MINC model. Besides the advantage of the matrix layers to handle various heterogeneities with excellent simulation accuracy and marginal errors, and hence a reliable performance and prediction.


Introduction
Naturally fractured carbonate reservoirs are very heterogeneous systems. The matrix heterogeneity in NFR's is related to the depositional environment and the subsequent diagenesis processes that alter the depositional properties of the rocks (e.g. (Aljuboori et al. 2020a)). Meanwhile, the pervasive occurring of fractures in most of the carbonate rocks creates additional heterogeneity to the system due to their multiscale and their various flow capacities (Agada et al. 2013;Nelson 2001). The compound heterogeneity of the two mediums makes the reservoir characterisation and fluid flow modelling even more complicated compared with the conventional unfractured reservoirs.
The dual-porosity model has been used widely to simulate the fluid flow behaviour in fractured reservoirs due to its convenience to represent fractures, besides it requires substantially less computational time compared with the fine-grid models (Gilman and Ozgen 2013;Ramirez et al. 2009). In the dual-porosity model, the flow occurs between the gridblocks in fractures only, while the matrix medium is discontinuous and considered as a source term to support the fractures (Gilman and Ozgen 2013). The interaction between these mediums, which is known as matrix-fracture transmissibility, fluid exchange or transfer function, has a significant impact on the ultimate oil recovery from the matrix medium (Aljuboori et al. 2019;Elfeel et al. 2013). The dual-porosity model assumes a quasi-steady-state flow between the matrix and fractures (Coats 1989;Warren and Root 1963). Therefore, the model is unable to represent the transient effect, and hence inaccurate recovery prediction that continues from the early time to the late time (Ahmed Elfeel et al. 2016). The model limitations have extended to inaccurate capturing of the fluid saturation and the moving front of the invaded fluids in the reservoir (Abushaikha and Gosselin 2008;Aljuboori et al. 2019).
In previous decades, numerous attempts were performed to improve the performance of the dual-porosity model by proposing several adjustments to the transfer function formula, the calculation of shape factor and the discretisation technique in the matrix block (e.g. (Alkandari 2002;Geiger et al. 2013;Gilman and Kazemi 1988;Heinemann and Mittermeir 2016;Kazemi and Gilman 1993;Lu et al. 2008Lu et al. , 2006Quandalle and Sabathier 1989;Sarma and Aziz 2006;Sonier et al. 1988;Su et al. 2013)). Furthermore, the outcomes of the proposed adjustments have frequently been compared and calibrated with a reference solution of fine-grid modelling. Despite the success in several solution approaches, the complexity and the solution stability are among the critical challenges that reduce the chances to apply them globally at a sector or reservoir scale.
Matrix block discretisation or refinement is one of the excellent solutions to enhance the accuracy of the dualporosity model, see Fig. 1(A and B) (Gilman and Kazemi 1988;Kazemi and Gilman 1993). However, this discretisation scheme often results in an overwhelmed number of cells that required significant computational power. The same idea has been applied using a different discretisation style that results in a minor increment in the matrix cells. This model is known as Multiple INteracting Continua (MINC), which was proposed by Pruess and Narasimhan (Pruess and Narasimhan 1982), and Pruess (Pruess 1992(Pruess , 1985. In this discretisation style, the matrix block is divided into overlapping shells called subdomains (or sub-blocks), as illustrated in Fig. 1C. The matrix subdomains provided the required mean to capture the transient effect through having their independent properties such as; pressure, saturation, chemical concentration, or any other characteristics, which help to improve the simulation accuracy of the dual-porosity model (Al-Rudaini et al. 2020

Matrix Sub
Fracture Block Fig. 1 The layout of a gridblock in the dual-porosity model illustrates A The matrix block and the fracture block B The refinement of the matrix block and the fracture block C The overlapping matrix shells and the fracture block 1 3 The number of subdomains in the MINC model depends on the process that is required to be simulated. Al-Rudaini et al. (2020) have investigated various scenarios of the matrix subdomains division to enhance the simulation accuracy of chemical-component transport by comparing the minimum error between MINC, the classical DP model and the Single Porosity (SP) model. Furthermore, a simpler MINC model of two matrix subdomains only has been used by Ahmed Elfeel et al. (2016) to simulate Water Alternate Gas (WAG) injection scenarios. The obtained results of these two sub-blocks models have referred to significant improvements in the simulation accuracy compared with the classical dual-porosity model. Moreover, simplifying the assumptions of the tested models might affect the comparison outcomes of the DP and MINC models with the reference solution using the SP model. Some of these simplifications comprised using an equal viscosity and density for the oil phase and the water phase to identify the impact of spontaneous imbibition solely without the gravity effect (e.g. (Al-Rudaini et al. 2020;Ahmed Elfeel et al. 2016)).
In this work, we adapted the illustrated MINC model in Fig. 1C using two matrix subdomains. An equal matrix volume was used for both matrix subdomains in addition to a correction term to adjust the calculation of water saturation for the outer matrix sub-block. Moreover, heterogeneous properties of the matrix sub-blocks have been utilised to include the variation in the matrix medium such as; permeability, porosity, initial water saturation, saturation exponents of the oil phase and water phase, and the capillary pressure exponent.
The dual-porosity model, evolution and limitations Barenblatt et al. (1960) were the first who described the concept of liquid transfer between the matrix medium and fractures. They had also referred to the pressure discontinuity between the two mediums (i.e. pressure differences between fractures and the adjacent matrix cells). This concept was later redefined by Warren and Root (1963), which they illustrated the calculation of the transfer function using Eq. (1) under quasi-steady conditions: where; τ Transfer function, σ Shape factor, µ Viscosity, k m Matrix permeability, P m , P f Pressure in the matrix and fractures, respectively.
In this flow condition, the transient effect occurs very fast for a single-phase, isothermal and slightly compressible fluid (Pruess 1992). However, the quasi-steady flow assumption is not valid for the two-phase flow, which usually occurs in the producing reservoirs. The transient effect in the two-phase (1) = k m P m − P f flow becomes significant, and it cannot be neglected. In such a case, the classical dual-porosity model fails to capture the actual reservoir behaviour such as fluid saturation, invaded fluid front, and recovery speed. The quasi-(semi-) steadystate condition assumed that the saturation front of the invading fluid reaches the block centre immediately. While the saturation in the fine-grid model takes the shape of the capillary pressure curve, in which the invaded fluid moves until its front reaches the centre of the sub-divided matrix block (Abushaikha and Gosselin 2008). The dual-porosity model using the transfer function of Gilman and Kazemi (1983) overestimates the exchange rate between the matrix blocks and fractures, hence exaggerating the oil recovery speed and amount from the matrix medium. The amount of fluid exchanged in an oil-water system between the matrix and fracture can be estimated using the given equation below (Kazemi and Gilman 1993): where; τ w Transfer function of water phase, V Matrix block volume, k m Matrix permeability, k rw Relative permeability of water phase, B w Water formation volume factor, P w m , P w f Pressure of water phase in the matrix and fractures respectively, Δρ Density differences between water and oil or gas, g gravitational constant.
The parameters (h w m ) and (h w f ) represent the saturation heights of water in the matrix and fracture respectively, which can be calculated using the following formula for both fracture and matrix:where; S w Water saturation, S wir irreducible water saturation, S or residual oil saturation, L z Matrix block hieght.
The performance of the dual-porosity model has been evaluated against a reference solution using several scenarios of various matrix block dimensions. Remarkable differences in oil recovery have been obtained, which highlights the critical limitation of the DP model in emulating the actual flow behaviour in fractured reservoirs, see Fig. 2. The oil recovery of various block dimensions was simulated using the DP model, and the outcome has compared with the oil recovery of fine-grid modelling, which has been utilised as a reference solution. The comparison outcome had referred to an increment in the differences when the matrix block volume increased.
It is worth mentioning that the highlighted differences in Fig. 2 are not constant, and the performance of the DP model is notably affected by the fluid characteristics and the petrophysical properties of the matrix medium. Therefore, a quick check comparison between the DP model and the fine-scale model using the actual properties of the studied reservoir would be necessary to understand the behaviour of the DP model in this specific reservoir.
Simulating the fluid flow behaviour in fractured reservoirs using the dual-porosity model required an in-depth evaluation to address the accompanying issues of the saturation front, recovery speed and the ultimate oil recovery. Therefore, two solution paths were often suggested to solve these problems. First, use a proper transfer function to accurately account for the fluid exchange rate between the matrix medium and fractures. Second, using a suitable discretisation scheme for the matrix block that allows a gradual movement of the invading fluid front towards the matrix block centre, which results in a smooth advance of the fluid saturation.
Defining a proper transfer function is a challenging task. The constant transfer function form of (Kazemi et al. 1976;Warren and Root 1963) and later (Kazemi et al. 1992), as explained in Eq. (4), is not able to address the saturation front movement and match the recovery curve of the reference solution of the fine-scale model.
where (σ) is the shape factor for the matrix block with the volume (V m ), and a fracture-matrix contact surface (A i ) with a distance from the contact surface to the matrix block centre (d i ). Several solutions have been proposed to improve the mathematical representation of the fluid exchange rate, capture the transient effect, and enhance the simulation accuracy of the DP model. The solutions include using a saturation-and pressure-dependent shape factor and multirate transfer function (e.g. (Lu et al. 2006;Lu and Blunt 2007)), a time-dependent shape factor (e.g. (Heinemann and Mittermeir 2016;Su et al. 2013)) and deriving the transfer function based on lab experiments (e.g. (Alkandari 2002)). However, the drawback of the above-mentioned suggestions is the applicability difficulties, stability issues, the matching parameters that are required to be modified frequently, and complicated derivation of the dynamic shape factor.
The second option is dividing the matrix block using the MINC subdomains principle instead of conventional matrix discretisation. The typical discretisation results in The oil recovery performance of the dual-porosity model and the fine-grid model results for various dimensions of equal matrix block dimensions of an oil-water system using the ECLIPSE black oil simulator a significant increment of gridblocks in the simulation model, hence higher CPU (Central Processing Unit) time is required compared with the classical dual-porosity model (e.g. (Gilman and Kazemi 1988)). In contrast, the MINC approach can obtain a better matching to the recovery curve of the fine-grid using a significantly less number of cells. Ahmed Elfeel et al. (2016) have provided an example of using two sub-blocks of the MINC approach with the calculation procedure based on the developed three-phase simulation model of Thomas et al. (1983). Furthermore, the mathematical formulations of their MINC model have been provided in the next section.

The MINC model, a simple case of two subdomains
Mathematical representation of the fluid flow in the matrix subdomains model is similar to the dual-porosity model, but the fluid exchange in the matrix occurs in two stages. The first sub-block represents the outer matrix shell, while the second sub-block represents the inner matrix. The conventional transfer function depicts the exchange rate between the outer matrix shell and fractures. Meanwhile, the second transfer function has been used to represent the interflow between the matrix sub-blocks, as illustrated in Fig. 3. The flowing fluid in this system is represented by Eq. (5): where (T w f ) is the transmissibility in fractures, as defined in Eq. (6) for the flow in the x-direction: Meanwhile, (τ w ) represents the total exchange rate of water between matrix sub-blocks and fractures, as illustrated in Eq. (7): where the individual exchange rate for each sub-blocks can be written in the illustrated forms in Eqs. (8 and 9) to represent the total exchanged rate of water in the double block model; The subscripts (1), and (2) in the above equations refer to the outer and inner sub-blocks of the matrix, respectively. The relative permeability, water pressure and water height of fractures are used as the upstream properties for Eq. (8). Meanwhile, the relative permeability, pressure, and water height of the outer matrix block are used as the upstream properties for Eq. (9). The volume of the matrix sub-blocks is an essential scaling factor, as the calculation of the total water saturation in each time step is scaled back to the total matrix block volume, as illustrated in Eq. (10): The thickness of the outer block matrix controls the subblocks volume, which influences the total water saturation. Ahmed Elfeel et al. (2016) have pointed out that the thickness of the outer matrix block depends on the matrix permeability, fluid viscosities, and dimensionless analysis, which can be estimated based on empirical calculations. Therefore, the outer matrix thickness can be utilised as a matching parameter that can be manipulated to match the recovery curve of the fine-grid model. However, using a constant thickness value (i.e. constant matrix fraction volume) can help to reduce the dependency of the solution on the matrix sub-volumes when the matrix or fluid properties changed.

Methodology and simulation workflow
In this work, we used two subdomains for the matrix block in the MINC model to evaluate the improvement in the simulation accuracy in fractured reservoirs. In the simulated scenarios, the following assumptions were applied: ig. 3 A A simplified case of the MINC model, which comprises two matrix subdomains, illustrating the transfer functions between the matrix sub-blocks and fractures B Inner and outer matrix blocks with their dimensions • It has been assumed that fractures are being filled by water instantaneously (i.e. endless support of water) to emulate the waterflooding process or the advancement of aquifer water. Meanwhile, the saturation of water in both matrix nodes (outer and inner) is calculated at each time step. • An equal volume of both outer and inner matrix blocks has been assumed to reduce the dependency of the simulation results on the matrix sub-block volumes (i.e. stop using matrix sub-volumes as a matching parameter). • Both fractures and matrix sub-blocks are under equilibrium conditions. However, the potential differences are due to the density differences besides capillary pressure forces, which provide the driving mechanisms for the fluid exchange between the matrix blocks and the surrounding fractures.
For two subdomains MINC model, the water saturation in the outer and inner matrix blocks is calculated using Eqs. (11 and 12). In Eq. (11), the fracture properties are used as the upstream flow source for the outer matrix block. Meanwhile, the characteristics of the outer block are the upstream flow source for the inner matrix block of Eq. (12).
In the first part of this work, we propose a correction term to Eq. (11) to adjust the water saturation in the outer matrix block. In the dual-porosity principle, the fluid exchange term is used, which means that the amount of flowing water from fractures to the matrix block is replacing the same amount of oil that flow back to these fractures (i.e. the matrix is a source term in the DP model). Applying the same fluid exchange principle on the inner matrix block of the MINC model results in correcting the water saturation ( S n+1 w1 ) in Eq. (11) to the amount of exchanged water in the previous time step ( ΔS n 2 = S n w2 − S n−1 w2 ) from the inner block. This water saturation amount is equal to the replaced oil from the inner matrix block towards the outer matrix block of the model, and hence, it should be subtracted from ( S n+1 w1 ) in each time step. Therefore, Eq. (11) can be adjusted and re-written to the illustrated form in Eq. (13) to account accurately for the exchanged water in each time step: S n+1 w2 = S n w2 + 0.001127 2 Δt Moreover, the calculation procedure can be summarised in the following steps:a. Provide all the required input data such as fluid densities, viscosities, relative permeability, capillary pressure, and matrix block dimensions, . b. Calculate the shape factors (σ 1 , σ 2 , σ Z1 , and σ Z2 ) for the inner and the outer matrix shells based on an equal volume fraction. c. Calculate the water saturation in the outer matrix layer using the upstream relative permeability of fractures d. Calculate the phase permeability for the outer matrix layer, which is going to be used as the upstream relative permeability for the inner matrix block. e. Calculate the water saturation in the inner matrix block f. Correct the water saturation in the outer shell from the previous step by substracting the term ΔS n 2 , then repeat the steps.
The shape factors (σ 1 , σ 2 , σ Z1 , and σ Z2 ) of both matrix sub-blocks in the above equations can be simply calculated using the illustrated formula in Eq. (4). However, the calculation of the gravity shape factor of the outer matrix (σ Z1 ) block is slightly different. Further calculation details of each shape factor can be found in (Ahmed Elfeel et al. 2016).
On the other hand, the second part of the workflow comprises using matrix heterogeneity in the MINC model, such as; permeability, porosity, initial water saturation, and the exponents of the flow curves (i.e. relative permeability and capillary pressure curves). The flow properties of the matrix are the primary controller parameters on the fluid exchange between the matrix and fractures (Babadagli and Al-Salmi 2004;Maslennikova 2013). The conventional permeability evaluation using Routine Core Analyses (RCA) and well logs were generally employed to populate the matrix block with the permeability. Meanwhile, in the RCA plug sampling process, the permeability of the micro-or nano-scale fractures are often neglected, as the recovered core has already parted in the fractured sections, and it becomes difficult to obtain a valid core plug sample for experimental purposes when such type of fractures exist. Hence, the flow characteristics of the fractured intervals are mostly underestimated in the laboratory assessment. Besides, natural fractures occur at various scales and sizes, and they are not limited to the outer boundary of the matrix blocks, but they extend to the edges of the matrix block itself in the form of discontinued fractures, as illustrated in Fig. 4 (indicated by red circles).
(13) S n+1 w1 = S n w1 + 0.001127 1 Δt In addition to the visible fractures, hairline fractures are also pervasive but difficult to observe in the outcrop or core samples (Nelson 2001). Despite discontinued fractures having low effective permeability compared to the main fracture network, they can increase the contact area between the matrix and fractures and boost the fluid exchange rate. The small-scale fractures at the block edges can significantly enhance the low matrix permeability into several orders. Furthermore, this improvement in the matrix permeability can be modelled in the outer matrix block of the MINC model, as shown in Fig. 5A, B, and C.
The variation of the matrix permeability has often averaged in the upscaling process using various techniques (e.g. arithmetic, geometric and harmonic average) to improve the representation of the permeability and to agree with the interpreted well test permeability. Nevertheless, the performance of the lumped matrix permeability into a single value for the dual-porosity model has a different performance from the real block permeability when the fine-grid model is used in the characterisation. Therefore, a set of the outer matrix block permeability has been used to evaluate the effect of permeability variations on the oil recovery using the classical DP model, two subdomains MINC model, and the finegrid model as a reference solution. This evaluation process is not limited to the matrix permeability, but it can be extended to include porosity, initial water saturation, relative permeability, and capillary pressure. Where the previously mentioned properties may vary based on the lithofacies type or depositional texture in carbonates, as detailed in Dunham classification (Dunham 1962;Lucia 2007).

Validate the suggested adjustments to the MINC model
The provided equations in the workflow have been coded in MATLAB environment using the provided calculation procedures by Thomas et al. (1983) and Ahmed Elfeel et al. (2016). The coded procedure allows using the MINC model as a classical DP model when the fractional volume is equal to zero or one. Then, a validation process was performed to QC the written code by comparing the oil recovery performance between the dual-porosity model using the MATLAB code and the dual-porosity model of the commercial simulator of "ECLIPSE". The comparison referred to an excellent matching with Root Mean Square Error (RMSE) of (0.44%), as shown in Fig. 6.
Next, the MINC model with the suggested adjustments has been tested against the original MINC model, the classical DP model, and the fine-grid model to assess its performance using various matrix block sizes. The results of these scenarios are illustrated in Fig. 7. Moreover, the accuracy of the adjusted MINC model has been evaluated by calculating RMSE using Eq. (14), and the results are listed in Table 1.  Fig. 4. B The discontinuous fractures at the block edges can be modelled in the outer matrix block of the MINC model (km 2 ). Meanwhile, the inner matrix block preserves the matrix flow properties (km 1 ). C A comparison of the oil recovery using several scenarios of permeability increment in the outer matrix block permeability increased using the MINC model

Adaption of matrix heterogeneity in the adjusted MINC model
The heterogeneity of the matrix properties has been evaluated using multiple scenarios for each property. Then, the RMSE of the oil recovery has been calculated between the DP model, the MINC model, and the result of the fine-grid model. This evaluation process helps to understand the impact of matrix heterogeneity on the performance of the DP model and MINC model. Besides, highlighting the goodness of these mathematical models in capturing the various heterogeneity of fractured reservoirs and the expected differences in their oil recoveries. In this part of the work and onward, we use the "MINC model" to refer to the adjusted MINC model for simplicity.
The matrix permeability is the first selected parameter to evaluate its impact on the ultimate oil recovery. Several scenarios of matrix block permeabilities have been used to achieve this purpose. The enhanced permeability of the outer matrix block due to discontinuous fractures can be modelled explicitly in the fine-grid model and the double block model. Meanwhile, an average value should be used for the dual-porosity model where a single value only can be assigned to the matrix cell. Various averaging techniques were tested, and the outcomes of the suggested scenarios were compared to understand the performance of each model at various levels of matrix heterogeneity. Two cases of matrix block dimensions were evaluated to understand the block dimensions effect, besides the contrast in the matrix permeability, as explained below:

Case 1: equal matrix block dimensions
A single matrix block, as shown in Fig. 5(A), was used to construct a single-cell model with equal dimensions of 10ft × 10ft × 10ft. The total number of the matrix block is 1, 2, and 1000 cells for the DP, MINC, and fine-grid models, respectively, as shown in Fig. 8.
The fracture network around the matrix has assumed to be fully saturated with water with endless water support, while the matrix blocks are saturated with oil and irreducible water saturation. A straight line relationship of relative permeability was used for fractures. These conditions were provided to mimic the water advances scenario into an oil zone due to the aquifer potential when the reservoir pressure dropped by production or the water advances due to water injection. Further input data are given in Table 2, which is almost similar to the data used by (Aljuboori et al. 2019): The variation in matrix permeability has been modelled explicitly using the fine-grid model and the MINC model, where the light cyan colour in Fig. 8 represents the inner region (k m1 ), and the light blue colour refers to the edges or the outer region (k m2 ), as illustrated in Fig. 8B and C. In contrast, an average value should be assigned to the DP model indicated by the light cyan area in Fig. 8A. Improved matrix permeability scenarios have been assumed at the edge of the matrix block region (k m2 region) using the permeability values (5, 10, 20, and 40 mD), while the inner matrix region preserves the matrix permeability of 1mD. As a single value should be used in the dual-porosity model, an average permeability was calculated by using the arithmetic, geometric, and harmonic method to examine the best average technique that best matches the oil recovery of the fine-grid model. The results of the above-mentioned scenarios are illustrated in Fig. 9.
The outcome of the MINC model has shown a better agreement with the result of the fine-grid model than the DP model, as illustrated in Fig. 9. Moreover, the DP model has overestimated the oil recovery regardless of the used average techniques for the matrix permeability. Nevertheless, the harmonic average of the matrix permeability provides a closer matching with the reference solution and less error percentage compared with other average techniques, especially when the differences in matrix permeability is significantly increasing. The RMSE percentage for the explained scenarios in Fig. 9 is shown in Table 3.

Case 2: non-equal matrix block dimensions
In nature, the matrix blocks may form at various sizes and shapes (Ahmed Elfeel et al. 2013;); see Fig. 4. The variations of the matrix block sizes are Fig. 6 Comparison of oil recovery between the dual-porosity model of the commercial ECLIPSE simulator and the dual-porosity model using the MATLAB code. The result referred to an identical performance with a negligible RMSE of 0. 44%  Fig. 7 Comparison of oil recovery between the dual-porosity model using the commercial ECLIPSE simulator, fine-grid modelling, original MINC model, and the adjusted MINC model using various matrix block dimensions. The obtained results referred to significant improvements in the simulation accuracy Table 1 Comparison of the illustrated oil recovery performance in Fig. 7   related to the fracture intensity and frequency besides their scale. The non-equal block dimensions have a direct impact on the calculated shape factor and the drainage height, and hence an influence on the fluid exchange rate between the matrix and fractures (cf. Geiger et al. 2013;). Therefore, a non-equal matrix block has been used in this assessment to explore its effect on oil recovery. The matrix block has a dimension of 5ft × 8ft × 12ft in the x-, y-, and z-direction, as illustrated in Fig. 10. The scenarios of the equal matrix block dimensions case have been repeated using the MINC model and the DP model. Then, the results were compared with the reference solution of the fine-grid model. The simulation results of Case 2 are presented in Fig. 11, while the evaluated RMSE is listed in Table 4. Similarly, additional scenarios were implemented using various porosity values and initial water saturation. These conditions were provided to represent additional matrix heterogeneities in the system to emulate a realistic variation of the petrophysical properties in carbonates due to various lithofacies or rock textures. Such heterogeneities can have a significant impact on the oil recovery and fluid saturation curves in the fractured reservoir model. The result of these scenarios is illustrated in Fig. 12.
Moreover, other flow properties (i.e. relative permeability and capillary pressure curves) can also be varied as they represent the rock properties of each lithofacies. Therefore, five sets of capillary pressure and relative permeability have been used to evaluate their impact on oil recovery. The used sets have the same endpoints and capillary pressure value at the irreducible water saturation, where the differences are the curvature path that controlled by Corey exponents (n o , n w , n p ), as shown in Eqs. (15, 16, and 17), (Ahmed 2019). The result of the above-described sensitivity scenarios is illustrated in Fig. 13.

Discussion
Capturing the transient effect is one of the well-known limitations of the dual-porosity model. During the early time of flow, the amount of exchanged oil between the matrix and fractures are overestimated due to the large exposed matrix volume for the fluid exchange. This effect is observed in Fig. 2, where the differences between the DP model and the reference solution decreased when the matrix block volume decreased. Quantitatively, in equal matrix block dimensions of 5ft, the RMSE between the DP and the fine-grid model is about 3%, while the difference increased to 5.5% when block dimensions become 50ft. Furthermore, the saturation of the invaded phase increased rapidly in the matrix block and lead  to an exaggeration of the capillary pressure contribution to the overall oil recovery, which resulted in an overestimated performance. The proposed correction of the water saturation is only applicable to the illustrated outer subdomains model in Fig. 1, except the innermost sub-block. Every outer matrix subdomains transfer the water towards the innermost cells and transfer the oil back towards the surrounding fractures. However, this principle is not valid for the match stick model or layered model as there is no innermost block, and all matrix sub-blocks are in contact with fractures and can exchange fluid directly.
MINC Model or any discretised form of the matrix block (e.g. the sub-face approach (Abushaikha and Gosselin 2009), match stick model, and layered model (Kazemi and Gilman 1993)) are mainly aimed to slow down the fluid exchange speed. This objective can be achieved by reducing the exposed matrix volume to the fluid exchange, and hence a lower oil recovery in the early time to improve the match with the reference solution. Adaption of any form of discretisation can enhance the matrix characterisation by including additional heterogeneities in the simulation model, which increases the simulation accuracy. Moreover, it can accurately capture the saturation value of the water phase in the matrix medium.
The heterogeneity in fractured reservoirs, such as matrix block shapes, volumes, and intrinsic properties, plays a significant role in controlling the recovery differences between the DP model and the fine-grid model. Lumping the inherent properties (i.e. permeability, porosity, initial water saturation, capillary pressure curves, and relative permeability curves) into one value or set is Fig. 9 Comparison between the oil recovery of the adjusted MINC model, DP model with an average matrix permeability (k m1 and k m2 ) using arithmetic, geometric and harmonic methods, in addition to the result of fine-grid modelling as a reference oversimplifying the actual matrix properties. In such a case, the model might not capture the actual potential of the fluid flow in the matrix blocks. Therefore, upscaling the matrix properties in the fractured reservoirs are most often failed to address the high variability in the fractured medium, hence unreliable recovery prediction.
The frequent attempts to modify the transfer function often had limited success to match the outcome of the finegrid model due to maintaining the same DP principle of using the lumped matrix properties. Therefore, introducing any level of heterogeneity to such models can lead to a significant deviation from the accurate flow performance, and hence unable to represent the actual behaviour of the matrix medium. An exclusion can be made for the multi-rate model (Di Donato et al. 2007), as they have considered the sub-grid heterogeneity in their calculations. Multiple Interacting Continua can provide sufficient space to capture the matrix heterogeneity. The tested scenarios have shown the advantage of the MINC model by using a minimum number of matrix subdomains (i.e. two subdomains only), which resulted in significant improvement in the simulation accuracy with a maximum RMSE of 1.73%. However, the maximum RMSE in the counterpart DP model is 12.5% when the arithmetic average permeability is used. The used case of MINC in this work, which consists of two matrix sub-blocks, has provided an excellent matching with the fine-grid model results; see the RMSE in Tables 1, 3, and 4. Furthermore, additional matrix heterogeneities can be adapted in the model and represented explicitly, such as porosity, water saturation, relative permeability, and capillary pressure, which increase the accuracy of the modelling process and obtain a reliable prediction.

Conclusion
1. The dual-porosity model using the Gilman and Kazemi model (Gilman and Kazemi 1988) tends to overestimate the ultimate hydrocarbon recovery due to overestimating the sweep efficiency (Di Donato et al. 2007). Moreover, the calculated error had increased when the matrix heterogeneity was introduced to the model, such as permeability. The calculated RMSE of the DP model increased up to 12.5%, while the maximum RMSE of the MINC model in the investigated scenarios is 1.73%. 2. Improper use of the DP can hide the actual improvements of various development scenarios, as the propagated error of the DP model might exceed the obtained improvement by a proposed scenario to improve oil recovery (e.g. (Aljuboori et al. 2020a)). 3. Although the MINC model required additional computational resources compared with the classical dualporosity model, excellent simulation accuracy can obtain with marginal RMSE compared with the outcome of fine-grid modelling with significantly less CPU time, see the comparison of the CPU simulation time in (Yan et al. 2016). 4. The calculated error was raised in the non-equal block dimensions when the dual-porosity model has used compared with the MINC model. Nevertheless, the harmonic average provides the best average technique compared with the arithmetic or the geometric average, as the harmonic method shows lower RMSE percentage, which gradually decreased when the outer matrix block permeability increased, see Table 4. 5. The matrix heterogeneity can be modelled explicitly in the MINC model, and their impact has been excellently captured as illustrated by the oil recovery curves. Moreover, the higher level of heterogeneity may require an additional partitioning layer of the model to capture  Yan et al. 2016)). 6. The matrix heterogeneity plays a significant role in various development scenarios. In low salinity waterflooding, for example, the alteration to the relative permeability and capillary pressure curves control the ultimate oil recovery and the process performance. Therefore, including additional petrophysical properties (e.g. initial water saturation, multiple sets of relative permeability and capillary pressure, see Fig. 13) in the simulation model can enhance the matrix characterisation, and hence increase the simulation accuracy of low salinity waterflooding (e.g. (Aljuboori et al. 2020a(Aljuboori et al. , 2020b) Fig. 13 Sensitivity scenarios of the capillary pressure exponent, water, and oil relative permeabilities exponents of the outer matrix block. The exponent of the oil relative permeability exhibited a higher impact on the ultimate oil recovery than both exponents of capillary pressure and relative permeability of the water phase