Regular Article
The MINC proximity function for fractured reservoirs flow modeling with nonuniform block distribution
Lebanese American University, P.O. Box 36, Blat, Byblos, Lebanon
^{*} Corresponding author: nicolas.farah@lau.edu.lb
Received:
1
February
2020
Accepted:
15
December
2020
Reservoir simulation is a powerful technique to predict the amount of produced hydrocarbon. After a solid representation of the natural fracture geometry, an accurate simulation model and a physical reservoir model that account for different flow regimes should be developed. Many models based on dualcontinuum approaches presented in the literature rely on the PseudoSteadyState (PSS) assumption to model the interporosity flow. Due to the low permeability in such reservoirs, the transient period could reach several years. Thus, the PSS assumption becomes unjustified. The numerical solution adopted by the Multiple INteracting Continua (MINC) method was able to simulate the transient effects previously overlooked by dualcontinuum approaches. However, its accuracy drops with increasing fracture network complexity. A special treatment of the MINC method, i.e., the MINC Proximity Function (MINC–PF) was introduced to address the latter problem. And yet, the MINC–PF suffers a limitation that arises from the existence of several gridblocks within a studied cell. In this work, this limitation is discussed and two possible solutions (transmissibility recalculation/adjusting the Proximity Function by accounting for nearby fractures) are put forward. Both proposed methods have demonstrated their applicability and effectiveness once compared to a reference solution.
© N. Farah & A. Ghadboun, published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Oil and gas production from unconventional resources (shale gas, tight oil) is, for the most part, uneconomical due to the lowpermeability and highheterogeneity reservoirs. These are characterized by high fracture irregularity which renders flow simulation to be very challenging. To overcome this problem, fluid flow simulation can be simplified by sacrificing fracture network complexity. Another challenge is the very low permeability in the tight matrix as compared to fractures which makes the main flow mechanisms (excluding gravity) primarily diffusive – from matrix to fracture and vice versa.
The literature has regarded modeling fractures in shale gas formations carefully (see Cipolla et al., 2010; Freeman et al., 2010; Li et al., 2011; Moridis et al., 2010; Rubin, 2010; Wu et al., 2012). However, few of the proposed approaches address the critical particularity of flow from unconventional oil/gas reservoirs with complex fracture networks (irregular fracture distribution). On real field scales, most of these models proved inadequate (Raterman et al., 2019). Additionally, most of the studies conducted on unconventional reservoirs used inhouse reservoir simulators that were initially developed for conventional reservoirs. This led to a limited understanding of unconventional reservoirs where flow from multiscaled fractured reservoirs remains a challenge to modeling approaches.
Earlier models such as singleporosity, dualcontinuum (Barenblatt et al., 1960; Chang, 1993; Coats, 1989; Kazemi et al., 1976; Lim and Aziz, 1995; Noetinger, 2015; Quintard and Whitaker, 1996; Van Heel and Boerrigter, 2006; Warren and Root, 1963) and multicontinuum models were all used for flow modeling in unconventional reservoirs.
An explicit discretized model with a singleporosity approach can accurately model regular fracture networks because fractures are identified by their spatial distributions (a detailed description of the governing equations is presented in Appendix A). As multiscale fractures become abundant in unconventional reservoirs post stimulation, single porosity models that use an LGR (Local Grid Refinement) technique (see Cheng, 2012; Cipolla et al., 2009; Ding et al., 2014) become CPUtime consuming due to the large number of grid cells used.
On the other hand, dualcontinuum models proved a worthy tool when modeling conventional fractured reservoirs (short transient periods) with high fracture density (see Kazemi, 1969; Wu et al., 2004). However, studies showed that modeling gas flow from low permeability fractured shale formations lies beyond the capability of these models (Ding et al., 2014; Farah et al., 2015). The long transient period caused by large matrix block size and low matrix permeability is not accounted for. Thus, the condition of a PSS flow cannot be assumed (Wu and Pruess, 1988). This issue was overcome by the Multiple INteracting Continua (MINC) model.
The MINC was proposed by Pruess and Narasimhan (1982) to model flow in porous media by subdividing the matrix grid block into nested volumes using analytical expressions (see Fig. 1a). The flexibility and ease of integration of the MINC model render it an effective choice when designing flow simulating models (Jiang et al., 2014; KarimiFard et al., 2006; Li et al., 2015; Tatomir et al., 2011). Moreover, in cases of multiphase (gas/water or oil/gas) flow simulations, a challenging task is to predict the GOR (GasOil Ratio). Due to the large grid blocks used by dualcontinuum models, such phenomena cannot be simulated. However, this problem could be accurately solved using the MINC method (see Delorme et al., 2016; Farah et al., 2018; Ricois et al., 2016). Recently, Farah and Delorme (2020) proposed a Unified Fracture Network Model (UFNM) for unconventional reservoirs modeling. Thus, the UFNM’s ability mimicking flow from irregular fracture networks while considering complex physical phenomena was analyzed and evaluated. Besides, the simulation of matrixfracture exchanges was highlighted.
Fig. 1 Schematic of the MINC concept for (a) a regular fractures network and (b) for an arbitrary/irregular fractures distribution; (after; Farah et al. (2018)). 
On the other hand, the MINC method becomes challenging in the presence of irregular fracture networks (see Fig. 1b). To overcome this problem, the MINC Proximity Function was proposed. The MINC–PF implements a fully numerical method that accounts for irregularity faced in unconventional reservoirs. This work is the development of the concept proposed in Farah (2016) Ph.D. work, in Delorme et al. (2016) and in Khvoenkova and Delorme (2011).
2 Framing the problem of nonuniformity
Figures 2a and 2b illustrate a regular and an irregular fractures distribution, respectively, where uniform subvolumes (Sv_{1}, Sv_{2}, Sv_{3} and Sv_{4} are equal) are present. The MINC–PF (see Appendix B: MINCPF Generation) was originally proposed by Pruess and Karasaki (1982). It was able to model fluid flow from regular/irregular fracture networks only for identical or similar blocks on each side of the fracture (Figs. 2a and 2b). In fact, based on the fracture network presented in Figures 2a and 2b, quarter of the total volume is found in each subvolume (Sv_{1}, Sv_{2}, Sv_{3} and Sv_{4}) present within the studied grid cell. Moreover, the distribution functions from each subvolume overlap as presented in Figure 2c. In fact, Figure 2c represents the distribution functions from each subvolume of the case presented in Figure 2b. In other words, the sum of all four distribution functions (overlapping) will add up to the total distribution function (dotted blue curve in Fig. 2c).
Fig. 2 Illustration of (a) a regular fracture distribution with uniform block distribution, (b) an irregular fractures distribution with uniform block distribution, and (c) the Proximity Function distribution regarding the diagonal case. 
When handling realistic cases, nonuniform subvolumes/blocks are present within a grid cell (see Figs. 3a and 3c), where a conforming unstructured matrix mesh is needed for an accurate representation of the fracture system (see Figs. 3b and 3d). In such cases, the MINCPF miscalculates the flow. Farah et al. (2018) proposed a DFM based on a MINC Proximity Function for unconventional reservoirs to simulate fluid flow in irregular fracture network. However, the presence of different subvolumes/blocks was not considered. This led to a miscalculation when modeling interporosity flow, thus causing a difference between the DFM results and the reference solution.
Fig. 3 Illustration of (a) a regular fracture network where different subvolumes (subblocks) are present (Sv1a, Sv2a, Sv3a, and Sv4a), (b) the MINC6–PF model, (c) an irregular fracture network distribution occurring in a matrix grid cell where different subvolumes (subblocks) are present (Sv_{1b}, Sv_{2b}, Sv_{3b}, and Sv_{4b}) and (d) its corresponding MINC6–PF model. 
Due to the nonuniformity of fractures presented in Figures 3a and 3c, nonuniform subvolumes (Sv_{1a/b}, Sv_{2a/b}, Sv_{3a/b} and Sv_{4a/b}) occur within the studied grid cell. In such case, different matrix refinements will occur within each subvolume since computing the MINC–PF depends on the distance to the nearest fracture. In other words, some subvolumes will exhibit a higher matrix volume and a larger distance to the fractures than in the others (highlighted in red in Figs. 3b and 3d). Thus, computing a MINC6–PF into the grid cell presented in Figures 3b and 3d leads to a presence of a 6^{th} subdivision in Sv_{4a}, Sv_{1b}, and Sv_{4b} only. In fact, when the contribution of a matrix subvolume ends in the total cumulative distribution curve, no more refinement can be done in this subvolume. So, a 6th subdivision could only be found in Sv_{4a}, Sv_{1b}, and Sv_{4b}.
Therefore, a connection between subdivision #5 and #6 should exist only in these subvolumes. Unfortunately, in the present MINC–PF model implemented in Farah et al. (2018) the 6th subdivision from these subvolumes exchanges flow with all subdivisions #5 from the other subvolumes (Sv_{1a}, Sv_{2a}, Sv_{3a} and Sv_{2b}, Sv_{3b}). Hence, these connections overestimate the exchange area between different subdivisions leading to nonaccurate flow modeling results when compared to a reference solution. So, in case of nonuniform subblocks equation (B.6) cannot be computed this way. Note that, this problem was first addressed by Farah (2016) Ph.D. work.
To solve this problem, three solutions could be envisaged:

Dividing each subvolume into 6 subdivisions.

Recalculating transmissibility values within the studied grid cell taking the assumption of noflow boundary at its borders. This solution is a weighted discretization of the matrix with adapted transmissibility when nonuniform block sizes exist within the studied grid cell.

Adjusting the proximity function taking into consideration the presence of nearby fractures in neighboring cells.
It should be mentioned that implementing one of these solutions to the DFM presented in Farah et al. (2018) will ensure a better flow modeling.
3 Solving the problem of nonuniformity
The first solution will be disregarded because it increases the simulation CPU time. Note that, computing a MINC6PF on the case presented in Figures 3a will result in 12 nodes (6 nodes for Sv_{1a}, Sv_{2a}, Sv_{3a} and 6 nodes Sv_{4a} given the difference of the pressure profiles) instead of 6 nodes. In fact, pressure profiles from Sv_{1a}, Sv_{2a} and Sv_{3a} are the same based on the distance to the fractures. So, nodes from subvolumes Sv_{1a}, Sv_{2a} and Sv_{3a} could be lumped. However, Sv_{4a} exhibits a higher matrix volume and a larger distance to the fractures compared to the others, resulting in additional new 6 nodes.
So, the other two solutions (2 and 3) suggested in the previous part will be studied to improve the MINC Proximity Function. Hereafter, the transmissibility correction, and the adjustment of the Proximity Function are presented on both a regular and an irregular case.
3.1 Transmissibility recalculation – nonuniform grid blocks
3.1.1 Method
For a steadystate flow, the total transmissibility through the subdivisions is described as follows:
(1)where N represents the number of matrix subdivisions and T _{ i,i+1} the transmissibility between #i and #i + 1 from different subblocks existing within the studied grid cell.
Moreover, the total flow between two continuums (matrixfracture or matrixmatrix) is equal to the sum of the flow from these continuums considering all existing subvolumes in a grid cell. So, the transmissibility between continuum #i and #i + 1 considering all the subvolumes are described as follows:
(2)where NV represents the number of subvolumes of the studied grid cell respectively and represents the transmissibility between #i and #i + 1 from different subblocks.
To overcome the earlier discussed error, the connection must be adjusted by correcting the transmissibility from the last subdivision to the fracture under a steadystate regime. The correction could be expressed as follows:

f represents the subdivision at which the transmissibility between #f and #f + 1 needs to be fixed.

represents the corrected/new transmissibility between subdivisions f and f + 1.

represents the transmissibility from the volumes requiring correction/adjustment indexed by CSV (Corrected SubVolume).
The presented work takes into account a regular/irregular fracture distribution. The identification and characterization of such irregularity is also a challenging task if a nonautomated general procedure is applied. In both cases presented earlier in Figures 3a and 3c, regular and irregular fracture networks, the 6th subdivision from subvolumes Sv_{4a}, Sv_{1b}, and Sv_{4b} exchanges flow with all subdivisions #5 from different subvolumes Sv_{1a}, Sv_{2a}, Sv_{3a} and Sv_{2b}, Sv_{3b}. This is not the case and would overestimate the gas production as shown in the following sections. The adjusted transmissibility for both cases could be obtained from the general equation (3).
3.1.2 Application
3.1.2.1 Regular fractures distribution
First, an example of two regular nonsymmetric orthogonal fractures is presented. The following example in Figure 4a consists of a matrix block of 164 ft in x and y directions with the presence of two hydraulic fractures (dashed lines). The fracture aperture and permeability are set at 0.004 ft and 2500 mD respectively. The matrix permeability is 0.0001 mD. This domain contains four matrix blocks with different sizes. The thickness of the block is 20 ft in z direction. The reference solution consists in an LGR technique (Fig. 4a) where fractures are explicitly discretized. On the other hand, the whole domain was discretized using the MINC Proximity Function model. An illustration of the MINC6–PF model is presented in Figure 4b.
Fig. 4 Illustration of (a) LGR reference solution, (b) the MINC6PF model, (c) the total cumulative matrix volume and per subvolume from Sv_{1a}, Sv_{2a}, Sv_{3a} and Sv_{4a}. 
Figure 4c presents the cumulative matrix volume of the distribution function for the total matrix grid cell and per subvolume function of the distance to fractures. As the fracture distribution is not symmetric in the matrix grid cell, the partitioning of the matrix subdivision is not the same around the fractures in each subvolume (Sv_{1a}, Sv_{2a}, Sv_{3a} and Sv_{4a}). However, the maximum distance found in Sv_{1a} is the same as that in Sv_{2a} and Sv_{3a}. A larger maximum distance to the fractures is found within Sv_{4a} due to its geometric form. Based on the distribution function (Fig. 4c) generated for this case, only a sixth matrix subdivision (continuum) could exist in Sv_{4a} due to its different block size inside the grid cell and the bigger volume found within this subvolume (red square as shown in Fig. 4b).
The corrected transmissibility must be calculated considering subdivision #5 and #6 (red square presented in Fig. 4b) only from Sv_{4a} in order to properly model the flow exchange between these two continua. An adjustment on the transmissibility calculation must take into consideration the right exchange surface between these two continua (#5 and #6 only from Sv_{4a}).
In our case, to find the adjusted transmissibility, equation (3) is solved for CSV = 4, f = 5 and N = 6. Hence, equations (4) and (5) are obtained:
The corrected transmissibility, being the only unknown value, could be obtained from equation (5).
3.1.2.2 Irregular fractures distribution
We consider a block of 164 ft in x and y directions containing three irregular fractures as shown in Figure 5a. The net thickness of the reservoir in z direction is 20 ft. In this study, the permeabilities are set at 10^{−4} mD in the matrix and 3 D in the fractures. The fracture aperture is fixed at 0.004 ft. The porosity of the matrix media is 0.05. Only a single phase (gas only) flow simulation is taken into consideration. The initial reservoir pressure is 3800 psi. The bottom hole well flowing pressure is 1000 psi. A reference solution with very fine matrix grid cells is used here (see the approach presented in Appendix C). The reference model consists of 10^{6} grid cells (see Fig. 5a) after discretizing the domain into 1000 grid cells in each x and y directions. An illustration of the MINC6–PF is presented in Figure 5b.
Fig. 5 Illustration of (a) an irregular fracture distribution where the reference solution is performed, (b) the MINC6PF model, (c) the total cumulative matrix volume and per subvolume from Sv_{1b}, Sv_{2b}, Sv_{3b} and Sv_{4b}. 
Figure 5c presents the cumulative matrix volume distribution function for the total matrix grid cell (red solid line) and per subvolume function of the distance to fractures. Based on this information, the maximum distances to the fractures found in each volume are 14 ft, 8 ft, 10 ft, and 38 ft; respectively for Sv_{1b}, Sv_{2b}, Sv_{3b} and Sv_{4b}. A larger distance is found from subvolumes Sv_{1b} and Sv_{4b} due to the fracture distribution and to the bigger volume presented by these subvolumes (Sv_{1b} and Sv_{4b}). As we applied a MINC6 model onto this grid cell, a sixth subdivision is only found in Sv_{1b} and Sv_{4b.} Unfortunately, in the present model the subdivisions #6 from subvolumes Sv_{1b} and Sv_{4b} exchange flow with all subdivisions #5 from different subvolumes (Sv_{1b}, Sv_{2b}, Sv_{3b} and Sv_{4b}).
To find the new/adjusted transmissibilityequation (3) is solved for f = 5 and N = 6. However, in this case the subvolumes needing adjustment are CSV = 1 and 4. Thus, equations (6) and (7) are obtained.
The corrected transmissibility can be obtained from equation (7) being the only unknown value. The transmissibility is calculated with all four subvolumes when the distance is smaller than 10 ft. (interaction length), but has to be computed with only subvolumes Sv_{1b} and Sv_{4b} when the distance is larger than 10 ft. No further refinement can be done into a matrix subvolume when it no longer contributes to the total cumulative distribution curve. However, when a distance is larger than 10 ft. further subdivisions could be introduced (see Fig. 5c).
3.1.3 Results and discussion
Hereunder, Figure 6 illustrates the results (cumulative gas production) for both cases presented in the previous part. Figure 6a presents a comparison of cumulative gas production for the case presented in Figure 4a with and without correcting the transmissibility. The MINC6–PF overestimates the gas production as compared to the reference solution. After the correction, the MINC6–PF can predict the gas production once compared to the explicit discretized model set as a reference solution.
Fig. 6 Illustration of the comparison of the cumulative gas production for (a) the explicit discretized model, the MINC6 Proximity Function with and without correction for 1000 days of production for the regular fractures distribution presented in Figures 4a and 4b and (b) the comparison of the cumulative gas production before and after the transmissibility correction for the irregular fractures distribution presented in Figure 5a. 
As for the irregular case, MINC–PF is performed using a MINC6 and a MINC8 model. The cumulative gas production for different simulation models is shown in Figure 6b for 5000 days of production. Based on Figure 6b, the MINC6–PF and MINC8–PF overestimate the gas production comparing to the reference solution. Note that, the MINC8 model amplifies the error as we expected due to the presence of higher matrix refinements which in turn let the solution diverge. In other words, more problematic surface exchanges occur when using a MINC8PF model. Clearly, as explained in the previous section when the fracture distribution is not symmetric in the studied matrix grid cell, the partitioning of the matrix subdivision is not the same around the fractures in each subvolume (Sv_{1b}, Sv_{2b}, Sv_{3b} and Sv_{4b}). Hence, the transmissibility from Sv_{1b} and Sv_{4b} must be corrected, where the new transmissibility value connecting subdivisions #5 and #6 is denoted presented by equation (7).
After transmissibility correction, the MINC–PF can predict and match the reference solution for the cumulative gas production when it is computed for a MINC6 and a MINC8 model (see Fig. 6b). Moreover, the MINCPF decreases the CPU time from 5 h for the reference solution to only 8 s.
At this early stage of development, the recalculation of the transmissibilities within the blocks needing correction was done manually. However, in our future work, and to test largescale cases and more complex fracture networks, these corrections must be detected and corrected automatically via our preprocessor. This would be possible after we integrate the recalculation method into the DFMMINC Proximity Function developed in Farah et al. (2018).
3.2 Adjusting the proximity function
This solution is characterized by simulating the matrixfracture flow exchange considering the presence of nearby fractures (outside the studied grid cell). In this part, this solution is applied on a nonuniform blocks distribution for regular fracture network and a singlephase flow (gas).
The example presented in Figure 7a relies on a regular fracture network. The fractures are represented by dashed blue lines. This example consists of a fracture spacing of 164 ft in x and y directions and the depth is 20 ft in z direction. The fracture aperture and permeability are set at 0.004 ft and 2500 mD, respectively. The matrix porosity and permeability are 0.05 and 0.0001 mD, respectively. This periodic infinite fracture network can be simulated with the mesh definition presented in Figure 7a by either scenario (1) or (2). However, the fractures are centered in scenario (1) while they are shifted in scenario (2) based on the mesh definition. Note that using an explicit discretized model (where the fracture network is explicitly discretized), the gas prediction from scenarios (1) and (2) should be the same for both models. An explicit discretized model and the MINC–PF are applied on each of the two presented structures. Figure 7b compares the cumulative gas production for the explicit model and the MINC–PF for scenarios (1) and (2).
Fig. 7 Illustration of (a) an infinite fracture network consisting in a regular distribution with a 164 ft of fractures spacing’s where; (1) represents a centered fractures scenario and (2) represents a shifted fracture network scenario compared to the mesh definition, (b) the comparison of the cumulative gas production for the explicit and the MINC proximity function model for scenario (1) and (2). 
As expected, based on Figure 7b, the cumulative gas productions of the explicit model (1) and (2) are the same. Explicit discretized models for scenarios (1) and (2) are set as the reference solutions. Besides, computing the MINC–PF based on a MINC6 method for scenario (1) can match the reference solution as the fractures are not shifted. However, computing the MINC6–PF on scenario (2) underestimates the cumulative gas volume after 5000 days of production.
This underestimation could be adjusted either by correcting the transmissibility for the needed subvolume as presented in Section 3.1 or by adjusting the proximity function taking into consideration nearby fractures. In the present approach only the presence of the fractures inside the studied grid cell for matrixfracture exchange is considered. In other words, this approach treats each grid cell as if a zero flux exists on the boundary i.e., there is no exchange between matrixes from different grid cells. No flow exchange with nearby fractures, no matter how close to the matrix cell, is present. This causes inaccuracy in flow simulations. A correct approach would account for the presence of nearby fractures around the studied grid cell (see Fig. 8a) to solve this problem.
Fig. 8 Illustration of (a) the MINC proximity function computed into the studied grid cell and (b) the cumulative matrix volume per subvolume for model (2) and (2)–NF and (c) comparison of the cumulative gas production for case (1), (2) and (2)–NF. 
In Figure 8a, the step preceding the application of the MINC–PF is computed in the studied cell, where 100 random points (green dots) are launched to plot the distribution function. The current MINCPF approach considered only the dashed red fractures inside the studied grid cell miscalculating the proximity function curves. To correct this distribution function, nearby fractures (green dashed lines) should also be taken into consideration if the distance from sample points (green dots) to nearby fractures (denoted by “d ptNF”) is smaller than the distance to the fractures inside the studied grid cell (denoted by “d ptfrac Cell”) as shown in Fig. 8a). Moreover, Figure 8b shows the cumulative distribution function for the case presented in Figure 8a with and without considering nearby fractures, where the notation “NF” refers to the consideration of nearby fractures.
Clearly, based on Figure 8b, the point distributions from Sv_{4a} as well as the total distribution are affected by the presence of nearby fractures. In fact, the sum of the four distribution curves resulting from Sv_{1a}, Sv_{2a}, Sv_{3a} and Sv_{4a} is the total matrix volume distribution function. Based on the results from Figure 8b, the maximum distances from a random sample point (green dot in Fig. 8a) to fractures with and without considering nearby fractures are 82 ft and 145 ft, respectively. Note that, the grid cell size is 164 ft.
Figure 8c presents the comparison of the cumulative gas production for the explicit discretized model for scenarios (1) and (2) and the MINC6–PF for scenario (2) with and without taking into consideration nearby fractures. Clearly, based on results from Figure 8c, considering the presence of nearby fractures our MINCPF accurately simulates matrixfracture exchange as it conserves the matrix block size. Hence, it adjusts the cumulative gas production for MINC–PF scenario “(2)–Nearby Fractures” to match the reference solution (Explicit Discretized scenario (2)).
Finally, the problem of the presence of nonuniform grid blocks (that leads to a miscalculation when modeling interporosity flow using the MINCPF model) is solved using one of the methods presented above. The solution consists in either transmissibility recalculation within the needed subvolumes or by adjusting the proximity function distribution taking neighboring fractures into consideration. Clearly, based on our results from Sections 3.1.3 and 3.2 both proposed methods have demonstrated their applicability and effectiveness once compared to a reliable reference solution. On the other hand, when either method is implemented additional cost in terms of preprocessor computing time is needed. This matter will be addressed within our future work once either method is performed on large/complex cases and when an automated procedure (concerning transmissibility recalculcation or proximity function adjustment) is implemented. Important to note is that the reservoir simulation CPU time is not affected when considering either solution presented above.
4 Conclusion
Unconventional reservoirs present multiscale fractures embedded in lowpermeability matrix. Neither an explicit discretized model using a singleporosity with an LGR technique nor a dualcontinuum model is suitable for practical flow modeling from unconventional reservoirs simulation. The first is too much CPU time consuming due to the large amount of grid cells needed to conform fractures distribution, while the second is not accurate due to the PSS regime assumption in matrixfracture flow exchange. On the other hand, due to the very low reservoir permeability faced in unconventional formations, the transient period could last several years before the PSS regime is reached. The MINC approach showed its ability to handle transient flow modeling from unconventional reservoirs while achieving considerable accuracy and reducing the CPU time compared to the reference solution (Farah et al. (2015). The MINC method can simulate matrixfracture interactions from unconventional reservoirs for different flow regimes (Delorme et al., 2016; Farah, 2016; Farah et al., 2015) and can be easily implemented in any Discrete Fracture Model (Farah et al., 2018).
However, the MINC method becomes challenging in the presence of irregular fracture networks. To this matter, the MINC–PF presents a solution to this problem and is efficient when uniform subvolumes/blocks within the studied grid cell are faced. However, due to the complex fracture network distribution and the presence of nonuniform subvolumes/blocks, the MINC–PF falls short of accuracy by miscalculating the flow.
The novelty of this work lies in the numerical method of the MINC–PF and its two possible improvements presented and discussed in Sections 3.1 and 3.2.
The first relies on recalculating the transmissibility within the studied grid cell. This weighted discretization was performed manually for this paper thus allowing for simple case testing only.
The second method consists in adjusting the proximity function in the studied grid cell to account for nearby fractures from neighboring cells.
In future work, the integration of either method into MINCPF will increase the accuracy of the DFMMINC Proximity Function developed in Farah et al (2018). Our future work will consist in testing more complex scenarios (largescale cases) while exploring both methods after introducing them into the DFMMINC Proximity function.
Appendix A
Governing equations
The general equation governing threephase, multicomponent, threedimensional flow for a single porosity model is expressed using the following equation:
(A.1)where, the subscript p represents the phase, ∅ corresponds to the matrix or fracture porosity, C _{cp} is the mass fraction of component c in phase p, S _{p} the saturation of component c in phase p, ρ _{p} is the mole density of phase p, v _{sg} corresponds to the gas sorption term, is the velocity of phase p and is the molecular diffusion and dispersion flux of component k in phase p. Also, q _{p} is the sink/source term of phase p per unit volume of formation.
The phase velocity in both media is expressed in both media using the following equation:
(A.2)where, Z is depth (positive, increasing downwards), g is the algebraic value of gravitational acceleration projection on Z axis, k is the absolute permeability tensor of the medium, P _{p} the pressure of phase p, μ _{p} the viscosity of phase p, k _{rp} the relative permeability of phase p.
The sorption terms in equation (A.1) is calculated using:
(A.3)where, V _{s} is the volume of adsorbed gas in standard conditions per unit mass of solid, ρ _{r} is the solid rock density, and ρ _{sc} is the gas mole density at standard condition.
We further reduce equation (A.1) into:
(A.4)Equation (A.4) is discretized in space using a controlvolume method, where time discretization is carried out using a backward, first order, fully implicit, finitedifference scheme. Explicit spatial discretization of the flow terms (transmissibility) is carried out. By neglecting the sorption terms, equation (A.4) becomes:
(A.5)where, the superscript n denotes the previous time step, and n + 1 the current time step to be solved; ∆t ^{ n+1} is the timestep size; V _{ i } is the volume of the cell _{ i }; N _{ i } contains the set of direct neighboring cells j of cell n; F _{ p,ij } is the flow term between cells i and j; and q _{ p,i } is the sink/source term in cell i.
When Darcy’s law is applied, flow terms in equation (A.5) are rewritten as:
(A.6)where, Г_{ ij } is the interface between cells i and j, λ _{ p } is the mobility term of phase p, k is the absolute permeability, Φ is the potential, and n is the normal direction at the interface Г_{ ij }. Considering a two point flux approximation scheme, equation (A.6) becomes:
(A.7)where, λ _{ I,ij } is calculated with an upstream scheme; T _{ ij } is the transmissibility between cells i and j.
Using the MINC method, the matrix is subdivided based on the distance from the fracture in order to build a grid of nested meshes. The flow regime is described by equation (A.4). The discretization of the flow is described by equation (A.5) with the flow exchange between two neighboring subcells described by equation (A.6). The approximation of the flow term F _{ ij } is described in equation (A.7) where the transmissibility T _{ ij } is given by:
(A.8)where, A _{ ij } is the area of the interface between these two subcells and D _{ ij } is the average distance between the two sub grids.
Appendix B
MINCPF generation
To construct the N subdivisions of MINC, we divide the considered grid cell G according to the distance to the fractures inside this grid cell. Let G(x, y, z) be the distance of a point (x, y, z) to the fractures. The domain with the distance to the fractures smaller than a given r is defined by:
(B.1)and its corresponding volume can be calculated by:
The N subdivisions of MINC can therefore be defined as follows:
(B.3)where r _{1} ≤ r _{2} ≤…≤ r _{ N } .
The corresponding volumes are given by:
Usually, the subdivisions are parametrized by a given percentage of the total volume, defined by the user. So, we can consider the volumes known. Unknowns are the distances which separate two consecutives subdivisions.
We notice that the volume of a subdivision can also be calculated in another way. Let A(r) be the area of the surface (inside the grid cell G) with a distance r to the fractures. The volume of the domain D with a distance smaller than r to the fractures can be computed by:
So, if the volume of the subdivision at a distance r to the fractures is known, the area A of its surface is just the derivative of the volume with respect to the distance:
(B.6)where, A _{ i+1} corresponds to exchange area between two consecutive volumes (i, i + 1), the distance D _{ i+1} corresponds to the average distance from volume V _{ i+1} to the fractures and D _{ i−1} is the average distance of V _{ i−1}. The volumes V _{ i−1} and V _{ i+1} correspond to the cumulative volume of grid cells i + 1 and i – 1. respectively.
Once the exchange area is known, the connection transmissibility between two consecutive matrixes subdivisions is calculated by equation (B.7), where k _{ m } is the matrix permeability:
(B.7)where, D _{ i,i+1} corresponds to the average distance between two successive matrix subdivisions and N corresponds to the number of matrix refinement.
In some cases, the volume of a subdivision V(r) can be expressed analytically, so the interface areas are obtained by an analytical expression (see Hajibeygi et al., 2011). However, in most cases, analytical expressions cannot be found, and numerical methods should be used instead to compute the subdivision volumes and their interface areas (Farah et al., 2018).
In most cases, we cannot find an analytical expression for the volume computation, and we can determine neither the distance r to separate two subdivisions of a MINC. A numerical solution is required. The grid cell G is first discretized in p subdomains of equal volumes ( with and for i ≠ j. Any integral over a domain inside G can be computed as the sum of the integral over the domain G _{ i }. Therefore, the volume of MINC subdivisions D _{ k } (k = 1, 2, …, N) can be calculated by:
On the other hand, we have:
For a very small domain G _{ i }, we assume that it is either inside or outside the domain D _{ k } by choosing a point pt _{ i } (the center of gravity or and randomly selected point) of G _{ i } and checking if it belongs to D _{ k }. This assumption is reasonable if the discretization is fine enough, that is, for p → ∞ and vol(G _{ i }) → 0. So, the volume of D _{ k } is approximated by:
(B.10)where, p _{ k } is the number of discretized subdomain G _{ i } inside the MINC subdivision D _{ k }.
In equation (B.10), approximation errors are committed around the boundary of D _{ k }, where a discretized subdomain G _{ i } is neither inside nor outside D _{ k }. But this kind of errors can be neglected due to small subdomain size of G _{ i } and the compensation of positive and negative errors in the summation in equation (B.8).
After the distance from each point to fracture is calculated, a relation between them and the corresponding volume can be established. This way, the probability density function of the distance to the fractures is calculated. The frequency of a distance in the distribution function is the surface of the isodistance to the fracture, normalized by the volume. The area under the distribution function between two distances corresponds to the volume proportion of the subdivision delimited them. The area of the interface between two subdivisions is calculated by the derivative of the cumulative function. For a more reliable probability density function, many samples are required.
Appendix C
Reference solution validation
The efficiency of any proposed model must be compared to a reliable solution. When fractures are parallel to the grid axes (Cross Fractures Case), performing a reference solution is done easily while using very fine grid cells for fractures aperture combined with an LGR technique in order to simulate properly the matrixfracture interaction. However, when an irregular fracture distribution occurs (Diagonal Fractures Case presented in Fig. 2b) this kind of technique becomes challenging due to matrix mesh generation that conforms to complex fractures distribution.
To provide a reference solution, the matrix media is meshed into very fine grid cells and the fractures are explicitly discretized considering the two connected fractures. The matrix cells exchange with fracture nodes are estimated as described in Delorme et al. (2016). It must be mentioned that our reference solution used in this work is reliable even if the fractures are not parallel to the grid axes. To prove it, let us consider a large regular fracture network as shown in Figure C1a. Thanks to the symmetrical geometry, the simulation can be limited on a small domain described by the cross or the diagonal fractures case (Fig. C1a). The volume of the domain regarding the diagonal case is two times the volume of that presented by the cross case. Hence, the gas production from the diagonal case should be theoretically twice of the cross fractures case.
Fig. C1 Illustration of (a) an infinite regular fracture network describing Case 2/3 and (b) the comparison of the cumulative gas production from the reference solution of Case 2/3. 
Figure C1b presents the comparison of the cumulative gas production between the reference solution of the diagonal case (solid blue line) and the Explicit Discretized Model of the cross fractures case (dashed red line). In fact, Case 2 the cross case produces 0.01 × 10^{9} cft, while the diagonal one produces the double, 0.02 × 10^{9} cft. Thus, multiplying the cumulative gas production of the cross case by a factor of 2, the cumulative gas production of the diagonal case is found perfectly as shown in Figure C1b (dotted red line). Finally, as the fractures presented by the cross case are parallel to the grid axes and fine grid cells associated with LGR around the fractures are used, the reference solution from the cross case is reliable. Hence, based on results from Figure C1b, we can consider that our reference solution used for the diagonal case is reliable and accurate. In general, we believe that discretizing explicitly the fracture network while using very fine matrix grid cells provides a reliable reference solution for any irregular fracture distribution. This kind of technique was used to perform a reference solution in Section 3.1.2.2.
Acknowledgments
The authors gratefully acknowledge the support of the Lebanese American University.
References
 Barenblatt G.I., Zheltov I.P., Kochina I.N. (1960) Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata], J. Appl. Math. Mech. 24, 5, 1286–1303. [CrossRef] [Google Scholar]
 Chang M.M. (1993) Deriving the shape factor of a fractured rock matrix (No. NIPER696), National Institute for Petroleum and Energy Research, Bartlesville, OK, USA. [CrossRef] [Google Scholar]
 Cheng Y. (2012) Impact of water dynamics in fractures on the performance of hydraulically fractured wells in gasshale reservoirs, J. Can. Petrol. Technol. 51, 2, 143–151. [CrossRef] [Google Scholar]
 Cipolla C.L., Lolon E.P., Erdle J.C., Rubin B. (2010) Reservoir modeling in shalegas reservoirs, SPE Reserv. Eva. Eng. 13, 4, 638–653. [CrossRef] [Google Scholar]
 Cipolla C.L., Lolon E.P., Erdle J.C., Tathed V. (2009) Modeling well performance in shalegas reservoirs, in: SPE/EAGE Reservoir Characterization & Simulation Conference, European Association of Geoscientists & Engineers, pp. cp170. [Google Scholar]
 Coats K.H. (1989) Implicit compositional simulation of singleporosity and dualporosity reservoirs, in: SPE Symposium on Reservoir Simulation, Society of Petroleum Engineers. [Google Scholar]
 Delorme M., BossieCodreanu D., BenGharbia I., Khebzegga O., Khebzegga N., Ricois O.M. (2016) Unconventional production forecast needs integration of field hydraulic stimulation data through fracture model calibration and optimized numerical scheme, in: SPE Argentina Exploration and Production of Unconventional Resources Symposium, Society of Petroleum Engineers. [Google Scholar]
 Ding D.Y., Wu Y.S., Farah N., Wang C., Bourbiaux B. (2014) Numerical simulation of low permeability unconventional gas reservoirs, in: SPE/EAGE European Unconventional Resources Conference and Exhibition, European Association of Geoscientists & Engineers, Vol. 2014, No. 1cp399. [Google Scholar]
 Farah N. (2016) Flow modelling in low permeability unconventional reservoirs, Doctoral dissertation, Paris 6. [Google Scholar]
 Farah N., Delorme M. (2020) Unified fracture network model (UFNM) for unconventional reservoirs simulation, J. Petrol. Sci. Eng. 195, 107874. https://doi.org/10.1016/j.petrol.2020.107874. [CrossRef] [Google Scholar]
 Farah N., Delorme M., Ding D.Y., Wu Y.S., Bossie Codreanu D. (2018) Flow modeling of unconventional shale reservoirs using a DFMMINC proximity function, J. Petrol. Sci. Eng. https://doi.org/10.1016/j.petrol.2018.10.014. [Google Scholar]
 Farah N., Ding D.Y., Wu Y.S. (2015) Simulation of the impact of fracturing fluid induced formation damage in shale gas reservoirs, in: SPE Reservoir Simulation Symposium, Society of Petroleum Engineers. [Google Scholar]
 Freeman C., Moridis G.J., Ilk D., Blasingame T. (2010) A numerical study of transport and storage effects for tight gas and shale gas reservoir systems, in: International Oil and Gas Conference and Exhibition in China, Society of Petroleum Engineers. [Google Scholar]
 Hajibeygi H., Karvounis D., Jenny P. (2011) A hierarchical fracture model for the iterative multiscale finite volume method, J. Comput. Phys. 230, 24, 8729–8743. [Google Scholar]
 Jiang J., Shao Y., Younis R.M. (2014) Development of a multicontinuum multicomponent model for enhanced gas recovery and CO2 storage in fractured shale gas reservoirs, in: SPE Improved Oil Recovery Symposium, Society of Petroleum Engineers. Ji, L., Settari, A. [Google Scholar]
 KarimiFard M., Gong B., Durlofsky L.J. (2006) Generation of coarsescale continuum flow models from detailed fracture characterizations, Water Res. Res. 42, 10. [CrossRef] [Google Scholar]
 Kazemi H. (1969) Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution, SPE J. 9, 4, 451–462. [Google Scholar]
 Kazemi H., Merrill L.S. Jr, Porterfield K.L., Zeman P.R. (1976) Numerical simulation of wateroil flow in naturally fractured reservoirs, SPE J. 16, 6, 317–326. [Google Scholar]
 Khvoenkova N., Delorme M. (2011) An optimal method to model transient flows in 3D discrete fracture network, in: IAMG Conference, Vol. 2011, pp. 1238–1249. [Google Scholar]
 Li J., Du C., Zhang X. (2011) Critical evaluations of shale gas reservoir simulation approaches: single porosity and dual porosity modeling, in: SPE Middle East Unconventional Gas Conference and Exhibition, Society of Petroleum Engineers. [Google Scholar]
 Li X., Zhang D., Li S. (2015) A multicontinuum multiple flow mechanism simulator for unconventional oil and gas recovery, J. Nat. Gas Sci. Eng. 26, 652–669. https://doi.org/10.1016/j.jngse.2015.07.005. [Google Scholar]
 Lim K.T., Aziz K. (1995) Matrixfracture transfer shape factors for dualporosity simulators, J. Petrol. Sci. Eng. 13, 3–4, 169–178. [CrossRef] [Google Scholar]
 Moridis G.J., Blasingame T.A., Freeman C.M. (2010) Analysis of mechanisms of flow in fractured tightgas and shalegas reservoirs, in: SPE Latin American and Caribbean Petroleum Engineering Conference, Society of Petroleum Engineers. [Google Scholar]
 Noetinger B. (2015) A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks accounting for matrix to fracture flow, J. Comput. Phys. 283, 205–223. [Google Scholar]
 Pruess K., Karasaki K. (1982) Proximity functions for modeling fluid and heat flow in reservoirs with stochastic fracture distributions. [Google Scholar]
 Pruess K., Narasimhan T.N. (1982) A practical method for modeling fluid and heat flow in fractured porous media. [Google Scholar]
 Raterman K.T., Liu Y., Warren L. (2019) Analysis of a drained rock volume: an eagle ford example, in: Unconventional Resources Technology Conference (URTeC), Denver, Colorado, 22–24 July 2019, Society of Exploration Geophysicists, pp. 4106–4125. [Google Scholar]
 Quintard M., Whitaker S. (1996) Transport in chemically and mechanically heterogeneous porous media. I: Theoretical development of regionaveraged equations for slightly compressible singlephase flow, Adv. Water Res. 19, 1, 29–47. [CrossRef] [MathSciNet] [Google Scholar]
 Ricois O.M., Gratien J., BossieCodreanu D., Khvoenkova N., Delorme M. (2016) Advantages of an unstructured unconventional fractured media model integrated within a multiphysics computational platform, International Petroleum Technology Conference. [Google Scholar]
 Rubin B. (2010) Accurate simulation of nonDarcy flow in stimulated fractured shale reservoirs, in SPE Western Regional Meeting, Society of Petroleum Engineers. [Google Scholar]
 Tatomir A.B., Szymkiewicz A., Class H., Helmig R. (2011) Modeling two phase flow in large scale fractured porous media with an extended multiple interacting continua method, Comput. Model. Eng. Sci. 77, 2, 81. [Google Scholar]
 Van Heel A.P., Boerrigter P.M. (2006) On the shapefactor in fractured reservoir simulation, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. [Google Scholar]
 Warren J.E., Root P.J. (1963) The behavior of naturally fractured reservoirs, SPE J. 3, 3, 245–255. [Google Scholar]
 Wu Y.S., Pruess K. (1988) A multipleporosity method for simulation of naturally fractured petroleum reservoirs, SPE Reserv. Eng. 3, 1, 327–336. [CrossRef] [Google Scholar]
 Wu Y.S., Wang C., Li J., Fakcharoenphol P. (2012) Transient gas flow in unconventional gas reservoir, in: SPE Europec/EAGE Annual Conference, Society of Petroleum Engineers. [Google Scholar]
 Wu Y.S., Liu H.H., Bodvarsson G.S. (2004) A triplecontinuum approach for modeling flow and transport processes in fractured rock, J. Contam. Hydrol. 73, 1–4, 145–179. [CrossRef] [PubMed] [Google Scholar]
All Figures
Fig. 1 Schematic of the MINC concept for (a) a regular fractures network and (b) for an arbitrary/irregular fractures distribution; (after; Farah et al. (2018)). 

In the text 
Fig. 2 Illustration of (a) a regular fracture distribution with uniform block distribution, (b) an irregular fractures distribution with uniform block distribution, and (c) the Proximity Function distribution regarding the diagonal case. 

In the text 
Fig. 3 Illustration of (a) a regular fracture network where different subvolumes (subblocks) are present (Sv1a, Sv2a, Sv3a, and Sv4a), (b) the MINC6–PF model, (c) an irregular fracture network distribution occurring in a matrix grid cell where different subvolumes (subblocks) are present (Sv_{1b}, Sv_{2b}, Sv_{3b}, and Sv_{4b}) and (d) its corresponding MINC6–PF model. 

In the text 
Fig. 4 Illustration of (a) LGR reference solution, (b) the MINC6PF model, (c) the total cumulative matrix volume and per subvolume from Sv_{1a}, Sv_{2a}, Sv_{3a} and Sv_{4a}. 

In the text 
Fig. 5 Illustration of (a) an irregular fracture distribution where the reference solution is performed, (b) the MINC6PF model, (c) the total cumulative matrix volume and per subvolume from Sv_{1b}, Sv_{2b}, Sv_{3b} and Sv_{4b}. 

In the text 
Fig. 6 Illustration of the comparison of the cumulative gas production for (a) the explicit discretized model, the MINC6 Proximity Function with and without correction for 1000 days of production for the regular fractures distribution presented in Figures 4a and 4b and (b) the comparison of the cumulative gas production before and after the transmissibility correction for the irregular fractures distribution presented in Figure 5a. 

In the text 
Fig. 7 Illustration of (a) an infinite fracture network consisting in a regular distribution with a 164 ft of fractures spacing’s where; (1) represents a centered fractures scenario and (2) represents a shifted fracture network scenario compared to the mesh definition, (b) the comparison of the cumulative gas production for the explicit and the MINC proximity function model for scenario (1) and (2). 

In the text 
Fig. 8 Illustration of (a) the MINC proximity function computed into the studied grid cell and (b) the cumulative matrix volume per subvolume for model (2) and (2)–NF and (c) comparison of the cumulative gas production for case (1), (2) and (2)–NF. 

In the text 
Fig. C1 Illustration of (a) an infinite regular fracture network describing Case 2/3 and (b) the comparison of the cumulative gas production from the reference solution of Case 2/3. 

In the text 