© N. Farah & A. Ghadboun, published by IFP Energies nouvelles, 2021

1 Introduction

Oil and gas production from unconventional resources (shale gas, tight oil) is, for the most part, uneconomical due to the low-permeability and high-heterogeneity 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 in-house reservoir simulators that were initially developed for conventional reservoirs. This led to a limited understanding of unconventional reservoirs where flow from multi-scaled fractured reservoirs remains a challenge to modeling approaches.

Earlier models such as single-porosity, dual-continuum (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 multi-continuum models were all used for flow modeling in unconventional reservoirs.

An explicit discretized model with a single-porosity 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 multi-scale 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 CPU-time consuming due to the large number of grid cells used.

On the other hand, dual-continuum 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; Karimi-Fard et al., 2006; Li et al., 2015; Tatomir et al., 2011). Moreover, in cases of multi-phase (gas/water or oil/gas) flow simulations, a challenging task is to predict the GOR (Gas-Oil Ratio). Due to the large grid blocks used by dual-continuum 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 matrix-fracture exchanges was highlighted.

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 non-uniformity

Figures 2a and 2b illustrate a regular and an irregular fractures distribution, respectively, where uniform sub-volumes (Sv₁, Sv₂, Sv₃ and Sv₄ are equal) are present. The MINC–PF (see Appendix B: MINC-PF 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 sub-volume (Sv₁, Sv₂, Sv₃ and Sv₄) present within the studied grid cell. Moreover, the distribution functions from each sub-volume overlap as presented in Figure 2c. In fact, Figure 2c represents the distribution functions from each sub-volume 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).

When handling realistic cases, non-uniform sub-volumes/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 MINC-PF 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 sub-volumes/blocks was not considered. This led to a miscalculation when modeling inter-porosity flow, thus causing a difference between the DFM results and the reference solution.

Due to the non-uniformity of fractures presented in Figures 3a and 3c, non-uniform sub-volumes (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 sub-volume since computing the MINC–PF depends on the distance to the nearest fracture. In other words, some sub-volumes 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 sub-volume ends in the total cumulative distribution curve, no more refinement can be done in this sub-volume. 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 sub-volumes. Unfortunately, in the present MINC–PF model implemented in Farah et al. (2018) the 6th subdivision from these sub-volumes exchanges flow with all subdivisions #5 from the other sub-volumes (Sv_1a, Sv_2a, Sv_3a and Sv_2b, Sv_3b). Hence, these connections overestimate the exchange area between different subdivisions leading to non-accurate flow modeling results when compared to a reference solution. So, in case of non-uniform sub-blocks 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:

1. Dividing each sub-volume into 6 subdivisions.

2. Recalculating transmissibility values within the studied grid cell taking the assumption of no-flow boundary at its borders. This solution is a weighted discretization of the matrix with adapted transmissibility when non-uniform block sizes exist within the studied grid cell.

3. 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 non-uniformity

The first solution will be disregarded because it increases the simulation CPU time. Note that, computing a MINC6-PF 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 sub-volumes 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 – non-uniform grid blocks

3.1.1 Method

For a steady-state 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 sub-blocks existing within the studied grid cell.

Moreover, the total flow between two continuums (matrix-fracture or matrix-matrix) is equal to the sum of the flow from these continuums considering all existing sub-volumes in a grid cell. So, the transmissibility between continuum #i and #i + 1 considering all the sub-volumes are described as follows:

(2)where NV represents the number of sub-volumes of the studied grid cell respectively and represents the transmissibility between #i and #i + 1 from different sub-blocks.

To overcome the earlier discussed error, the connection must be adjusted by correcting the transmissibility from the last subdivision to the fracture under a steady-state regime. The correction could be expressed as follows:

(3)where,

• 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 Sub-Volume).

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 non-automated general procedure is applied. In both cases presented earlier in Figures 3a and 3c, regular and irregular fracture networks, the 6th subdivision from sub-volumes Sv_4a, Sv_1b, and Sv_4b exchanges flow with all subdivisions #5 from different sub-volumes 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 non-symmetric 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.

Figure 4c presents the cumulative matrix volume of the distribution function for the total matrix grid cell and per sub-volume 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 sub-volume (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 sub-volume (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:

(4)

(5)

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⁻⁴ 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⁶ 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.

Figure 5c presents the cumulative matrix volume distribution function for the total matrix grid cell (red solid line) and per sub-volume 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 sub-volumes Sv_1b and Sv_4b due to the fracture distribution and to the bigger volume presented by these sub-volumes (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 sub-volumes Sv_1b and Sv_4b exchange flow with all subdivisions #5 from different sub-volumes (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 sub-volumes needing adjustment are CSV = 1 and 4. Thus, equations (6) and (7) are obtained.

(6)

(7)

The corrected transmissibility can be obtained from equation (7) being the only unknown value. The transmissibility is calculated with all four sub-volumes when the distance is smaller than 10 ft. (interaction length), but has to be computed with only sub-volumes Sv_1b and Sv_4b when the distance is larger than 10 ft. No further refinement can be done into a matrix sub-volume 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.

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 MINC8-PF 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 sub-volume (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 MINC-PF 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 large-scale cases and more complex fracture networks, these corrections must be detected and corrected automatically via our pre-processor. This would be possible after we integrate the recalculation method into the DFM-MINC Proximity Function developed in Farah et al. (2018).

3.2 Adjusting the proximity function

This solution is characterized by simulating the matrix-fracture flow exchange considering the presence of nearby fractures (outside the studied grid cell). In this part, this solution is applied on a non-uniform blocks distribution for regular fracture network and a single-phase 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).

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 sub-volume 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 matrix-fracture 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.

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 MINC-PF 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 pt-NF”) is smaller than the distance to the fractures inside the studied grid cell (denoted by “d pt-frac 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 MINC-PF accurately simulates matrix-fracture 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 non-uniform grid blocks (that leads to a miscalculation when modeling inter-porosity flow using the MINC-PF model) is solved using one of the methods presented above. The solution consists in either transmissibility recalculation within the needed sub-volumes 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 pre-processor 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 multi-scale fractures embedded in low-permeability matrix. Neither an explicit discretized model using a single-porosity with an LGR technique nor a dual-continuum 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 matrix-fracture 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 matrix-fracture 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 sub-volumes/blocks within the studied grid cell are faced. However, due to the complex fracture network distribution and the presence of non-uniform sub-volumes/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 MINC-PF will increase the accuracy of the DFM-MINC Proximity Function developed in Farah et al (2018). Our future work will consist in testing more complex scenarios (large-scale cases) while exploring both methods after introducing them into the DFM-MINC Proximity function.

Appendix A

Governing equations

The general equation governing three-phase, multi-component, three-dimensional 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 control-volume method, where time discretization is carried out using a backward, first order, fully implicit, finite-difference 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 ⁿ⁺¹ 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 sub-cells 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 sub-cells and D _ij is the average distance between the two sub grids.

Appendix B

MINC-PF 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:

(B.2)

The N subdivisions of MINC can therefore be defined as follows:

(B.3)where r ₁ ≤ r ₂ ≤…≤ r _N .

The corresponding volumes are given by:

(B.4)

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:

(B.5)

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 sub-domains 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:

(B.8)

On the other hand, we have:

(B.9)

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 sub-domain 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 sub-domain G _i is neither inside nor outside D _k . But this kind of errors can be neglected due to small sub-domain 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 iso-distance 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 matrix-fracture 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.

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⁹ cft, while the diagonal one produces the double, 0.02 × 10⁹ 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

All Figures