Modeling of gas flow in confined formations at different scales

Gas flow in fractured nano-porous shale formations is complicated by a hierarchy of structural features, ranging from nanopores to hydraulic fractures, and by several transport mechanisms that differ from standard viscous flow used in reservoir modeling. The use of accurate simulation techniques that honor the physical complexity of these reservoirs and capture the associated dynamics of nanopores is required. However, these simulations often necessitate a large amount of computational resources for field scale models and therefore require upscaling. Usually, the upscaling techniques are based on idealizations that do not reflect the discrete features of the reservoir. In this work, we first incorporate the physics model that describe dynamics of shale gas into a numerical Discrete Fracture and Matrix (DFM) model. The formulation of our DFM model applies an unstructured control volume finite difference approach with a two-point flux approximation. We then propose to upscale these detailed descriptions using two different techniques, with the major difference in their coarse-grid geometry. The first approach, referred to as Embedded DFM upscaling, relies on a structured Cartesian coarse grid. The second method, which we call the Multiple Sub-Regions (MSR) upscaling, introduces a flow based coarse grid to replicate the diffusive character of the pressure in the matrix. The required parameters for the coarse-scale model in both methods and the geometry of the subregions in the second method are determined using numerical homogenization of the underlying discrete fracture model. An accurate comparison with the fine-scale representation indicates an existence of an additional transient phenomenon at coarse scale. To account for this effect, the transmissibility of both types of coarse models is related to the pressure in our approach. Both upscaling methods are applied to simulate a shale-gas flow in 2D fractured reservoir models and are shown to provide results in close agreement with the underlying fine-scale model and with a considerable reduction in the


Introduction
The abundance of shale gas in the world has caught the attention of many countries looking for an alternative to the declining conventional resources of natural gas. The U.S Energy Information Administration ( http://www.eia.gov/ EIA) estimates that known shale-gas deposits worldwide add 47% to the global technically recoverable natural gas resources. As shale-gas reservoirs are characterized by a very low permeability, gas production can only be achieved by stimulating the shale formation through hydraulic fracturing. Combined with horizontal drilling, it maximizes extraction by allowing multiple fractures along the shale bed, which in turn allows the gas to be commercially extracted. The large development cost of shale-gas reservoirs makes it a necessity to develop accurate and reliable modeling tools for uncertainty quantification and risks assessment.
The behavior of fluid flow in a shale matrix also deviates from that in a conventional gas reservoir/rock matrix as matrix porosity in shale consists of pores within the nanometer scale [9,23]. As the size of the confining pore space reduces to nano-scale, the validity of the standard approach based on the Navier-Stokes equation with no-slip boundary condition diminishes [16]. Researchers (e.g. Javadpour et al. [22], Darabi et al. [10] as well as others) identified the main transport mechanisms in shale gas as viscous flow and self-diffusion due to gas expansion. Additionally, as pore diameter becomes of the order of the molecular mean free path, the molecule-wall collisions become more pronounced, also known as Knudsen diffusion [40]. The storage of gas molecules in organic-rich shale sediments also differs from conventional resources. The gas is stored as a compressed gas in the pores, but also as an adsorbed gas to the pore walls and as a dissolved gas in the solid organic matter of the shale (i.e., kerogen and clays) [22,37]. The majority of studies investigated the shale-gas dynamics by means of an apparent permeability model [22,23,8,39] that combine the different transport mechanisms written in a format compatible with Darcy formulation. Other modeling approaches, such as molecular dynamics [6], direct simulation Monte Carlo [28] and Lattice-Boltzmann [19] can be used to model gas flow in nanopores. However, these modeling techniques are computationally demanding and cannot be used for systems larger than a few microns. One promising approach is a direct incorporation of confinement effects into the thermodynamic description proposed in [41].
Furthermore, the analysis of gas production from shale formation is complicated by a hierarchy of structural features induced by the multistage hydraulic fracturing and the consequent micro-seismic events which creates a connected network of secondary fractures within the reservoir. The characteristics and properties of these fracture networks play an important role in shale-gas reservoir performance. Therefore, it is significantly important to accurately characterize and represent these features in the modeling of these formations. This is often done by means of highly detailed discrete fracture representations that recent measurements and modeling techniques are able to generate. In Discrete Fracture and Matrix (DFM) modeling approaches, each fracture is modeled explicitly using, in most of the cases, highly resolved unstructured grids [25]. This allowed the simulation of fine-scale geological models with complex and various fracture geometries. For these reasons, DFM is considered as the most accurate representation of fracture networks but with the disadvantage of high computational cost as a large number of grid cells are generally involved; representing an obstacle for the application of DFM models at the field scale, therefore the necessity of upscaling. Hence, our objective in this study is to provide a coarse-scale representation of flow in shale systems without losing the important characteristics related to the structural features and transient effects observed in shale-gas reservoirs.
Conventionally, the modeling of fractured reservoirs is handled by means of multiple continua approaches such as the dual porosity method. The idea, proposed by Barenblatt and Zheltov [3] and later introduced to the oil industry by Warren and Root [45], is founded on the subdivision of the system into two separate continua, the fracture network and the matrix, and to model the exchange between these two media using a transfer function, also called a shape factor. These multiple continua representations were the basic foundation for many of the upscaling models for fractured reservoirs, where equivalent permeability for the coarse block is determined based on either homogenization approach [1,36] or local single phase flow simulations over the fine-scale model [5,14].
Although these techniques are quite practical, they rely on a simplified representation of a fracture network and often limit the representation of the geological and flow complexity by a single parameterthe shape factor. Lee et al. [30] proposed a hierarchical fracture upscaling that is meant to reduce the error brought by homogenization when fracture length scale distribution is non-uniform or the network is poorly connected. In their approach, large-scale fractures were modeled explicitly and the effective permeability contribution from smaller fractures was determined analytically. The EDFM model in Li and Lee [31] can also be used for upscaling of fracture networks. This method represents an interesting development to the dual continuum approaches, where fracture networks are discretely connected to matrix blocks by a series of source terms. It also has the advantage that fractures and matrix are represented by independent grid domains. However, difficulties still might arise when evaluating the transfer function to describe the exchange of fluid between matrix and fracture.
Gong et al. [15] and Karimi-Fard et al. [27] presented a systematic multi-subregion (MSR) upscaling approach based on the integration of DFM into a general multiple continuum representation. The method was developed as an effort to include spatial variability within a local matrix region, considering that most of the dual-porosity implementations model the pressure and saturation as constant within the matrix. They succeed to resolve spatial variation within the matrix with a novel flow-based sub-griding technique using the solution of a local discrete fracture flow problem over each coarse grid-block. Application of the method to simulate 2D and 3D fracture models, with viscous, bouyancy and capillary forces is also shown in their work.
In our study, we investigate the upscaling of shale-gas reservoir with the objective of accurately reproducing important characteristics of flow in nano-scale porosity, and to honor geological descriptions of the fracture networks in full details. We introduce a high fidelity model that serves as the base case for our upscaling approach. In the high fidelity model, we implement the formulations that describe the physics of gas transport in the shale matrix and accurately describe the interactions with the structural features in stimulated shale reservoirs. This is done by resolving the fractures of various scales and geometries using a general unstructured grid, and then solving for the flow equations using the finite volume DFM approach. The coarse representation is then constructed by using agglomeration to either structured (EDFM upscaling) or unstructured grids (MSR upscaling). The transmissibilities between control volumes of upscaled models are extracted from the solution of the underlying fine-scale DFM model. In our coarse models, we focus on accounting of the strong transient effects induced by the shale formations. This entails an adjustment of the effective transmissibilities in coarse models by relating them to pressure.
We implemented all approaches in the Automatic Differentiation General Purpose Research Simulator (ADGPRS) [48,43,42,47], which supports a generally unstructured representation of computational grid. In addition, the ADGPRS simulator also offers the flexibility in coupling the flow to geomechanics [12], to account for the influence of stresses on the flow and the conductivity of the fracture network. An effective upscaling of flow in shale gas systems in combination with geomechanic effects will be a focus of our future work.

Gas dynamics in shale
For simplicity, we restrict our approach to a single-phase flow and assume that the gas phase is represented as a pure gas component (methane) which is behaving ideally. To describe the physics of gas transport in the nano-pore scale, we introduce an idealized fluid flow model based on a recent study made by Lunati and Lee [33]. In this study, they proposed a conceptual bundle of dual tube model to statistically predict the performance of shale-gas reservoirs.
The study of flow and transport phenomenon in shale systems entails the analysis of conservation of mass. The general mass balance equation is used to describe the flow across a single tube. For an element (control volume) within the medium, it is defined as: The equation depicts the accumulation of mass, the mass flux through the system and the physico-chemical interactions, where ρ is the density, j is the total flux in shale matrix and q is the source term.
The total flux refers to the dominating physical phenomena. In the modeling process, this is an advective flux due to mean gas velocity and a diffusive flux due to the gradient of density, The mean gas velocity u ( ) is defined as proportional to pressure gradient similarly to Darcy's law as follow: where k is the absolute permeability tensor and μ is the viscosity of the gas. Next, we will briefly summarize the model proposed by Lunati and Lee [33] based on the elementary kinetic theory for an isothermal system [18].

Knudsen number
First we present the definition of the Knudsen number and the mean free path. Knudsen number, K N , is used in order to determine the appropriateness of the continuum model. It is a widely recognized dimensionless parameter, defined as the ratio of the gas mean free path l and the pore diameter d, where l the molecular mean free path is defined as the average distance the molecules of gas travels between two successive collisions with other molecules. Using kinetic elementary theory, and assuming a Maxwell-Boltzmann distribution of the velocity, the mean free path becomes: where m is the molecular mass and σ is the cross-sectional area of collision. For methane = = σ πδ nm 0.42 2 2 , using a methane molecular diameter =°δ A 3.8 . Here we assume that the Knudsen diffusion is dominated which may be affected by the sorption in some shales [4].

Viscosity and permeability
The viscosity defines the ability of intermolecular collisions to transfer momentum and based on the elementary kinetic theory of gases, it is defined as: 1 is the Boltzmann constant. Note that viscosity is only a function of temperature; hence it is considered as constant for isothermal processes.
Based on the Hagen-Poiseuille equation for viscous flow in a pipe, permeability in the longitudinal direction can be expressed as: where κ is a dimensionless constant that is related to the configuration of the flow-paths (1/32 for a circular tube and 1/12 for planar fractures), d is the diameter of a pore or the aperture of a plane fracture.
As pore size d becomes very small and comparable to the mean free pathl, the no-slip condition at the solid boundary is no more applicable and the equation has to be modified. Brown et al. [7] proposed a correction to account for the slippage effect as: is the tangential-momentum accommodation coefficient that indicates the fraction of molecules that are diffusively reflected by the wall. Assuming we have tubes of rough surfaces that reflects all molecules diffusively, then = σ 1 v and substituting with the definition of mean free path in Eq. (5), the equation simplifies to:

Effective diffusivity
The molecular self-diffusion coefficient describes the mass transfer due to molecular collisions, and therefore is proportional to the thermal velocity and to the mean free path: When the size of the pore d is comparable to or smaller than the mean free path of the gas molecules l, the interactions with the solid walls of the pores become more significant than the intermolecular interactions. The mass transfer due to the effects of a collision with the wall is described by the Knudsen diffusion coefficient, obtained by replacing l by d in Eq. (10): Knudsen diffusion is more likely to be prevailing at lower pressure, while molecular self-diffusion in nano-pores is dominant at a higher pressure as intermolecular collisions are more likely. To consider the contribution of these two mechanisms, the effective diffusion coefficient introduced by Lunati and Lee [31] is used: which describes their combined effect and tends to D m when

Adsorption and desorption
In shale-gas systems, the gas is stored as a compressed gas in pores but also as an adsorbed gas to the pore walls and as a soluble gas in solid organic materials [22]. As the pressure within the reservoir is depleted by production, the adsorbed gas gets released into the nanopores followed by the diffusion of the dissolved gas to the surface of the pores.
The volume of adsorbed gas in shale formations can be significant; up till now its significance to production is the subject of many studies and analyses for different types of shales (e.g. [46,44]), and its contribution is mostly accounted when the pressure has declined to very low values [17]. For simplicity, we will only include the contribution of free gas trapped in the rock pores into our simulations which makes the proposed model representative for an early production from shale-gas reservoirs.

Discrete fracture model
Discrete fracture models, in which the fractures are represented individually and constrained by geological information, are considered as one of the most accurate techniques to model fracture networks. The use of unstructured grids allows the modeling of non-ideal fracture geometries, such as non-orthogonal and non-planar fracture orientations.

Different DFM approaches
Various procedures to solve the flow equations for systems using unstructured grids can be found in the literature. Using the finite element approach (FEM) such as the earliest work of Baca et al. [2] as they proposed to solve for 2D single-phase flow with heat and solute transport. Juanes et al. [24] presented a more general approach with FEM for 2D and 3D for single-phase flow in fractured porous media. These early methods were then extended to handle incompressible two-phase fluid flow including capillary pressure effect such as in the work of Kim and Deo [29] and Karimi-Fard and Firoozabadi [26], compositional multicomponent flow in Hoteit and Firoozabadi [20], and three-phase flow in Fu et al. [11]. Matthäi et al. [35] presented a control-volume finite-element (CVFE) approach to accurately quantify two-phase flow simulation in fractured rock masses using 3D hybrid meshes. Their work was then expanded to include applications for compressible threephase flow in Matthäi et al. [34] and in Geiger-Boschung et al. [13].
Finite element procedures are generally more computationally expensive than the standard finite-volume discretization, the latter being the most popular choice among the majority of existing reservoir simulation techniques. Karimi-Fard et al. [25] offered a simplified DFM model based on a finite volume approach. The method is applicable to discretization based on the connection list [32] and offers significant improvement in the efficiency of DFMs using unstructured grids.
In this work, we follow the approach proposed by Karimi-Fard et al. [25] where an unstructured control volume finite-difference technique with a two-point flux approximation is used. For a 2D problem, the matrix blocks are represented by polygons while fractures are represented by segments. The fracture thickness is not represented in the grid domain but only in the computational domain for flow rate evaluation, which consequently simplifies the griding of the fractured domain. The mean properties of the grid block as well as the evaluated variables are defined in nodes, at the centroid of each corresponding control volume which is representative of the entire grid block. In addition, the transmissibility at a fracture intersection is evaluated in a way that eliminates intermediate control volumes at the intersection. While for multiple fractures intersecting at a point, a star-delta connectivity transformation is used. This improves the numerical stability and time step size for the simulation.

Discretized equations
The mass balance equation for single phase, single component shale system after substituting with all the derived coefficients corresponding to each of the transport phenomena in shale gas can be described as: As we integrate the partial differential equation (13) over a finite control-volume, we introduce porosity (ϕ) into the flux term to upscale the equations from a cylindrical nanotube to nano-porous media. This allows the flux to take place across a fraction of the grid-block interface (ϕ A . ) equivalent to a stack of (n) tubes, where n represents a certain number of tubes and A ij is the interface between the grid blocks i and j.
The discretized flow term is then represented in terms of a list of connected control volumes, where the flow rate is related to the pressure gradient by a transmissibility term. For each control volume (V), we express the flow rate across each one of its interfaces with the neighboring cells as: which represents the flow rate from V i to V j . The cell to cell transmissibility presented here accounts for the coefficients corresponding to the advective flux and diffusive flux, T A and T D respectively. The term λ is the ratio ρ μ and along with the porosity ϕ, they are both computed using upstream weighting (upwind). The transmissibilities are defined at the interface using a harmonic average of the properties within the connected blocks, The half transmissibility, α i , is then defined for each of the flow regimes: • for advection . , while α j is computed in the exact same manner but using the corresponding subscript j. Both indexes are obtained from the connection list.
In equations (17) and (18),A ij is the area of the interface between the neighboring V i and V j , D i is the distance between the centroid of the interface and the centroid of V i , ⎯→ ⎯ n i is the unit normal to the interface inside V i and → f i is the unit vector along the direction of the line joining the control volume centroid to the centroid of the interface as illustrated in Fig. 1.
To evaluate the transmissibility for a connection between two fractures, an intermediate control volume V ( ) 0 is introduced at their intersection, Fig. 2. The purpose of this node is to redirect the flow and to allow for thickness variation, without the need to estimate variables at this location. This is presumed from the assumption that such intermediate control volume will typically have similar properties to the adjacent cells and has a much smaller volume.
As for systems made of more than two fractures intersecting at a point, the transmissibility computations follow an analogy to electrical circuits "star-delta" transformation for networks of resistors proposed in Karimi-Fard et al. [25]. Therefore, transmissibility for each connection pair of n intersecting fractures is computed by: Integrating Eq.(13) over each control volume and expressing flow       rates using eq.(15) provides the discrete form of the flow equations which are solved to obtain the pressure profile for the fully resolved model and for the determination of the upscaled model parameters.

Upscaling
In this section, we elaborate on the upscaling models that aims to approximately reproduce the same flow behavior of the fine-scale DFM model with a lower number of degrees of freedom (grid-blocks) for cases where transient effects are important. We propose two different methods to achieve that objective, in a procedure that is based on two distinct steps. The first step concerns the building of the coarse grid and the second step is to determine the upscaled transmissibility information from the solution of the DFM model as explained in the previous section.

EDFM upscaling
The upscaling approach developed here is an extension to the work of Li and Lee [31], where they propose an embedded discrete fracture modeling (EDFM) technique. Based on this approach, the matrix domain is discretized separately with a structured Cartesian grid, while the fractures intersecting these matrix blocks are discretized by the matrix cell boundaries, see Fig. 3. The coupling of these two domains was made through a connectivity index based on the concept of the wellbore productivity index (PI). In their work, Li and Lee prescribe an analytical approach to approximate this index.
In our study, we generate the coarse-model grid using Cartesian blocks for the matrix. However, we use a numerical approach (upscaling) for characterizing flow between the matrix and the fracture. The fractured domain is first resolved based on an unstructured grid, established using a standard Delaunay triangulation scheme Shewchuk [38]. The coarse grid is then defined over the detailed model using a conformal aglomeration of fine-scale degrees of freedom. The fully resolved DFM model is then run up to a steady-state condition. The numerical integration over the coarse-grid blocks is used to provide the transmissibilities for the coarse model.
In general terms, transmissibility expresses the block to block mass flow rate in terms of the difference in pressure between the two blocks: From the solution of the fine-scale DFM model, we have matrix pressure p i k and fracture pressure p j k corresponding to each block k and the flow rate between them Q i j k , . For pressures this is accomplished using a pore volume weighted average of the pressures within the fine-scale cells associated with coarse block k, see Fig. 4.
The mass flow rate Q i j , is determined from the sum of fluxes crossing the fractures interface within the coarse block k as shown in Fig. 5. Matrix-matrix and fracture-fracture transmissibility, on the other hand, are computed in the exact same way used for DFM with the only difference that we are using a Cartesian grid to represent the matrix. Note that for simplicity we are considering a homogeneous matrix and therefore we do not deal with permeability upscaling in our EDFM model.

MSR upscaling
The use of the structured grid for the coarse model has certainly a practical aspect. However, it masks some of the flow dynamics and features observed in the fine-scale flow solution. Therefore, we propose another upscaling technique that tackles the stated problem and also portrays the extended transient effects prevalent in tight formations more accurately. The proposed method suggests the use of a flow-based gridding technique to capture the spatial variability within the matrix in a similar way to the work of Karimi-Fard et al. [27].
Based on the fact that pressure variation inside the matrix behaves like a diffusion process, the matrix is divided into regions using isopressure curves obtained from the pressure solution of a discrete fracture model. As in the previous upscaling technique, we obtain the solution of the fine-scale model using the shale-gas model. Next, we use this pressure solution for the construction of the multiple-subregion model as shown in Fig. 6; the shapes of the iso-pressure curves are strongly dependent on the fracture geometries. Note that the pressure inside the fractures is not changing significantly due to their high permeability.
The number of subregions in the upscaled model depends on the pressure levels selected from the DFM solution. Once the coarse-grid geometry is defined, the next step entails the determination of the connectivity map that will be used for our coarse-scale simulation. Fig. 7 shows a sample of the upscaled domain, indicating that the setup of these regions implies a one-dimensional character which allows the representation of the connections using a linear sequence as shown.
Next, we can perform a simulation with the coarse-scale model using a similar finite volume scheme. Therefore we need to determine the upscaled transmissibility to describe the exchange of flux from one subregion to another (Matrix-Matrix) or into a fracture (Matrix-Fracture). The average pressure of each adjacent region and the flux across the interface is also captured from the fine-scale DFM solution.

Numerical examples
In this section, we present applications of the proposed upscaling approaches. The two methods are applied to a single phase compressible shale-gas flow in two dimensions. The general properties of the reservoir domain are summarized in Table 1. For simplicity, in all the cases, the matrix rock is considered to be isotropic and homogeneous. All fractures are fully open and considered to be identical to simplify the comparison.
The simulation scenario in the fully resolved model is obtained by drainage of the reservoir domain through a producing well perforated in the fracture network. Since the system is isolated by no-flow boundaries, the pressure of the domain will get depleted with time and, after a transient period, will reach a pseudo-steady state. Fig. 8 depicts an example of the matrix-fracture transmissibility over time for a typical coarse block, indicating the pseudo-steady state Darcy flow at which the upscaled transmissibility is recorded.
The use of pseudo-steady state transmissibility is generally acceptable for conventional Darcy problems, where permeability is assumed as an intrinsic porous medium property. However, the flow physics in ultra-low permeability formations imposes a non-linear character on the relationships between the flow rate and the pressure drop. Therefore, we introduce a modification to account for the prevailing transient effects by relating the upscaled transmissibility to pressure. This is done by simply fitting a linear relationship to the same transmissibility we have extracted. Here we consider transmissibility as a function of the upwinding pressure in each coarse block as shown in Fig. 9.
We demonstrate the capabilities of the upscaled models presented in this work first on a small-scale synthetic fracture network (100 × 100 m 2 ) to highlight the inaccuracy that results from using a pseudo-steady-state transmissibility in shale-gas problems. Then we represent the application of a realistic case on a fracture network obtained from the outcrop analysis of the Whitby Mudstone Formation (WMF) (5500 × 7000 m 2 ), the exposed counterpart of the Posidonia Shales buried in the Dutch subsurface and a possible target for unconventional gas [21].

EDFM upscaling
The fracture domain is resolved with an unstructured grid made of 5458 cells (5223 triangles for the matrix) that will be the base for our fine-scale solution. A coarse grid for the EDFM upscaled model is made up of 42 blocks (6 × 7).
To extract the matrix-fracture connectivity information, we run the DFM simulation on the fine-scale problem and we record the transmissibilities for each coarse block. In our first upscaling attempt, we used a constant transmissibility factor, which is representative of the late time pseudo-steady state regime. The obtained upscaled solution, Fig. 12, demonstrates a large error which indicates that results from using a steady-state transmissibility are not consistent due to transient effects introduced by the shalegas dynamics. Notice, that in our upscaled model based on EDFM, the transient effects described in Section 2 are present in the evaluation of matrix-matrix flux. However, they are not sufficient to represent the transient regime at a larger (reservoir) scale.
To improve the upscaling results with EDFM, we introduced the modification to account for the prevailing transient effects in low permeability shale-gas reservoirs. This is done by relating the extracted transmissibility to the upwinding pressure associated to each coarse block, Fig. 13. The results of an upscaled simulation with pressure-dependent transmissibility demonstrate a significant reduction in the error against the fine-scale results. This example indicates that the introduction of the shale-gas dynamic, described in Section 2 or in other models (e.g. [23]), is not sufficient to fully resolve the transient effect typical for shale systems at a large scale.

MSR upscaling
In this approach, we follow the upscaling workflow based on the Multiple-Sub-Region approach. The fine-scale model is run for a relatively short period of time (about 1000 days) in order to extract the regions. We use 4 pressure levels to generate 41 regions in total. Once the coarse model is set up, we proceed on the extraction of the transmissibilities to represent the matrix-matrix and matrix-fracture exchange by running the DFM model for a considerable amount of time until the pseudo-steady state is reached. Fig. 14 shows the results of the upscaling solution using the constant pseudo-steady-state transmissibility at two different time steps compared to the DFM solution. We also perform the upscaling using transmissibility taken as a function of the upwinded pressure and we compare it again to the fine-scale solution as in Fig. 15. Notice, that in our MSR upscaled model, the flow in the matrix is assumed to follow a Darcy law unlike the EDFM model discussed before. Fig. 10 represents the L2 errors, computed over time, for the different upscaling we performed. The graph highlights the reduction in error as we introduce the transmissibility modification to account for the shale-gas transients; it also indicates that the flow-based grids are more accurate in capturing these effects. In general, the MSR model is capable to accurately capture the shale-gas dynamics at coarse scale even without additional modifications for the transient regimes. That indicates the importance of capturing the small-scale heterogeneity patterns even without any modification of dynamic models at a large scale.

Large scale fracture network
We now consider an application of the upscaling models on the larger complex network. This model is discretized using 19,727 unstructured cells (15,972 triangles for the matrix). The coarse model for the EDFM upscaling was defined using 100 blocks (10 × 10 square blocks), while for the MSR upscaling it was defined using 106 regions. Similarly to the previous example, the defined workflow is followed to obtain the upscaled models that correctly accounts for the shale-gas dynamics. The simulation results are shown in Fig. 16 and 17 while the L2 error plots are shown in Fig. 11.
We notice that the MSR model is generally more accurate than the EDFM upscaling. The error plots in Fig. 11 indicates that in EDFM upscaling, the error is overall higher and more dispersed compared to the MSR upscaling. In MSR model, the error can be higher in certain regions, especially with smaller volume, as these regions are more sensitive to the linear fit imposed on the extracted transmissibilities. On the other hand, EDFM upscaling require less work to implement as only the matrix-fracture connectivity data is extracted from the fine-scale solution while matrix-matrix connectivity can be easily computed analytically for Cartesian grids. We also observe that when fitting a linear relationship to flow data from the fine-scale solution to compute the upscaled transmissibilities, a generally better fit is obtained when flow based coarse grids (MSR) are employed compared to the Cartesian blocks (EDFM).

Concluding remarks
In this work, shale-gas dynamics were implemented in Discrete Fracture and Matrix (DFM) model in order to capture the highly detailed geological features of shale formations. The approach employs a finite volume discretization on a generally unstructured grid, which can efficiently resolve complex fracturenetworks. We thoroughly investigated and proposed the application of two systematic upscaling techniques that honor detailed fracture characterizations and applicable to flow in shale-gas systems. Moreover, the proposed methods are linked to the results of accurate DFM simulation for the extraction of upscaling parameters.
The first upscaling approach can be seen as an adaptation of the EDFM method of Li and Lee [31], as the coarse scale is made of structured grid blocks and the fractures are embedded within. In this approach, the governing dynamic of shale-gas system is applied to both fine-scale DFM and upscaled EDFM models without any modification. This method has an advantage of being easier to apply as we propose a semi-analytical approach for defining the coarse-model transmissibility and only the matrix-fracture connections are upscaled from the DFM solution.
Next, we propose a second upscaling model, which uses a flowbased gridding technique and can be seen as an adaptation of the MSR upscaling technique [27,15] for shale-gas dynamics. The regions are obtained from the DFM flow simulation, which depicts the pressure diffusion character of shale gas under production. Then, we extract the upscaling parameters, in the case of matrix-matrix and matrix-fracture connections, similarly to the original approach. In this method, the flow Fig. 16. Comparison between EDFM upscaled solution (100 coarse block) for the large-scale fracture system using shale-gas formulation and a pressure dependent transmissibility, on the left side the results correspond to 10,000 days of production and on the right side the results for the 50,000 days simulation, (a) and (b) represents the upscaling results, (c) and (d) are the DFM averaged results.
in the upscaled blocks is following simple Darcy assumptions without introducing an additional shale-gas dynamic model in the coarse-scale MSR blocks. Using the flow based regions to form the upscaled blocks offers an improved accuracy in replicating the flowing profile of the DFM solution with a significant reduction in the number of grid cells.
Using numerical examples of naturally fractured media, we identified an additional transient phenomenon related to the numerical homogenization used in the construction of both coarse models. We have shown that to depict the transient character of low permeability shale gas, one needs to consider the pressure dependency in the parameters defining the mass transfer across the regions. For that, we correlated the transmissibility to the upwinding pressure and obtained more accurate results, which we validated against the fine-scale DFM solution. While gas flow in EDFM upscaling is following a modified shale-gas dynamic, it is still not enough to capture accurately the finescale flow regime. On the other hand, the MSR approach only employs the classic Darcy flow and still is capable to capture the fine-scale dynamic with better accuracy.
In this work, we used a global upscaling approach that is most convenient for shale-gas field simulation. The practical development of shale-gas fields can be divided into separate isolated blocks each of which drained by an individual well; as the low permeability of these formations limits to almost none the interference between wells. This aspect allows us to run global DFM flow simulation for each isolated block and can be considered an achievable objective for determining the upscaled parameters. The use of a global flow solution with no flow boundaries is the key optimization and advantage over many of the available upscaling techniques that rely on local flow problems with specified boundaries to obtain the upscaling parameters. Fig. 17. Comparison between MSR upscaled solution (106 region) for the large-scale fracture system using shale-gas formulation and a pressure dependent transmissibility, on the left side the results correspond to 10,000 days of production and on the right side the results for the 50,000 days simulation, (a) and (b) represents the upscaling results, (c) and (d) are the DFM averaged results.