Next Article in Journal
Engineering-Geological Features Supporting a Seismic-Driven Multi-Hazard Scenario in the Lake Campotosto Area (L’Aquila, Italy)
Next Article in Special Issue
Performance Analysis of Multi-Task Deep Learning Models for Flux Regression in Discrete Fracture Networks
Previous Article in Journal
Saturated Hydraulic Conductivity Measurements in a Loam Soil Covered by Native Vegetation: Spatial and Temporal Variability in the Upper Soil Layer
Previous Article in Special Issue
Experimental Reproducibility and Natural Variability of Hydraulic Transport Properties of Fractured Sandstone Samples
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Open-Source Code for Fluid Flow Simulations in Unconventional Fractured Reservoirs

1
Craft and Hawkins Department of Petroleum Engineering, Louisiana State University, Baton Rouge, LA 70803, USA
2
Dipartimento di Ingegneria dell’Innovazione, Università del Salento, 73100 Lecce, Italy
*
Author to whom correspondence should be addressed.
Submission received: 5 January 2021 / Revised: 4 February 2021 / Accepted: 12 February 2021 / Published: 22 February 2021
(This article belongs to the Special Issue Quantitative Fractured Rock Hydrology)

Abstract

:
In this article, an open-source code for the simulation of fluid flow, including adsorption, transport, and indirect hydromechanical coupling in unconventional fractured reservoirs is described. The code leverages cutting-edge numerical modeling capabilities like automatic differentiation, stochastic fracture modeling, multicontinuum modeling, and discrete fracture models. In the fluid mass balance equation, specific physical mechanisms, unique to organic-rich source rocks, are included, like an adsorption isotherm, a dynamic permeability-correction function, and an Embedded Discrete Fracture Model (EDFM) with fracture-to-well connectivity. The code is validated against an industrial simulator and applied for a study of the performance of the Barnett shale reservoir, where adsorption, gas slippage, diffusion, indirect hydromechanical coupling, and propped fractures are considered. It is the first open-source code available to facilitate the modeling and production optimization of fractured shale-gas reservoirs. The modular design also facilitates rapid prototyping and demonstration of new models. This article also contains a quantitative analysis of the accuracy and limitations of EDFM for gas production simulation in unconventional fractured reservoirs.

1. Introduction

Unconventional gas reservoirs are gaining interest as decades of oil and natural gas production have resulted in extensive use of conventional resources. As a consequence, new technologies, like horizontal well drilling and hydraulic fracturing, are constantly being introduced, driving substantial growth of the production of these unconventional resources [1].
The exploitation of such resources is technologically demanding and expensive. There is a need for high-quality predictive models for more efficient and economic design and management of the production process. Unlike the conventional reservoirs, unconventional reservoirs are characterized by complex transport mechanism and the occurrence of multiscale fractures [2], therefore, the well-established models for flow and transport in porous media are not usable [3]. In fact, with specific reference to fractured shale, the rock mass is comprised of a matrix of dense, tiny pores of sizes down to a few nanometers, cracks, and microfractures, giving intrinsic heterogeneity and anisotropy [2,4], and superposed natural/artificial long or relatively long fractures; given these features, an accurate and reliable prediction of a well’s performance is challenging (Figure 1).
The methods proposed in the scientific literature for the prediction of gas flow in such reservoirs can be categorized into analytical, semianalytical, and numerical. The analytical methods date back to the 70s; solutions of the linearized real-gas equation were derived for simple fracture geometries [5,6]. The semianalytical methods were based on the Boundary Element Method (BEM) and found to be effective for more complex fracture networks and shale gas storage, including adsorption [7,8,9,10,11,12,13,14,15]. Although the analytical or semianalytical methods are, when applicable, fast and accurate, the numerical methods are generally preferred [15,16]. Discrete Fracture Networks (DFNs) are frequently used as geometrical/hydromechanical input for simulations of fluid flow in fractured rock masses that have multiscale fractures [17]. Usually, the DFN is limited to the larger features and is embedded into a continuum of given hydraulic properties, condensing the matrix and the network of cracks and tiny fractures [18,19].
The coupling of the DFN mesh with the continuum mesh is crucial, and unstructured meshing/gridding with Local Grid Refinement (LGR) is generally used, however, the operation is challenging, especially when dealing with intricate DFNs. With a conforming mesh, it may be difficult to efficiently simulate the sharp pressure gradients around the fractures [20]. In order to resolve such a problem, the Embedded Discrete Fracture Model (EDFM) has been recently proposed. In an EDFM, the DFN is embedded into the continuum (matrix, cracks and tiny fractures) without conforming the meshes. A structured mesh/grid is employed for the continuum; the required precision can be comfortably reached by reducing the element size. In addition, an EDFM can be easily incorporated in existing, well-established reservoir simulators [15]. In Table 1, the EDFM is compared to other methods for the simulation of gas transport in unconventional reservoirs. EDFM has been widely adopted for gas production simulation in unconventional fractured reservoirs by considering various flow mechanisms and models [21,22,23,24,25].
To the authors’ knowledge, almost all the existing codes for shale gas reservoirs are implemented in in-house simulators or commercial simulators [2]. As a better understanding of the flow mechanisms in unconventional reservoirs is gained, many emerging flow models and technologies will be developed for the effective development of unconventional reservoirs [26,27,28,29]. However, the prediction of gas transport in shale is still a hard task considering the complexity of the involved phenomena, specifically with reference to the well performance; in this respect, more flexible and open-source tools may be useful to the scientific collective for better understanding and rapid prototyping and demonstration of new simulation models. Further, although EDFM is an efficient method to solve the fluid flow problem in fractured reservoirs, Tene et al. [30] show that it may introduce large errors when dealing with low conductivity or sealing faults/fractures. The effect of such errors on the prediction of well performance is hard to define.
In this note, a numerical code (ShOpen) for the simulation of nonlinear gas flow in unconventional reservoirs with multiscale fractures is described. It is an efficient and flexible framework (ShOpen) including the MATLAB Reservoir Simulation Toolbox (MRST), an open-source reservoir simulation toolkit, and an EDFM. ShOpen can handle discrete hydraulic fractures of arbitrary geometry and location. The code also encompasses the geomechanical effects on the fractures of the DFN. The typical results of the code are compared to those of a commercial simulator and of an in-house reservoir simulator, employing unstructured grids to simulate gas desorption, gas slippage, and diffusion in shale with nonplanar hydraulic fractures. The advantages and limitations of EDFM in predicting the well performance of unconventional reservoirs are then discussed. Finally, some practical applications are reported.
The note is structured as follows: the physicomathematical statements underpinning the code and the numerical model are described in Section 2 and Section 3, respectively. In Section 4, the validation of the code is reported; in Section 5, the practical examples of application are illustrated; in Section 6, concluding remarks are provided.

2. Physicomathematical Statements

A domain Ω m dense of a porous matrix is crossed by fractures of a domain Ω f . In two dimensions (plane flow), the equations ruling the isothermal single-component single-phase gas flow in both Ω m and Ω f is as follows [31]:
t b ϕ α + ( 1 ϕ α ) m a α ρ 0 + · b k α μ p α = b q + b Ψ α ,
where α is equal to ‘m’ or ‘f’ for matrix and fracture, respectively; t is time; b is inverse formation volume factor, equal to the ratio ρ / ρ 0 , with ρ , ρ 0 gas densities at the reservoir conditions and at standard conditions, respectively; ϕ α is porosity; m a α is gas mass adsorption (negligible for the fracture); f i is the permeability correction factor for the i-th gas transport mechanism; k α is the absolute permeability; μ is the gas dynamic viscosity; p α is pressure; q is gas volumetric sink/source; and Ψ α is the leakage term between the two domains ( Ψ m = Ψ f ). The permeability k α derives from the initial absolute permeability k 0 α after correction by using factors f 1 , f 2 , f 3 , f 4 (i.e., k α = f 1 f 2 f 3 f 4 k 0 α ), condensing, respectively, the effects of slippage and diffusion, deviation from Darcy law and indirect hydromechanical coupling, and f 3 for micro/meso-cracks/fractures embedded into Ω m , f 4 for discrete fractures.
Note that in Equation (1), the gravity is neglected.
In what follows, the terms of Equation (1) are explicated.

2.1. Gas Density

For the density ρ of a natural gas, modified ideal gas law applies:
ρ = p M Z R T ,
where p is gas pressure (irrespective of the domain); M is the gas molar mass; R is the ideal gas constant (8.314 JK 1 mol 1 ); T is the reservoir temperature; and Z is the compressibility factor, a function of p and T.
For factor Z, the implicit Peng–Robinson Equation-Of-State (PR-EOS) [32] can be used, or, alternatively, empirical explicit equations can be applied. Elliot and Lira [33] propose for Z a cubic equation:
Z 3 + a 2 Z 2 + a 1 Z + a 0 = 0 ,
whose coefficients are as follows: a 0 = ( A B B 2 B 3 ) , a 1 = A 3 B 2 2 B , a 2 = B 1 , with A = 0.457235 ω p p r / T p r 2 and B = 0.0777961 p pr / T pr , where ω is the acentric factor for the species; p pr is the pseudo-reduced pressure; and T pr is the pseudo-reduced temperature, respectively, equal to p / p pc and T / T pc , where p pc and T pc are the pseudo-critical pressure and pseudo-critical temperature of a natural gas. The pseudo-critical values are weighted averages (by considering the molar fractions) of the critical values of the gases constituting the mixture.
Mahmoud (2013) [34] instead provides the following equation:
Z ( p , T ) = 0.702 exp ( 2.5 T pr ) p pr 2 5.524 exp ( 2.5 T pr ) p pr + 0.044 T pr 2 0.164 T pr + 1.15 .
In this note, an analytical solution is used for solving Equation (3) [33].
In Figure 2, an estimation of Z-factor for methane by using Equations (3) and (4) is reported.

2.2. Gas Viscosity

The Lee–Gonzalez–Eakin empirical correlation [35] is used herein for the density-dependent viscosity of a natural gas:
μ = 10 7 A 1 exp ( A 2 ρ A 3 ) ,
where A 1 = [ ( 9.379 + 0.01607 M ) T 1.5 ] / ( 209.2 + 19.26 M + T ) , A 2 = 3.448 + 986.4 / T + 0.01009 M , Y = 2.447 0.2224 A 3 , with μ in Pa·s, M in g/mol, and T in R. In Figure 3, an estimation of μ for methane using Equation (5) is reported.

2.3. Adsorption and Transport

With the development of unconventional tight reservoirs in recent years, enormous effort has been spent to simulate gas adsorption and transport in shale (Figure 1). The related mechanisms and corresponding models [36,37,38,39,40,41] are summarized in Table 2, with indications of the relative domain of application (matrix or fracture).
In ShOpen, all the adsorption and transport models of Table 2 can be easily implemented by introducing the corresponding formulations for the nonlinear gas adsorption m a and the permeability correction functions f i . In what follows, examples of the models implemented in ShOpen are described.
For the adsorption, the amount of gas molecules adsorbed in the pore wall of Kerogen in shale reservoir can be estimated by using a monolayer Langmuir isotherm and a multilayer BET isotherm. In formula, respectively [40]:
m a = ρ s ρ 0 p V L p + P L ,
m a = ρ s ρ 0 V m C p r 1 p r 1 ( n + 1 ) p r n + n p r n + 1 1 + ( C 1 ) p r C p r n + 1 ,
where V L is the Langmuir volume, i.e., the maximum adsorbed gas volume at a given temperature for an infinite pressure; P L is the Langmuir pressure, i.e., the pressure at which the adsorbed gas volume is equal to V L /2; ρ s is the bulk density of the rock matrix (rock density in what follows); V m is the BET adsorption volume; C is the BET adsorption constant; n is the number of adsorption layers; p r is the reduced pressure, equal to p / P s ; and P s is the BET pseudo-saturation pressure, equal to exp [ 7.7437 1306.5485 / ( 19.4362 + T ) ] .
In Figure 4, adsorption isotherms obtained by using Equations (6) and (7) are reported.
For the effect of slippage flow and diffusion, the permeability in the low-pressure regions around the fractures increases. In ShOpen, the correction factor f 1 of Florence et al. [37] is implemented. In formula:
f 1 = 1 + α r K n 1 + 4 K n 1 + K n ,
where K n is the Knudsen number, equal to
K n = μ 2.8284 p g π R T 2 M ϕ k 0 ,
and α r is the rarefaction parameter, equal to
α r = 128 15 π 2 tan 1 4 K n 0.4 .
In Figure 5, an estimation of f 1 for methane by using Equations (8)–(10) is reported, where the limits of the flow fields, Darcy, slip, transition, and free molecular are reported.
In case of a high Forchheimer number F oc (>0.11) in a fracture, the Darcy law is no longer applicable [42]. The corresponding correction factor f 2 is expressed as follows [43]:
f 2 = 2 1 + 1 + 4 ρ β k 0 μ 2 | p | ,
where β is the empirical Forchheimer coefficient; for propped open hydraulic fractures, β is [44]:
β = 3.2808 1.485 × 10 9 k 0 × 10 15 1.021 .

2.4. Indirect Hydromechanical Coupling

As shown in Figure 1b, a shale reservoir may have multiscale fractures. The fracture conductivity may decrease/increase with the production time for the closing/opening of the fracture in the near field of a well. The effect descends from the so-called indirect hydromechanical coupling. Specifically, the pumping/injection-induced pressure change Δ p in a fracture gives rise to a corresponding normal effective stress change Δ σ n that, in turn, drives to a normal deformation of the fracture and a change of the conductivity [2,45,46,47,48]. The effect can be properly estimated by resorting to full simulations of the hydromechanical coupling, where both fluid flow and stress/strain regime are included [49,50]. However, under the assumption of total stresses unchanged ( Δ σ n = 0 ), the decoupling can be applied, thus, Δ σ n = Δ p and relations for the correction of the fracture permeability k f can be easily applied. So, a fully hydromechanically coupled solution is not pursued, rather, these relations are introduced to incorporate the effect of effective stress changes on the fluid flow.
For micro/meso-cracks/fractures, a correction factor f 3 is furnished by Gangi [51]:
f 3 = 1 σ n E a m 3 ,
where σ n is the normal (with respect to the mean plane of the fracture) effective stress; E a is the effective modulus of the asperities, which is equal to the rock modulus multiplied by the ratio of the area of the asperities and fracture area (from one-tenth to one-hundredth of the rock bulk modulus); and m is related to roughness of the fracture walls. The stress σ n is defined as σ n = σ n α B p , where σ n is the total normal stress and α B is the Biot coefficient for a fracture. In this work, we assume the vertical overburden stress σ o is the total normal stress applied to fractured rocks.
In Figure 6, an estimation of f 3 by using Equation (13) is reported using the value suggested by [52] for unconventional rock samples.
For propped fractures, Alramahi and Sundberg (2012) [53] performed a series of experiments to measure the effect of the normal (to the fracture plane) effective stress σ n on the conductivity C f of fractures in shale having different stiffness (fracture conductivity C f is the product of the fracture permeability k f with the fracture width w, unit md-ft). The authors reported an empirical relationship of the type:
C f = 10 c 1 σ n + c 2 ,
where the units of C f , σ n are md-ft and psia, respectively, and c 1 , c 2 are coefficients related to the stiffness of the shale as follows: stiff shale: −0.00011, −0.0971; medium shale: −0.00035, 0.2396; soft shale: −0.00064, −0.4585.
Wu et al. (2019) [54] performed similar experiments on unpropped fractures. The type of the derived relations is
C f = 10 c 1 ln σ n + c 2 ,
where the units of C f , σ n are again md-ft and psia, respectively, and coefficients c 1 , c 2 are as follows: stiff shale: −0.793, 4.5618; medium shale: −0.89, 5.0725; soft shale: −1.04, 6.0216.
In Figure 7, C f , normalized to the initial value of C f 0 , is plotted against σ n for propped and unpropped fractures.
In case of an artificial (hydraulic, H) fracture, the direction of propagation is normal to the minimum total horizontal stress σ h , min . The orientations of a natural (N) fracture is instead random, even though they cluster around the mean values pertaining to the fracture sets. The total normal stress on a fracture can be approximated by the average total normal stress. Therefore, σ n can be defined as follows, respectively:
σ nH = σ h , min p σ nN = ( σ h , min + σ h , max ) / 2 p ,
where σ h , min , σ h , max may coincide with S h , min , S h , max horizontal stresses at the site if the fractures are supposed to be oriented sub-vertically. Note that α B is set to 1.
In ShOpen, a dynamical correction factor f 4 for H and N k 0 f fracture permeabilities is introduced:
f 4 = C f ( p ) C f ( p 0 ) ,
where p 0 is the initial reservoir pressure.
Typical trends of the correction factor f 4 versus p for both H and N fractures are reported in Figure 8.
In ShOpen, the three categories of fractures, i.e., micro/meso-cracks/fractures, N fractures, and H fractures, are defined according to length l f and aperture e f as follows, respectively: l f /2 < 1 m, e f < 0.1 mm; l f /2 = 1–20 m, e f = 0.1–1 mm; l f /2 = 20–100 m, e f > 1 mm.

3. The Numerical Model

ShOpen includes the following modules: the automatic differentiation module (ad-core, ad-props), the black-oil module (ad-blackoil), and the hierarchical fracture model (hfm), implemented into the open-source MATLAB Reservoir Simulation Toolbox (MRST) [55]. The Two-Point-Flux Approximated Finite Volume Method (TPFA-FVM) is applied for the solution of the governing Equations (1). A fully implicit first-order backward scheme is applied for the time discretization, in which the Jacobian matrix of the nonlinear system is calculated by the automatic differentiation modules. All the nonlinear functions for adsorption, transport, and hydromechanically coupled effects are separately handled. The large H and N fractures are explicitly modeled. The micro/meso-crack/fractures are embedded into the matrix.

3.1. Numerical Discretization

By using an implicit scheme and backward finite differences, and by subdividing Ω m in cells of volume V, Equation (1) can be expressed as follows:
ϕ V Δ t ( b n + 1 b n ) + ( 1 ϕ α ) V ρ 0 Δ t ( m a α p α n + 1 m a α p α n ) · i f i n + 1 k 0 α V μ n + 1 b n + 1 p α n + 1 V b n + 1 q n + 1 V b n + 1 Ψ α n + 1 = 0 ,
where there is no mechanical storage and ϕ α is factored out of the temporal derivative. In a two-dimensional Cartesian mesh with uniform formation thickness h, the cell volume is defined as V = A h , with A cell area.
The leakage terms Ψ α are coupling terms (as illustrated later), therefore, for both domains in Equation (18), p m and p f coexists. The matrix form of Equation (18) is as follows:
A mm A mf A fm A ff p m p f = Q m Q f ,
where Q m and Q f condense the gas volumetric sources/sinks.
Note that the correction factors f i all depend on p α . In order to handle the nonlinearities of Equation (19), Newton’s iteration method is used. The residual function R of variables x = ( x 1 , , x m ) at each iteration n can be expressed as follows:
R ( x n ) = ( x n + 1 x n ) J ( x n ) ,
where the Jacobian matrix J i j ( x n ) , equal to R i ( x n ) / x j , is calculated by automatic differentiation in MRST.

3.2. EDFM

As shown in Figure 9, in an EFDM, the matrix mesh is not necessarily conforming with the fracture elements. There are three Non-Neighbor Connections (NNCs) in an EDFM: (1) fracture–matrix, (2) fracture–fracture, (3) fracture–well. For a fracture–matrix NNC, the following equation applies [56]:
ψ m = ψ f = T f m p f n + 1 p m n + 1 ,
where T f m is a fracture-to-matrix transfer coefficient. Being that i is a node of the matrix mesh/grid and k is a neighbor node of a fracture intersecting the cell of i, the transfer coefficient is
T f m ( i k ) = k 0 , i k μ i k A i k d i k ,
where A i k is the area of the intersection between the fracture and grid cell. For a 2D grid, A i k is the product of the linear intersection with the uniform formation thickness D Z . Note that the harmonic average and upwind scheme are used for the permeability and viscosity, respectively. The quantity d i k is the volumetric averaging normal distance between matrix cell and fracture plane, which can be calculated as follows:
d i k = d i k d v V i .
For two-dimensional structured grids, an analytical solution is available for d i k [57].
For fracture–fracture NNC, star–delta transformation is used to compute the fracture–fracture transfer coefficient T f f ( j k ) for a fracture intersection [58]:
T f f ( j k ) = t j t k m = 1 N int t m , t m = A f , m 0.5 h f , m k 0 , m μ m ,
where t m is a transfer coefficient for a fracture cell m in a fracture intersection; A f , m is the cross-section area of the fracture in m, which can be calculated by multiplying the fracture aperture e f with the formation thickness D Z ; h f , m is the fracture cell length in m; and N int is the total number of fracture–fracture connections.
For a fracture–well NNC, the effective wellbore index W I w f between a fracture cell and a well and the equivalent radius r e can be expressed, respectively, as follows [56]:
W I w f = 2 π k f w f ln r e / r w + s ,
r e = 0.14 h f 2 + D Z 2 ,
where s is the wellbore skin factor.

4. Validation

For the validation of ShOpen, the results of numerical simulations for a methane (CH 4 ) reservoir are compared to the results furnished by the CMG© simulator (Computer Modelling Group Ltd., Calgary, AB, Canada, www.cmgl.ca (accessed on 23 February 2021)) (validation examples VE1A, VE1B, VE1C) and an in-house simulator with unstructured mesh (validation example VE2) [59]. For the simulations, the following data are fixed: rock density ρ s 2500 kg/m 3 ; CH 4 molecular weight, critical pressure, critical temperature, acentric factor are, respectively, 0.01604 kg/mol, 4.60 MPa, 190.6 K, 0.01142, well radius 0.1 m.

4.1. Comparison with the Commercial Simulator

With reference to Figure 10, the size of the reservoir domain is 606.6 × 606.6 m 2 . A well is located at the center of a 200 m-long hydraulic fracture parallel to the y axis. As the commercial simulator does not support all the physics presented in Section 2, only the Langmuir adsorption (Equation (6)) is considered. Compressibility factor Z and natural gas viscosity μ values are derived through interpolation of the values suggested in CMG©. Specific input data are reported in Table 3. For the well, the bottom-hole pressure (BHP) is given. The production time is 30 years. The domain is discretized with 201 n x and 65 n y cells.
EDFM is an efficient method to solve fluid flow problem in fractured reservoirs, however, Tene et al. [30] show that it may introduce large errors when dealing with low-conductivity or sealing faults/fractures. Additionally, LGR is usually required to capture the sharp pressure gradients near the fractures. To evaluate the performance of EDFM method for practical shale gas production problem, three simulation scenarios (VE1A, VE1B, and VE1C) are performed.
In VE1A, three scenarios with fracture conductivity equal to 10,000 mD-ft, 50 mD-ft, and 5 mD-ft, respectively, are simulated. In Figure 11, the results in terms of gas flow rate and cumulative production are plotted against time for ShOpen (EDFM and explicit EFM) and CMG.
The results of ShOpen/EFM are in good agreement with the results provided by CMG, whereas the results of ShOpen/EDFM consistently deviate with a significant difference of even 11%. It is apparent that the production in a reservoir with fractures of low to medium permeability cannot be efficiently simulated with EDFM, as previously stated by Tene et al. (2017) [30].
For unconventional reservoir simulations, LGR is usually required to capture the sharp pressure gradients near the fractures. In VE1B, the sensitivity of the solution provided by both ShOpen/EDFM and ShOpen/EFM to the type of grid is explored.
With reference to Figure 12, the three grid types are LGR with logarithmic refinement; a standard uniform EDFM grid [30]; and a combined EDFM + LGR grid, with an additional EDFM fracture cell superposed to the LGR grid. Note that the number of all cells is equal for all the grids (499 × 61).
Fracture conductivity C f is set to 10,000 mD-ft. All other parameters are maintained equal to VE1A. The results are compared to the results offered by CMG. In Figure 13 and Figure 14, it is shown that the results in terms of gas flow rate and cumulative gas production provided by ShOpen are in good agreement with the results of CMG. Again, with ShOpen/EFM, the obtained results match with those of CMG, whereas with ShOpen/EDFM, a 3.3% difference for the cumulative gas production is experienced. It is demonstrated therefore that EDFM can capture the sharp pressure gradient near a fracture only when LGR is applied.
In consideration of the abovementioned, in VE1C, the effect of grid refinement of fractures on the accuracy of the solution for ShOpen/EDFM is also investigated. With reference to Figure 15, six 116.74 m-long natural fractures are added in the domain; the grids of the fractures are built with and without resorting to LGR. The fracture conductivity C f of these fractures is 5 mD-ft, whereas it is kept to 10,000 mD-ft in one of the hydraulic fractures.
In Figure 16, it is shown that by using ShOpen/EDFM, a significant difference (almost 17%) with respect to CMG results is obtained. Furthermore, an underestimation of the well performance by using EDFM without LGR is experienced (3.4% difference). Again, it is further demonstrated that EDFM without LGR is not adequate for modeling reservoir models with low-permeability fractures.
To sum up, EDFM, without LGR, converges to the reference solution only when fracture conductivities are very high. In practical cases, it could be misleading to assume that all the natural fractures or faults have high permeability values, given that many of these features may be no longer active under the existing stress state in the rock of the reservoir. The use of LGR and an improved EDFM algorithm [30] for handling low-permeability fractures in the context of EDFM are all but mandatory. To relieve the limitations of EDFM in the following simulations, an empirical skin-factor and uniform grid refinement are adopted.

4.2. Comparison with an In-House Simulator

ShOpen is further validated by comparing the results with those provided by an in-house simulator [59], with reference to a reservoir with a complex networks of embedded fractures (simulation VE2). For the simulation, an unstructured mesh with LGR is used for capturing the sharp pressure gradients near the fractures. The performance of a well is considered in a reservoir with a simple “regular” fracture network (Figure 17). To reproduce the results in [59], adsorption (Equation (6)) and slippage/diffusion effects (Equation (8)) are considered in this example. With reference to the adsorption and transport mechanisms, four scenarios are considered: with adsorption and slippage/diffusion, adsorption only, slippage/diffusion only, and no mechanisms. The input data for VE2 are reported in Table 4. The production time is 10,000 days.
In Figure 18, the simulated gas flow rate and cumulative gas production are reported for both ShOpen and the in-house simulator, and for the considered scenarios. The results provided by the two codes match almost perfectly.

5. Application Examples

In order to illustrate the applicability of ShOpen to practical problems, two of the unconventional reservoirs with complex fracture networks are simulated (simulations AE1 and AE2) and the relative results presented.

5.1. Barnett Shale Reservoir

The history of the field production data of the Barnett shale reservoir is simulated with Shopen (simulation AE1). The input data are excerpted from the literature [60] and reported in Table 5. The scheme of the simulation and the mesh are shown in Figure 19. A 1100 × 290 m 2 rectangular reservoir with a thickness of 90 m is discretized by using 148 × 39 cells. Twenty-eight 94.4 m-long parallel H fractures are situated in the middle. The mutual spacing is 30.5 m. The fractures all have an aperture of 0.003 m and a permeability of 100 mD. Langmuir adsorption (Equation (6)) is considered. The depth of the reservoir is 5463 m.
In Figure 20, the pressure contours, derived from the ShOpen simulation, after 1600 days are shown.
In Figure 21a, a comparison of simulated gas production rate with the field data is reported. The results of ShOpen are evidently in good agreement with the measurements. A numerical prediction of the gas production rate is displayed in Figure 21b.

5.2. Reservoir with a Stochastic DFN

To illustrate the ShOpen capability of modular design and fast prototyping, a simulation of a reservoir having the geometry and properties of the Barnett shale reservoir, also including a 248-N-fractures stochastic DFN, is described (simulation AE2, see Figure 22). The DFN is generated by resorting to the open-source fracture generator ADFNE [61]. In addition, some of the H fractures are curvilinear. The indirect hydromechanical coupling is assumed (Equations (16)–(18)). The geomechanical parameters for shale reservoir are reported in Table 6. The remaining input data refers to Table 2 and Table 5.
To evaluate the run-time performance of our code, the domain is discretized by using 5001 × 31 cells. For a computer with a Core I9-9820X CPU, the EDFM preprocessing time and computational time of this case required 405 s and 361 s, respectively. The results of a first simulation for a scenario with the H fractures are only compared to those for a scenario including the DFN. In Figure 23, the pressure contours after 3.75 years for the two scenarios are reported. As expected, with the inclusion of the DFN, a larger Stimulated Reservoir Volume (SRV) is experienced.
The effect of the HM coupling on well performance is investigated by implementing Equations (16)–(18). In Figure 24, a comparison is given between the results of a scenario with curvilinear HFs, DFN, and HM coupling and those of the scenario with linear HFs only. In the short term, the occurrence of HM coupling leads to a reduction in gas production; however, in the long term, the discrepancy between the two scenarios is negligible. Thus, the proper consideration of the HM coupling is relevant when dealing with the short-term performance.

6. Conclusions

In this work, a numerical code and an open-source framework ShOpen for simulations of unconventional gas reservoirs are described. The model encompasses up-to-date algorithms for fluid flow and transport mechanisms. Comparisons with the results provided by a commercial simulator and an in-house simulator are also reported for the validation of the code. Application examples are also illustrated. Some statements follow:
  • the algorithms of the code refer to gas adsorption, gas slippage and diffusion, non-Darcy flow, and hydromechanical coupling;
  • with the aid of EDFM, automatic differentiation, and a modular-designed framework, the use of ShOpen can be extended to shale gas reservoirs with embedded natural and/or artificial fractures of arbitrary fracture geometries and mutual connections;
  • EDFM algorithms, implemented in ShOpen, can efficiently and accurately model stochastic DFNs, however, when dealing with low-permeability (natural) fractures and hydraulic fractures with strong gradients in the near-field, substantial errors are experienced, thus, the resort to LGR and an improved EDFM algorithm [30] are recommended;
  • shale gas adsorption and transport mechanisms have a significant impact on well performance; less dramatic is the impact of HM coupling and the complexity of the DFN;
  • ShOpen is an efficient and flexible tool for research investigations and practical applications for the implementation of nonlinearities and the fast handling of fracture networks.

Author Contributions

B.W.: conceptualization, code implementation, numerical analyses; C.F.: conceptualization (revision), editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

We would like to acknowledge discussions with Olufemi Olorode for the assistance in the implementation of EDFM. Source code and example cases are available at https://github.com/BinWang0213/MRST_Shale/tree/ShOpen (accessed on 23 February 2021).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bowker, K. Barnett Shale gas production, Fort Worth Basin: Issues and discussion. AAPG Bull. 2007, 91, 523–533. [Google Scholar] [CrossRef]
  2. Akkutlu, I.Y.; Efendiev, Y.; Vasilyeva, M.; Wang, Y. Multiscale model reduction for shale gas transport in poroelastic fractured media. J. Comput. Phys. 2018, 353, 356–376. [Google Scholar] [CrossRef]
  3. Gensterblum, Y.; Ghanizadeh, A.; Cuss, R.J.; Amann-Hildenbrand, A.; Krooss, B.M.; Clarkson, C.R.; Harrington, J.F.; Zoback, M.D. Gas transport and storage capacity in shale gas reservoirs—A review. Part A: Transport processes. J. Unconv. Oil Gas Resour. 2015, 12, 87–122. [Google Scholar] [CrossRef]
  4. Chen, L.; Jiang, Z.; Liu, K.; Gao, F. Quantitative characterization of micropore structure for organic-rich Lower Silurian shale in the Upper Yangtze Platform, South China: Implications for shale gas adsorption capacity. Adv. Geo-Energy Res. 2017, 1, 112–123. [Google Scholar] [CrossRef] [Green Version]
  5. Gringarten, A.C.; Ramey, H.J.; Raghavan, R. Unsteady-state pressure distributions created by a well with a single infinite-conductivity vertical fracture. Soc. Pet. Eng. J. 1974, 14, 347–360. [Google Scholar] [CrossRef] [Green Version]
  6. Cinco Ley, H.; Samaniego, F.; Dominguez, A.N. Transient pressure behavior for a well with a finite-conductivity vertical fracture. Soc. Pet. Eng. J. 1976, 18, 253–264. [Google Scholar] [CrossRef] [Green Version]
  7. Zuo, L.; Yu, W.; Wu, K. A fractional decline curve analysis model for shale gas reservoirs. Int. J. Coal Geol. 2016, 163, 140–148. [Google Scholar] [CrossRef]
  8. Yu, W.; Wu, K.; Sepehrnoori, K.; Xu, W. A comprehensive model for simulation of gas transport in shale formation with complex hydraulic-fracture geometry. Soc. Pet. Eng. J. 2017, 20. [Google Scholar] [CrossRef]
  9. Chen, Z.; Liao, X.; Sepehrnoori, K.; Yu, W. A semianalytical model for pressure-transient analysis of fractured wells in unconventional plays with arbitrarily distributed discrete fractures. Soc. Pet. Eng. J. 2018, 23. [Google Scholar] [CrossRef]
  10. Li, Y.; Zuo, L.; Yu, W.; Chen, Y. A fully three dimensional semi-analytical model for shale gas reservoirs with hydraulic fractures. Energies 2018, 11, 436. [Google Scholar] [CrossRef] [Green Version]
  11. Wang, B.; Feng, Y.; Du, J.; Wang, Y.; Wang, S.; Yang, R. An embedded grid-free approach for near-wellbore streamline simulation. SPE J. 2018, 23, 567–588. [Google Scholar] [CrossRef]
  12. Wang, B.; Feng, Y.; Pieraccini, S.; Scialò, S.; Fidelibus, C. Iterative coupling algorithms for large multidomain problems with the boundary element method. Int. J. Numer. Methods Eng. 2019, 117, 1–14. [Google Scholar] [CrossRef] [Green Version]
  13. Chen, Z.; Liao, X.; Zhao, X.; Dou, X.; Zhu, L. Performance of horizontal wells with fracture networks in shale gas formation. J. Pet. Sci. Eng. 2015, 133, 646–664. [Google Scholar] [CrossRef]
  14. Chen, Z.; Liao, X.; Zhao, X.; Lv, S.; Zhu, L. A semianalytical approach for obtaining type curves of multiple-fractured horizontal wells with secondary-fracture networks. SPE J. 2016, 21, 538–549. [Google Scholar] [CrossRef]
  15. Olorode, O.M.; Akkutlu, I.Y.; Efendiev, Y. Compositional reservoir-flow simulation for organic-rich gas shale. Soc. Pet. Eng. J. 2017, 22. [Google Scholar] [CrossRef]
  16. Cipolla, C.; Lolon, E.; Erdle, J.; Rubin, B. Reservoir modeling in shale-gas reservoirs. SPE Reserv. Eval. Eng. 2013, 13, 638–653. [Google Scholar] [CrossRef]
  17. de Dreuzy, J.R.; Davy, P.; Bour, O. Hydraulic properties of two-dimensional random fracture networks following power-law distributions of length and aperture. Water Resour. Res. 2002, 38, 12-1–12-9. [Google Scholar] [CrossRef]
  18. Lee, S.H.; Lough, M.F.; Jensen, C.L. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resour. Res. 2001, 37, 443–455. [Google Scholar] [CrossRef]
  19. Karimi-Fard, M.; Gong, B.; Durlofsky, L.J. Generation of coarse-scale continuum flow models from detailed fracture characterizations. Water Resour. Res. 2006, 42. [Google Scholar] [CrossRef]
  20. Karimi-Fard, M.; Durlofsky, L. A general gridding, discretization, and coarsening methodology for modeling flow in porous formations with discrete geological features. Adv. Water Resour. 2016, 96, 354–372. [Google Scholar] [CrossRef]
  21. Miao, J.; Yu, W.; Xia, Z.; Zhao, W.; Xu, Y.; Sepehrnoori, K. An Easy and Fast EDFM Method for Production Simulation in Shale Reservoirs with Complex Fracture Geometry. In Proceedings of the 2nd International Discrete Fracture Network Engineering Conference, Seattle, WA, USA, 20–22 June 2018; American Rock Mechanics Association: Richardson, TX, USA, 2018. [Google Scholar]
  22. Xu, J.; Chen, B.; Sun, B.; Jiang, R. Flow behavior of hydraulic fractured tight formations considering Pre-Darcy flow using EDFM. Fuel 2019, 241, 1145–1163. [Google Scholar] [CrossRef]
  23. Xue, X.; Yang, C.; Onishi, T.; King, M.J.; Datta-Gupta, A. Modeling Hydraulically Fractured Shale Wells Using the Fast-Marching Method With Local Grid Refinements and an Embedded Discrete Fracture Model. SPE J. 2019, 24, 2590–2608. [Google Scholar] [CrossRef]
  24. Wan, X.; Rasouli, V.; Damjanac, B.; Yu, W.; Xie, H.; Li, N.; Rabiei, M.; Miao, J.; Liu, M. Coupling of fracture model with reservoir simulation to simulate shale gas production with complex fractures and nanopores. J. Pet. Sci. Eng. 2020, 193, 107422. [Google Scholar] [CrossRef]
  25. Bai, Y.; Liu, L.; Fan, W.; Sun, H.; Huang, Z.; Yao, J. Coupled compositional flow and geomechanics modeling of fractured shale oil reservoir with confined phase behavior. J. Pet. Sci. Eng. 2021, 196, 107608. [Google Scholar] [CrossRef]
  26. Cai, J.; Lin, D.; Singh, H.; Wei, W.; Zhou, S. Shale gas transport model in 3D fractal porous media with variable pore sizes. Mar. Pet. Geol. 2018, 98, 437–447. [Google Scholar] [CrossRef]
  27. Cai, J.; Lin, D.; Singh, H.; Zhou, S.; Meng, Q.; Zhang, Q. A simple permeability model for shale gas and key insights on relative importance of various transport mechanisms. Fuel 2019, 252, 210–219. [Google Scholar] [CrossRef]
  28. Zhang, T.; Sun, S.; Song, H. Flow mechanism and simulation approaches for shale gas reservoirs: A review. Transp. Porous Media 2019, 126, 655–681. [Google Scholar] [CrossRef] [Green Version]
  29. Ning, X.; Feng, Y.; Wang, B. Numerical simulation of channel fracturing technology in developing shale gas reservoirs. J. Nat. Gas Sci. Eng. 2020, 83, 103515. [Google Scholar] [CrossRef]
  30. Ţene, M.; Bosma, S.B.; Al Kobaisi, M.S.; Hajibeygi, H. Projection-based Embedded Discrete Fracture Model (pEDFM). Adv. Water Resour. 2017, 105, 205–216. [Google Scholar] [CrossRef]
  31. Yang, R.; Huang, Z.; Yu, W.; Li, G.; Ren, W.; Zuo, L.; Tan, X.; Sepehrnoori, K.; Tian, S.; Sheng, M. A comprehensive model for real gas transport in shale formations with complex non-planar fracture networks. Sci. Rep. 2016, 6, 36673. [Google Scholar] [CrossRef] [Green Version]
  32. Peng, D.Y.; Robinson, D.B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15, 59–64. [Google Scholar] [CrossRef]
  33. Elliott, J.R.; Lira, C.T. Introductory Chemical Engineering Thermodynamics; Prentice-Hall International Series in the Physical and Chemical Engineering Sciences; Prentice Hall: Upper Saddle River, NJ, USA, 2012. [Google Scholar]
  34. Mahmoud, M. Development of a new correlation of gas compressibility factor (Z-factor) for high pressure gas reservoir. J. Energy Resour. Technol. 2013, 136, 012903. [Google Scholar] [CrossRef]
  35. Lee, A.L.; Gonzalez, M.H.; Eakin, B.E. The viscosity of natural gases. J. Pet. Technol. 1966, 18, 997–1000. [Google Scholar] [CrossRef]
  36. Klinkenberg, L.J. The permeability of porous media to liquids and gases. In Drilling and Production Practice; American Petroleum Institute: New York, NY, USA, 1941. [Google Scholar]
  37. Florence, F.A.; Rushing, J.; Newsham, K.E.; Blasingame, T.A. Improved permeability prediction relations for low permeability sands. In Proceedings of the SPE Rocky Mountain Oil & Gas Technology Symposium, Denver, CO, USA, 16–18 April 2007; Society of Petroleum Engineers: Tyler, TX, USA, 2007. [Google Scholar]
  38. Javadpour, F.; Fisher, D.; Unsworth, M. Nanoscale gas flow in shale gas sediments. J. Can. Pet. Technol. 2007, 46. [Google Scholar] [CrossRef]
  39. Civan, F. Effective correlation of apparent gas permeability in tight porous media. Transp. Porous Media 2010, 82, 375–384. [Google Scholar] [CrossRef]
  40. Yu, W.; Sepehrnoori, K.; Patzek, T.W. Modeling gas adsorption in Marcellus shale with Langmuir and BET isotherms. Soc. Pet. Eng. J. 2016, 21, 589–600. [Google Scholar] [CrossRef] [Green Version]
  41. Shen, W.; Li, X.; Cihan, A.; Lu, X.; Liu, X. Experimental and numerical simulation of water adsorption and diffusion in shale gas reservoir rocks. Adv. Geo-Energy Res. 2019, 3, 165–174. [Google Scholar] [CrossRef]
  42. Zeng, Z.; Grigg, R. A criterion for non-Darcy flow in porous media. Transp. Porous Media 2006, 63, 57–69. [Google Scholar] [CrossRef]
  43. Barree, R.D.; Conway, M.W. Beyond beta factors: A complete model for Darcy, Forchheimer, and trans-Forchheimer flow in porous media. In Proceedings of the SPE Annual Technical Conference and Exhibition, Houston, TX, USA, 26–29 September 2004; Society of Petroleum Engineers: Tyler, TX, USA, 2004. [Google Scholar]
  44. Rubin, B. Accurate simulation of non-Darcy flow in stimulated fractured shale reservoirs. In Proceedings of the SPE Western Regional Meeting, Anaheim, CA, USA, 27–29 May 2010; Society of Petroleum Engineers: Tyler, TX, USA, 2010. [Google Scholar]
  45. Hu, X.; Wu, K.; Li, G.; Tang, J.; Shen, Z. Effect of proppant addition schedule on the proppant distribution in a straight fracture for slickwater treatment. J. Pet. Sci. Eng. 2018, 167, 110–119. [Google Scholar] [CrossRef]
  46. Wang, H.F. Theory of Linear Poroelasticity; Princeton University Press: Princeton, NJ, USA, 2000. [Google Scholar]
  47. Rutqvist, J.; Stephansson, O. The role of hydromechanical coupling in fractured rock engineering. Hydrogeol. J. 2003, 11, 7–40. [Google Scholar] [CrossRef] [Green Version]
  48. Cammarata, G.; Fidelibus, C.; Cravero, M.; Barla, G. The hydro-mechanically coupled response of rock fractures. Rock Mech. Rock Eng. 2007, 40, 41–61. [Google Scholar] [CrossRef]
  49. Fidelibus, C. The 2D hydro-mechanically coupled response of a rock mass with fractures via a mixed BEM–FEM technique. Int. J. Numer. Anal. Methods Geomech. 2007, 31, 1329–1348. [Google Scholar] [CrossRef]
  50. Hu, M.; Wang, Y.; Rutqvist, J. Fully coupled hydro-mechanical numerical manifold modeling of porous rock with dominant fractures. Acta Geotech. 2017, 12, 231–252. [Google Scholar] [CrossRef] [Green Version]
  51. Gangi, A.F. Variation of whole and fractured porous rock permeability with confining pressure. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1978, 15, 249–257. [Google Scholar] [CrossRef]
  52. Shi, J.Q.; Durucan, S. Near-exponential relationship between effective stress and permeability of porous rocks revealed in Gangi’s phenomenological models and application to gas shales. Int. J. Coal Geol. 2016, 154, 111–122. [Google Scholar] [CrossRef] [Green Version]
  53. Alramahi, B.; Sundberg, M.I. Proppant embedment and conductivity of hydraulic fractures in shales. In Proceedings of the 46th ARMA US Rock Mechanics/Geomechanics Symposium, Chicago, IL, USA, 24–27 June 2012; Society of Petroleum Engineers: Tyler, TX, USA, 2012. [Google Scholar]
  54. Wu, W.; Zhou, J.; Kakkar, P.; Russell, R.; Sharma, M.M. An experimental study on conductivity of unpropped fractures in preserved shales. Soc. Pet. Eng. J. 2019, 34, 280–296. [Google Scholar] [CrossRef]
  55. Lie, K.; Krogstad, S.; Ligaarden, I.S.; Natvig, J.R.; Nilsen, H.M.; Skaflestad, B. Open-source MATLAB implementation of consistent discretisations on complex grids. Comput. Geosci. 2012, 16, 297–322. [Google Scholar] [CrossRef] [Green Version]
  56. Xu, Y. Implementation and application of the Embedded Discrete Fracture Model (EDFM) for Reservoir Simulation in Fractured Reservoirs. Master’s Thesis, The University of Texas at Austin, Austin, TX, USA, 2015. [Google Scholar]
  57. Ţene, M.; Al Kobaisi, M.S.; Hajibeygi, H. Algebraic multiscale method for flow in heterogeneous porous media with embedded discrete fractures (F-AMS). J. Comput. Phys. 2016, 321, 819–845. [Google Scholar] [CrossRef] [Green Version]
  58. Karimi-Fard, M.; Durlofsky, L.J.; Aziz, K. An efficient discrete fracture model applicable for general purpose reservoir simulators. In Proceedings of the SPE Reservoir Simulation Symposium, Houston, TX, USA, 3–5 February 2003; Society of Petroleum Engineers: Tyler, TX, USA, 2003. [Google Scholar]
  59. Jiang, J.; Younis, R.M. Numerical study of complex fracture geometries for unconventional gas reservoirs using a discrete fracture-matrix model. J. Nat. Gas Sci. Eng. 2015, 26, 1174–1186. [Google Scholar] [CrossRef]
  60. Cao, P.; Liu, J.; Leong, Y.K. A fully coupled multiscale shale deformation-gas transport model for the evaluation of shale gas extraction. Fuel 2016, 178, 103–117. [Google Scholar] [CrossRef]
  61. Fadakar Alghalandis, Y. ADFNE: Open source software for discrete fracture network engineering, two and three dimensional applications. Comput. Geosci. 2017, 102, 1–11. [Google Scholar] [CrossRef]
Figure 1. Shale gas production: (a) gas transport, (b) multiscale features.
Figure 1. Shale gas production: (a) gas transport, (b) multiscale features.
Geosciences 11 00106 g001
Figure 2. Natural gas Z-factor for methane estimated by Peng–Robinson Equation-Of-State (PR-EOS) [32] and empirical equation [34]; T = 352 K, T c = 191 K, p c = 4.64 MPa, R = 8.314 JK 1 mol 1 .
Figure 2. Natural gas Z-factor for methane estimated by Peng–Robinson Equation-Of-State (PR-EOS) [32] and empirical equation [34]; T = 352 K, T c = 191 K, p c = 4.64 MPa, R = 8.314 JK 1 mol 1 .
Geosciences 11 00106 g002
Figure 3. Methane viscosity μ (in cp, 1 cp = 1 mPa·s) estimated by the Lee–Gonzalez–Eakin empirical correlation [35] with M = 16.04 g/mol and T = 633.6 R.
Figure 3. Methane viscosity μ (in cp, 1 cp = 1 mPa·s) estimated by the Lee–Gonzalez–Eakin empirical correlation [35] with M = 16.04 g/mol and T = 633.6 R.
Geosciences 11 00106 g003
Figure 4. Langmuir and BET isotherms for V L = 0.0031 m 3 /kg, P L = 7.89 MPa, T = 327.59 K, P s = 53.45 MPa, V m = 0.0015 m 3 /kg, C = 24.56, and n = 4.46.
Figure 4. Langmuir and BET isotherms for V L = 0.0031 m 3 /kg, P L = 7.89 MPa, T = 327.59 K, P s = 53.45 MPa, V m = 0.0015 m 3 /kg, C = 24.56, and n = 4.46.
Geosciences 11 00106 g004
Figure 5. Permeability correction factor f 1 plotted against the Knudesen number K n for methane (T = 191 K, k 0 = 1 × 10 10 m 2 , ϕ = 0.1) with flow fields.
Figure 5. Permeability correction factor f 1 plotted against the Knudesen number K n for methane (T = 191 K, k 0 = 1 × 10 10 m 2 , ϕ = 0.1) with flow fields.
Geosciences 11 00106 g005
Figure 6. Permeability correction factor f 3 versus gas pressure with m = 0.5, α B = 0.5, E a = 180 MPa, σ n = 38 MPa, and α B = 0.5.
Figure 6. Permeability correction factor f 3 versus gas pressure with m = 0.5, α B = 0.5, E a = 180 MPa, σ n = 38 MPa, and α B = 0.5.
Geosciences 11 00106 g006
Figure 7. Normalized fracture conductivity C f / C f 0 plotted against the effective normal stress σ n on a fracture for (a) propped fractures and (b) unpropped fractures.
Figure 7. Normalized fracture conductivity C f / C f 0 plotted against the effective normal stress σ n on a fracture for (a) propped fractures and (b) unpropped fractures.
Geosciences 11 00106 g007
Figure 8. Permeability correction factor f 4 for hydraulic fractures and natural fractures with pressure p 0 = 34.5 MPa, σ h , min = 29 MPa, and σ h , max = 34 MPa.
Figure 8. Permeability correction factor f 4 for hydraulic fractures and natural fractures with pressure p 0 = 34.5 MPa, σ h , min = 29 MPa, and σ h , max = 34 MPa.
Geosciences 11 00106 g008
Figure 9. Nodes in an EDFM for the matrix, a hydraulic fracture, and a natural fracture.
Figure 9. Nodes in an EDFM for the matrix, a hydraulic fracture, and a natural fracture.
Geosciences 11 00106 g009
Figure 10. Scheme for the comparison with the CMG (Computer Modelling Group Ltd.) simulator (simulations VE1A, VE1B).
Figure 10. Scheme for the comparison with the CMG (Computer Modelling Group Ltd.) simulator (simulations VE1A, VE1B).
Geosciences 11 00106 g010
Figure 11. Results of the VE1A simulation: (a) gas flow rate and (b) cumulative gas production versus time t for fracture conductivity C f equal to 5 (green lines), 50 (blue lines) and 10,000 (red lines) mD-ft. Solid lines—ShOpen/EDFM; dashed lines—ShOpen/EFM; circles—CMG.
Figure 11. Results of the VE1A simulation: (a) gas flow rate and (b) cumulative gas production versus time t for fracture conductivity C f equal to 5 (green lines), 50 (blue lines) and 10,000 (red lines) mD-ft. Solid lines—ShOpen/EDFM; dashed lines—ShOpen/EFM; circles—CMG.
Geosciences 11 00106 g011
Figure 12. Grids for VE1B sensitivity analysis simulations: (a) Local Grid Refinement (LGR), (b) EDFM, (c) LGR+EDFM; the fracture cell is magnified 10 times.
Figure 12. Grids for VE1B sensitivity analysis simulations: (a) Local Grid Refinement (LGR), (b) EDFM, (c) LGR+EDFM; the fracture cell is magnified 10 times.
Geosciences 11 00106 g012
Figure 13. Comparison of the results provided by ShOpen and CMG in terms of (a) gas flow rate and (b) cumulative production, for VE1B with C f 10,000 mD-ft.
Figure 13. Comparison of the results provided by ShOpen and CMG in terms of (a) gas flow rate and (b) cumulative production, for VE1B with C f 10,000 mD-ft.
Geosciences 11 00106 g013
Figure 14. Comparison of the results provided by ShOpen and CMG in terms of (a) gas flow rate and (b) cumulative production, for VE1B with C f 5 mD-ft.
Figure 14. Comparison of the results provided by ShOpen and CMG in terms of (a) gas flow rate and (b) cumulative production, for VE1B with C f 5 mD-ft.
Geosciences 11 00106 g014
Figure 15. Grids for VE1C simulation: (a) LGR, (b) EDFM, (c) LGR+EDFM.
Figure 15. Grids for VE1C simulation: (a) LGR, (b) EDFM, (c) LGR+EDFM.
Geosciences 11 00106 g015
Figure 16. Comparison of the results provided by using the different grids in terms of (a) gas flow rate and (b) normalized cumulative gas production, for VE1C.
Figure 16. Comparison of the results provided by using the different grids in terms of (a) gas flow rate and (b) normalized cumulative gas production, for VE1C.
Geosciences 11 00106 g016
Figure 17. Scheme for the comparison with the in-house simulator (simulation VE2).
Figure 17. Scheme for the comparison with the in-house simulator (simulation VE2).
Geosciences 11 00106 g017
Figure 18. Gas flow rate (a) of VE2 and cumulative gas production (b) for adsorption and slippage/diffusion, adsorption only, slippage/diffusion only, and no mechanisms; comparison between ShOpen and the in-house simulator [59].
Figure 18. Gas flow rate (a) of VE2 and cumulative gas production (b) for adsorption and slippage/diffusion, adsorption only, slippage/diffusion only, and no mechanisms; comparison between ShOpen and the in-house simulator [59].
Geosciences 11 00106 g018
Figure 19. AE1 simulation scheme (top) and relative EDFM grid with 28 linear HFs (bottom).
Figure 19. AE1 simulation scheme (top) and relative EDFM grid with 28 linear HFs (bottom).
Geosciences 11 00106 g019
Figure 20. Pressure contours after 1600 days for the Barnett shale reservoir (simulation AE1).
Figure 20. Pressure contours after 1600 days for the Barnett shale reservoir (simulation AE1).
Geosciences 11 00106 g020
Figure 21. History matching (a) and prediction (b) of the gas production rate in the Barnett shale reservoir (simulation AE1).
Figure 21. History matching (a) and prediction (b) of the gas production rate in the Barnett shale reservoir (simulation AE1).
Geosciences 11 00106 g021
Figure 22. Scheme of simulation AE2 showing the 28 HFs and the 248 NFs.
Figure 22. Scheme of simulation AE2 showing the 28 HFs and the 248 NFs.
Geosciences 11 00106 g022
Figure 23. Pressure contours at 3.75 years for the Barnett shale reservoir (simulation AE2) for the scenario with curvilinear HFs (top) and for the scenario with the DFN added (bottom).
Figure 23. Pressure contours at 3.75 years for the Barnett shale reservoir (simulation AE2) for the scenario with curvilinear HFs (top) and for the scenario with the DFN added (bottom).
Geosciences 11 00106 g023
Figure 24. Comparison of gas flow rate (a) and cumulative production (b) between the scenario with linear HFs and the scenario with curvilinear HFs and DFN and, in addition, consideration of HMC (simulation AE2).
Figure 24. Comparison of gas flow rate (a) and cumulative production (b) between the scenario with linear HFs and the scenario with curvilinear HFs and DFN and, in addition, consideration of HMC (simulation AE2).
Geosciences 11 00106 g024
Table 1. Comparison of methods for the prediction of gas transport in unconventional reservoirs; 1 : nonlinear gas transport and storage, multiphase and multicomponent flow. DFN—Discrete Fracture Networks, EDFM—Embedded Discrete Fracture Model.
Table 1. Comparison of methods for the prediction of gas transport in unconventional reservoirs; 1 : nonlinear gas transport and storage, multiphase and multicomponent flow. DFN—Discrete Fracture Networks, EDFM—Embedded Discrete Fracture Model.
AnalyticalSemianalyticalStructured GridUnstructured GridEDFM
accuracy++++++++++++
handling nonlinear mechanisms 1 +++++++++++
handling rock heterogeneity+++++++++++
quality of DFN mesh+++++++++++
preprocessing efficiency++++++++++++++
computational efficiency+++++++++++
Table 2. Adsorption and transport, corresponding models and domain for the gas flow; A—adsorption, T—transport.
Table 2. Adsorption and transport, corresponding models and domain for the gas flow; A—adsorption, T—transport.
MechanismModelTypeDomain
adsorptionLangmuir, BETAmatrix
slip flow/diffusionKlinkenberg [36], Florence et al. [37], Javadpour et al. [38], Civan [39]Tmatrix
non-Darcy flowDarcy–ForchheimerTfracture
Table 3. Input data of the simulations for the comparison with CMG. BHP—bottom-hole pressure.
Table 3. Input data of the simulations for the comparison with CMG. BHP—bottom-hole pressure.
PropertyUnitValue
domain dimensions Δ x × Δ y m 2 606.6 × 606.6
formation thickness D Z m45.72
initial reservoir pressure p 0 MPa34.47
reservoir temperature TK327.60
Langmuir pressure P L MPa8.96
Langmuir volume V L m 3 /kg0.0041
matrix porosity ϕ m -0.07
matrix compressibility1/Pa1.45 × 10 10
matrix permeability k m nD500
fracture permeability k f mD0.5–1000
fracture width wm0.003
fracture half-length ½ l f m106.68
fracture conductivity C f mD-ft5–10,000
well BHPMPa3.45
Table 4. Key reservoir and simulation parameters of VE2.
Table 4. Key reservoir and simulation parameters of VE2.
PropertyUnitValue
domain dimensions Δ x × Δ y m200,140
formation thickness D Z m10
initial reservoir pressure p 0 MPa16
reservoir temperature TK343.15
Langmuir pressure P L MPa4
Langmuir volume V L m 3 /kg0.018
matrix porosity ϕ m -0.1
matrix compressibility1/Pa1.0 × 10 9
matrix permeability k m nD100
fracture porosity ϕ f -1.0
fracture permeability k f D1
fracture width wm 10 3
well BHPMPa4
wellbore skin factor s-43
Table 5. Input data for Barnett shale reservoir simulation AE1 [60]; other input in Table 2.
Table 5. Input data for Barnett shale reservoir simulation AE1 [60]; other input in Table 2.
PropertyUnitValue
domain dimensions Δ x × Δ y m1200,300
formation thickness D Z m90
initial reservoir pressure p 0 MPa20.34
reservoir temperature TK352
rock density ρ s kg/m 3 2500
Langmuir pressure P L MPa4.47
Langmuir volume V L m 3 /kg0.00272
matrix porosity ϕ m ϕ m -0.1
matrix compressibility1/Pa1.0 × 10 9
matrix permeability k m nD200
fracture porosity ϕ f -0.03
fracture permeability k f D1.0 × 10 8
fracture width wm0.003
fracture half-length ½ l f m47.2
fracture conductivity C f mD-ft1
well bottom-hole pressure BHPMPa3.69
wellbore skin factor s-19
Table 6. Geomechanical parameters of Barnett shale for simulation AE2; the other input data are shown in Table 2 and Table 5.
Table 6. Geomechanical parameters of Barnett shale for simulation AE2; the other input data are shown in Table 2 and Table 5.
PropertyUnitValue
Biot coefficient α B -0.5
overburden stress σ 0 MPa38
maximum horizontal stress S h , max MPa34
minimum horizontal stress S h , min MPa29
effective modulus of the asperities E a MPa180
Gangi exponential constant m-0.5
N fracture permeability k f mD10
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wang, B.; Fidelibus, C. An Open-Source Code for Fluid Flow Simulations in Unconventional Fractured Reservoirs. Geosciences 2021, 11, 106. https://0-doi-org.brum.beds.ac.uk/10.3390/geosciences11020106

AMA Style

Wang B, Fidelibus C. An Open-Source Code for Fluid Flow Simulations in Unconventional Fractured Reservoirs. Geosciences. 2021; 11(2):106. https://0-doi-org.brum.beds.ac.uk/10.3390/geosciences11020106

Chicago/Turabian Style

Wang, Bin, and Corrado Fidelibus. 2021. "An Open-Source Code for Fluid Flow Simulations in Unconventional Fractured Reservoirs" Geosciences 11, no. 2: 106. https://0-doi-org.brum.beds.ac.uk/10.3390/geosciences11020106

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop