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

: In this article, an open-source code for the simulation of ﬂuid ﬂow, 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 ﬂuid mass balance equation, speciﬁc 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 ﬁrst 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.


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 Sections 2 and 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.

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 singlephase gas flow in both Ω m and Ω f is as follows [31]: 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.

Gas Density
For the density ρ of a natural gas, modified ideal gas law applies: 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: whose coefficients are as follows: with A = 0.457235ωp pr /T 2 pr and B = 0.0777961p 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: 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.

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. Table 2. Adsorption and transport, corresponding models and domain for the gas flow; A-adsorption, T-transport.

Mechanism
Model Type Domain adsorption Langmuir, BET A matrix slip flow/diffusion Klinkenberg [36], Florence et al. [37], Javadpour et al. [38], Civan [39] T matrix non-Darcy flow Darcy-Forchheimer T fracture 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]: 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: where K n is the Knudsen number, equal to and α r is the rarefaction parameter, equal to 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]: where β is the empirical Forchheimer coefficient; for propped open hydraulic fractures, β is [44]:

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]: 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: where the units of  In Figure 7, C f , normalized to the initial value of C f0 , 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: 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 0f fracture permeabilities is introduced: 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.

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.

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: 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 = Ah, 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: 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: where the Jacobian matrix J ij (x n ), equal to ∂R i (x n )/∂x j , is calculated by automatic differentiation in MRST.

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]: 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 where A ik is the area of the intersection between the fracture and grid cell. For a 2D grid, A ik is the product of the linear intersection with the uniform formation thickness DZ. Note that the harmonic average and upwind scheme are used for the permeability and viscosity, respectively. The quantity d ik is the volumetric averaging normal distance between matrix cell and fracture plane, which can be calculated as follows: For two-dimensional structured grids, an analytical solution is available for d ik [57]. For fracture-fracture NNC, star-delta transformation is used to compute the fracturefracture transfer coefficient T f−f (jk) for a fracture intersection [58]: 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 DZ; 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]: where s is the wellbore skin factor. matrix H N Figure 9. Nodes in an EDFM for the matrix, a hydraulic fracture, and a natural fracture.

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.

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.
x y well Figure 10. Scheme for the comparison with the CMG (Computer Modelling Group Ltd.) simulator (simulations VE1A, VE1B). 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 lowconductivity 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 Figures 13 and 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. LGR with NFs EDFM with NFs LGR+EDFM w.NFs

Explicit Frac Cell EDFM Frac Cell
LGR with NFs EDFM with NFs LGR+EDFM with NFs 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.

Comparison with an In-House Simulator
ShOpen is further validated by comparing the results with those provided by an inhouse 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.

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.

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.  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.

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 Tables 2 and 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. Non-planar hydraulic fracture + Natural fractures 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 shortterm performance.

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: 1.
the algorithms of the code refer to gas adsorption, gas slippage and diffusion, non-Darcy flow, and hydromechanical coupling; 2.
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; 3.
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; 4.
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; 5.
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.