Hanyi Wang
email: hanyi@utexas.edu

Discrete Fracture Networks Modeling of Shale Gas Production and Revisit Rate Transient Analysis in Heterogeneous Fractured Reservoirs*

Keywords: Discrete Fracture Network (DFN), Shale Gas Production, Rate Transient Analysis (RTA), Complex Fracture Networks, Naturally Fractured Reservoirs, Heterogeneous Reservoirs

Horizontal wells with multiple hydraulic fractures are necessary stimulation technique for economically developing tight and shale gas reservoirs. In such reservoirs, the conventional well-test techniques are not suitable because of ultralow formation permeability. Rate transient analysis (RTA) is the widely used tool for analyzing these reservoirs for the purpose of reserves estimation, hydraulic fracture stimulation optimization, and development planning. However, the conventional rate transient analysis is based on the models that were derived from idealistic assumptions for homogenous reservoirs. In this article, we first review the industry's common practice for rate transient analysis and discuss why the idealized conceptual model may not be adequate for analyzing production data from shale gas reservoirs. Then, a unified shale gas reservoir model based on Discrete Fracture Networks (DFN) is presented to investigate how each mechanism influences shale gas production and the corresponding rate transient behavior. It is found that shale gas production and rate transient behavior are significantly impacted reservoir heterogeneity, fracture networks, non-Darcy flow, gas adsorption and completion efficiency. Short earlytime linear flow with long transitional flow period is an indication of either existence of abundant complex fracture network or heterogeneous completion with unevenly distributed hydraulic fractures. Consider the nature of non-unique results of RTA, information from other independent sources is required to achieve a consistent and holistic interpretation.

Introduction

The matrix permeability of shales is generally in the range of nano-Darcy, an enormous conductive surface area is required between the wellbore and the shale matrix to attain commercial production rates. To achieve this surface area, massive multistage hydraulic fracture treatments are used to create complex networks of fractures connected to the well. The geometry, areal extent, conductivity, and typology of these propped/un-propped fracture networks, which dictate shale gas production rate and its decline trend [START_REF] Wang | What Factors Control Shale Gas Production and Production Decline Trend in Fractured Systems: A Comprehensive Analysis and Investigation[END_REF], are generally difficult to quantify. So it is a challenge to diagnose production behavior and evaluate completion efficiency in these reservoirs. In conventional reservoirs, pressure transient analysis (PTA) is commonly used to estimate reservoir properties and post-stimulation productivity, but it has limited application in tight and shale reservoirs because the shut-in period required to is often too long to be viable. Rate transient analysis (RTA) affords the longterm testing of wells without shutting them in and allows for the estimation of key reservoir properties, which are essential to obtain reliable production forecasts and reserves estimates and to improve field development strategies.

For production forecast, the decline curve analysis (DCA) is probably the most frequently used production forecasting tool for shale gas reservoirs due to its relative simplicity and speed. The common methods used to estimate oil and gas reserves rely on a set of empirical production decline curves based on the following hyperbolic function [START_REF] Arps | Analysis of Decline Curves[END_REF]:

𝑞 𝑡 = 𝑞 𝑖 (1 + 𝑏𝐷 𝑖 𝑡) -1 𝑏 (1)
where 𝑞 𝑡 is the production rate at time t, 𝑞 𝑖 is the initial production rate at time t = 0, 𝐷 𝑖 and 𝑏 are two constants (the former is the initial rate of decline in production and 𝑏 is the rate of change in 𝐷 𝑖 over time, which control the curvature of the decline trend).The Arps equation was designed for conventional reservoirs where the boundary-dominated flow is the norm. However, shale gas reservoirs are characterized by transient production behavior and in general, boundary-dominated flow only occurs in later times. The flaws in Arps model has led to the development of many new DCA models for predicting estimated EUR in the shale gas wells, such as the power law exponential model [START_REF] Ilk | Exponential vs. Hyperbolic Decline in Tight Gas Sands: Understanding the Origin and Implications for Reserve Estimates Using Arps' Decline Curves[END_REF], logistic growth analysis [START_REF] Clark | Production Forecasting with Logistic Growth Models[END_REF]) and Duong's model [START_REF] Duong | Rate-Decline Analysis for Fracture-Dominated Shale Reservoirs[END_REF], etc. Even though all these models were formulated differently, they are all empirical equations and lack of underlying support of physics, so the same production data may lead to different estimation of production decline trend, with different practice of tuning parameters. In addition, empirical models cannot be used to analyze what factors cause the shift of production decline curve in field cases with different practices, which makes it impossible to assign value to one design over another and equally impossible to optimize the treatment for whichever goal is sought, either acceleration of recovery or increase in reserves.

Similar to pressure transient analysis, rate transient analysis starts with the conceptual modeling of wellbore and fracture geometry, then identify the flow regimes by plotting production data on a diagnostic plot. Fig. 1 shows a typical evolution of macroscopic flow regimes during production from a horizontal wellbore with multiple transverse hydraulic fractures in a homogeneous reservoir. For low permeability gas formations, fracture flow capacity is normally large enough to be treated as infinite conductivity and hence bi-linear flow is normally absent, and transient linear flow is often the first flow regime we encounter. Depending on fracture spacing and matrix permeability, this period can last for months or even years. When the pressure disturbance generated by multiple hydraulic fractures start to interfere with each other, virtue no-flux boundaries start to emerge between fractures and isolate each fracture to only deplete fluid in its own compartmentalized domain. This period is often referred as quasi-steady state flow or boundary dominated flow (The mode of boundary dominated flow seen in conventional reservoirs results from the pressure transient investigating all of the surrounding no-flow boundaries in the system, this is unlikely to occur in shale gas reservoirs because the matrix permeability is too low to enable investigation of large areas. In this article, the term "boundary dominated flow" specifically refers to the interference between the adjacent hydraulic fractures when the production pulse reaches the no-flow boundary). In very late time when most recoverable hydrocarbons have been depleted inside each compartmentalized domain, fluid from the far field that beyond the penetration of hydraulic fracture starts to contribute to production, and infinite acting, linear flow ensues. If production time is long enough and without the interference of nearby wells, the pseudo-radial flow may finally emerge. For each flow regime, special plots that based on the assumption of the underlying dominating mechanism can be used to estimate reservoir parameters, such as the drainage area, effective fracture surface area, average permeability, etc. Once these key parameters are required, we can predict the production decline trend and assess the effectiveness of the completion and stimulation design. Because the infinite acting, linear flow and pseudo-radial flow regimes may only occur at the very end of production life, so early-time transient linear flow and boundary dominated flow provide the most valuable data to analyze.

Fig. 1 Top view of typical macroscopic flow regimes for hydraulic fractured horizontal wells during production

To differentiate macroscopic flow regimes, rate normalized pressure (RNP) is often used [START_REF] Economides | Petroleum production systems[END_REF]). For gas reservoirs, pressure and rate transient responses need to be analyzed in terms of pseudopressure, 𝑚(𝑃), which can be defined as (Al-Hussainy and Ramey 1966):

𝑚(𝑃) = 2 ∫ 𝑃 𝜇 𝑔 𝑍 𝑃 𝑃 𝑟𝑒𝑓 𝑑𝑃 (2)
where 𝑃 𝑟𝑒𝑓 is some arbitrary reference pressure, 𝜇 𝑔 is gas viscosity and Z is Gas deviation factor. The rate normalized pressure and its derivative are computed as:

RNP = 𝑚(𝑃 𝑖 ) -𝑚(𝑃 𝑤𝑓 ) 𝑞(𝑡) (3) RNP′ = d RNP d ln (𝑡 𝑒 ) (4)
where 𝑃 𝑖 is the initial reservoir pressure, 𝑃 𝑤𝑓 is the wellbore flowing pressure, 𝑞(𝑡) is the production rate and 𝑡 𝑒 is the material balance time, which is calculated with the cumulative production 𝑄(𝑡) as:

𝑡 𝑒 = 𝑄(𝑡) 𝑞(𝑡) (5)
Plot RNP and RNP′ data on a log-log plot against material balance time, we can identify the flow regimes based on the slope of RNP′ (e.g., a half-slope indicates linear flow, a unit slope designates boundary dominated flow and a zero slope reveals radial flow). Once the flow regimes are identified, one can use specialized plots to estimate key reservoir parameters. For instance, the square root of time plot, RNP versus √𝑡, is probably the single most important plot to analyze data from ethe arlytime transient linear flow. Based on the early-time transient linear flow solution with the assumption of wellbore constant pressure, the RNP and √𝑡 follows a linear relationship (Wattenbarer et al. 1998;[START_REF] El-Banbi | Analysis of Linear Flow in Gas Flow Production[END_REF]:

RNP = 1 𝐴 𝑓 √𝑘 40.925 𝑇 √ 𝜙 𝑚 𝜇 𝑔 𝑐 𝑡 √𝑡 (6)
where 𝐴 𝑓 is the total fracture surface area and 𝑘 in formation permeability. The slope of the linear portion of RNP versus √𝑡 data can be used to estimate 𝐴 𝑓 √𝑘. If boundary dominated flow can be clearly observed right after early-time transient linear flow, then the termination time of transient linear flow can be used to estimate the distance of investigation (DOI), therefore, the fracture spacing, and hydrocarbon pore volumes (HCPV) can be determined based on simple volumetric calculations [START_REF] Anderson | Analysis of production data from fractured shale gas wells[END_REF].

Conventionally, RTA is based on the common assumptions that the reservoir is homogeneous and hydraulic fractures are uniformly placed along the horizontal wellbore. This may not be the case if the fracture spacing design is not optimized (e.g., a large number of tightly spaced perforation clusters are simultaneously propagating fractures), the stress interference may prohibit some fractures from growing [START_REF] Shin | Factors Controlling the Simultaneous Propagation of Multiple Competing Fractures in a Horizontal Well[END_REF] and promote some fractures to coalescence [START_REF] Wang | Numerical Investigation of Fracture Spacing and Sequencing Effects on Multiple Hydraulic Fracture Interference and Coalescence in Brittle and Ductile Reservoir Rocks[END_REF]). Field study [START_REF] Minner | Rose Field: Surface Tilt Mapping Shows Complex Fracture Growth in 2500' Laterals Completed with Uncemented Liners[END_REF]) also indicates that 80% fracture volume created at the heel and toe of a horizontal well and only 20% fracture volume created at the mid-lateral, due to poor stimulation design. If rock property is heterogeneous along the horizontal wellbore, fractures may only grow in brittle sections and leave the ductile sections unstimulated, because lower pressure is needed to initiate and propagate in brittle rocks [START_REF] Wang | Poroelastic and Poroplastic Modeling of Hydraulic Fracturing in Brittle and Ductile Formations[END_REF]. [START_REF] Wu | Paper SPE 184861 presented at the SPE Hydraulic Fracturing Technology Conference and Exhibition[END_REF] conducted numerical simulations of proppant distribution among multiple perforation clusters in a horizontal well, their study shows that proppant concentration in the toe-side clusters can be several times higher than the injected concentration, which increases the screenout risk of the toe-side clusters. Recent field study [START_REF] Palisch | Far-Field Proppant Detection Using Electromagnetic Methods -Latest Field Results[END_REF]) using electromagnetic methods to detect electrically conductive proppants distribution reveals that the proppants are unevenly distributed among each cluster and the effective hydraulic fracture length/propped fracture surface area can differ from cluster to cluster in a given stage. Field study [START_REF] Raterman | Sampling a Stimulated Rock Volume: An Eagle Ford Example[END_REF]) by examining core samples through the upper section of fractures along a horizontal wellbore indicates that fractures are not evenly distributed spatially and thus reservoir drainage may be non-uniform. So sub-optimal completion, reservoir heterogeneity and premature screen-out may lead to large reservoir volume unstimulated and in turn, impact the rate transient behavior.

In addition, gas can be stored as compressed fluid inside the pores or it can be adsorbed by the solid matrix in shale formations. The amount of adsorbed gas varies from 35-58% (Barnett Shale, USA) up to 60-85% (Lewis Shale, USA) of total gas initial in-place [START_REF] Darishchev | On Simulation of Flow in Tight and Shale Gas Reservoirs[END_REF]. Fig. 2 shows laboratory measurement of gas adsorption capacity from a shale sample at different pressures and temperatures. It can be observed that adsorbed gas can be released by decreasing reservoir pressure or increasing formation temperature. Without the aid of non-isothermal stimulation to alter gas desorption behavior (Wang et al.2014), how much desorption gas can be produced depends on the dynamics of pressure transient behavior in the reservoir. The coupled nature between pressure and gas desorption certainly can affect shale gas production and rate transient behavior.

Fig. 2 Gas adsorption capacity from a shale sample at different temperatures (Yue et al. 2015)

Darcy's law, which models pressure driven viscous flow, is another inherent assumption that underpins typical RTA. However, in shale reservoirs which have pore throat radii in the range of nanometers, Darcy's law may break down and gas molecules follow a somewhat random path while still maintaining a general flow direction governed by the pressure gradient. [START_REF] Wang | Impact of Shale Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas-Slippage, Desorption and Geomechanics[END_REF] presented a comprehensive literature review on gas transport in nanopores and the evolution of matrix permeability during pressure depletion and proposed a unified matrix permeability model which incorporates the mechanisms of non-Darcy flow/Gas-Slippage, the release of the adsorption gas layer and geomechanical effects into a coherent global model, as shown in Fig. 3. Besides pressure-dependent matrix apparent permeability, the existence of fracture networks can further complicate rate transient behavior. Now it is a well-known fact that most brittle shale formations have massive pre-existing natural fractures that are generally sealed by precipitated materials weakly bonded with mineralization. Such poorly sealed natural fractures can interact heavily with the hydraulic fractures during the injection treatments, serving as preferential paths for the growth of complex fracture networks [START_REF] Cipolla | A Practical Guide to Interpreting Microseismic Measurements[END_REF][START_REF] Zakhour | Real-Time Use of Microseismic Monitoring for Horizontal Completion Optimization Across a Major Fault in the Eagle Ford Formation[END_REF]. In addition, the permeability/conductivity of the fractures is a function of confining stress, fluid and proppant type for a specific reservoir [START_REF] Fredd | Experimental Study of Hydraulic Fracture Conductivity Demonstrates the Benefits of Using Proppants[END_REF][START_REF] Ghassemi | Sustaining Fracture Area and Conductivity of Gas Shale Reservoirs for Enhancing Long-Term Production and Recovery[END_REF], and fracture surface asperities/roughness play a critical role in maintaining residual fracture width of un-propped fractures [START_REF] Wang | A New Variable Compliance Method for Estimating In-Situ Stress and Leak-off from DFIT data[END_REF]2018). All these coupled mechanisms indeed pose a challenge for accurate interpretation of production data and reservoir characterization.

To overcome the limitations of conventional RTA of shale gas reservoirs, many authors have proposed some modifications for data interpretation. Nobakht and Clarkson (2012) studied the linear flow for constant flowing pressure, they found that the square root of time plot does not account for changing gas viscosity and compressibility, which leads to overestimation of 𝐴 𝑓 √𝑘. To address this issue, they proposed to use a corrected pseudo-time in which the gas viscosity and gas compressibility are evaluated at the average pressure in the region of influence [START_REF] Al-Hussainy | Application of real gas flow theory to well testing and deliverability forecasting[END_REF]. [START_REF] Clarkson | Production Analysis of Tight-Gas and Shale-Gas Reservoirs Using the Dynamic-Slippage Concept[END_REF] modified the pseudo-variables used for analyzing gas reservoirs in production data analysis to account for slippage/non-Darcy flow. They demonstrated that if the effect of slippage is not considered, it leads to noticeable errors in reservoir characterization. Nobakht et al. (2012) proposed a methodology to analyze the production data from a fractured well in a shale gas reservoir producing under a constant flowing pressure in the presence of desorption and slippage. This method uses a new pseudo-time definition that leads to a better estimation of 𝐴 𝑓 √𝑘. [START_REF] Clarkson | Incorporating Geomechanical and Dynamic Hydraulic-Fracture-Property Changes Into Rate-Transient Analysis: Example From the Haynesville Shale[END_REF] investigated the impact of both stress-dependent matrix permeability and fracture-conductivity changes on rate-transient signatures by use of modified pseudo-pressure and pseudo-time, and the fracture-conductivity changes are approximated by applying a time-dependent skin effect. They demonstrated that uncorrected data appear to have a very short early-time transient linear flow period, followed by apparent boundary dominated flow. [START_REF] Behmanesh | Impact of Distance-of-Investigation Calculations on Rate-Transient Analysis of Unconventional Gas and Light-Oil Reservoirs: New Formulations for Linear Flow[END_REF] proposed a new method to estimate DOI that can be used to calculate average pressure and evaluate the corrected pseudo-time. Their work shows the improved DOI and pseudo-variable calculations increases accurathe cy of linear-flow parameter estimation for unconventional reservoirs. Despite significant efforts have been made to improve RTA for unconventional reservoirs, currently, only early-time transient linear flow has been thoroughly investigated. In addition, all previous studies assume uniform fracturing spacing and homogenous stimulated reservoir volume (SRV) properties. However, in reality, gas flow in shale reservoirs, differs significantly from that in conventional reservoirs, is subject to more nonlinear, coupled mechanisms, including adsorption/desorption, non-Darcy flow, and strong rock-fluid interaction, and rock deformation within nano-pores, coexisting with complex fractures and multi-scaled heterogeneity. How these fully coupled physics will impact our identification of flow regimes and the corresponding RTA have not been investigated in a rigorous and systematic manner and no literature has attempted to do so. In this study, we are going to present a discrete fracture network (DFN) model which incorporates real gas transport, nano-flow mechanisms and geomechanics into fractured shale systems. This model is used to predict shale gas production and reservoir depletion under different scenarios and investigate how each mechanism influences rate transient behavior. The reservoir heterogeneity presented in this study refers to the heterogeneous properties on the scale of natural fractures or hydraulic fractures, not refers to grain size heterogeneity. The intention of this paper is not to seek a general formulation/methodology that can accurately analyze and history match production data under all circumstances, but rather to enhance our understanding of reservoir depletion and rate transient behavior in heterogeneous fractured shale gas reservoirs. More in-depth knowledge regarding the effects of factors controlling the rate transient behavior can help us develop more reliable models to better characterize reservoir properties, forecast shale gas production trend, optimize completion strategy and conduct uncertainty analysis.

Fully Coupled Discrete Fracture Network Modeling and Base Case Simulation

Darcy's law with apparent permeability correction is used to model the gas flow rate within the shale matrix:

𝒒 𝑔 = - 𝒌 𝑎 𝜇 𝑔 • 𝛻𝑃 (7)
where 𝒒 𝑔 is the gas velocity vector, 𝑃 is reservoir pressure and 𝒌 𝑎 is the matrix apparent permeability tensor, which is pressure and temperature dependent. The continuity equation within shale gas formation can be written as:

𝜕𝑚 𝜕𝑡 + ∇ • (𝜌 𝑔 𝒒 𝑔 ) = 𝑄 𝑚 (8)
where 𝑚 is the gas mass content per unit volume, 𝜌 𝑔 is gas density, 𝑄 𝑚 is the source term and 𝑡 is the generic time. The gas mass content 𝑚 is obtained from two contributions: 𝑚 = 𝜌 𝑔 𝜙 𝑚 + 𝑚 𝑎 (9)

𝜙 𝑚 is matrix porosity and 𝜌 𝑔 𝜙 𝑚 is the free gas mass in the shale pore space per unit volume of formation, while 𝑚 𝑎 is the adsorbed gas mass per unit volume of formation, which can be determined from the Langmuir isotherm [START_REF] Langmuir | The Constitution and Fundamental Properties of Solids and Liquids[END_REF]:

𝑚 𝑎 = 𝜌 𝑚 𝜌 𝑔𝑠𝑡 𝑉 𝐿 𝑃 𝑃 + 𝑃 𝐿 ( 10
)
where 𝜌 𝑚 is the shale matrix density, 𝜌 𝑔𝑠𝑡 is the gas density at standard conditions, 𝑉 𝐿 is the Langmuir volume and 𝑃 𝐿 is Langmuir pressure. Discrete Fracture Network with tangential derivatives can be used to define the flow along the interior boundary representing (natural and hydraulic) fractures within the porous medium. Flow behavior inside fracture is governed by:

𝒒 𝑓 = - 𝒌 𝑓 𝜇 𝑔 • 𝑑 𝑓 ∇ 𝑇 𝑃 (11)
where 𝒒 𝒇 is the gas volumetric flow rate vector per unit height in the fracture, 𝒌 𝒇 is the fracture permeability tensor, 𝑑 𝑓 is fracture width and ∇ T 𝑃 is the pressure gradient tangent to the fracture surface. The fracture permeability can be stressdependent, and their relationship can be obtained from laboratory experiment or empirical correlations, such as a power law relationship [START_REF] Cho | Pressure-Dependent Natural-Fracture Permeability in Shale and Its Effect on Shale-Gas Well Production[END_REF]:

𝑘 𝑓 = 𝑘 𝑓,𝑖 𝑒 -B𝜎 𝑚 (12)
where 𝑘 𝑓 is the fracture permeability, 𝑘 𝑓,𝑖 is the fracture permeability at initial reservoir conditions, B is a fracture compaction parameter that can be determined from experimental data, and 𝜎 𝑚 is the mean effective stress, which is the mean total stress minus pore pressure. The continuity equation along the fracture reflects the generic material balance within the fracture:

𝑑 𝑓 𝜕𝜙 𝑓 𝜌 𝑔 𝜕𝑡 + ∇ 𝑇 • (𝜌 𝑔 𝒒 𝑓 ) = 𝑄 𝑓 ( 13
)
where 𝜙 𝑓 is the fracture porosity, and 𝑄 𝑓 is the mass source term and it can be calculated by adding the mass flow rate per unit surface area from two fracture walls (left and right):

𝑄 𝑓 = 𝑄 𝑙𝑒𝑓𝑡 𝑓 + 𝑄 𝑟𝑖𝑔ℎ𝑡 𝑓 (14) 𝑄 𝑙𝑒𝑓𝑡 𝑓 = -𝜌 𝑔 𝒌 𝑎 𝜇 𝑔 𝜕𝑃 𝑙𝑒𝑓𝑡 𝜕𝒏 𝒍𝑒𝑓𝑡 (15) 𝑄 𝑟𝑖𝑔ℎ𝑡 𝑓 = -𝜌 𝑔 𝒌 𝑎 𝜇 𝑔 𝜕𝑃 𝑟𝑖𝑔ℎ𝑡 𝜕𝒏 𝒓𝒊𝒈𝒉𝒕 ( 16
)
where 𝒏 is the vector perpendicular to fracture surface. The in-situ gas density is calculated according to the real gas law:

𝜌 𝑔 = 𝑃𝑀 𝑍𝑅𝑇 (17)
T is the reservoir temperature, M is the gas average molecular weight and R is the universal gas constant. The Z-factor can be estimated by solving Equation of State (EOS) or using correlations for the gas mixtures. In this study, the Z-factor is calculated using an explicit correlation [START_REF] Mahmoud | Development of a New Correlation of Gas Compressibility Factor (Z-Factor) for High Pressure Gas Reservoirs[END_REF] based on the pseudo-reduced pressure (p pr ) and pseudo-reduced temperature (T pr ):

𝑍 = (0.702𝑒 -2.5𝑇 𝑝𝑟 )(𝑝 𝑝𝑟 2 ) -(5.524𝑒 -2.5𝑇 𝑝𝑟 )(𝑃 𝑝𝑟 ) + (0.044𝑇 𝑝𝑟 2 -0.164𝑇 𝑝𝑟 + 1.15) (18)
The advantage of using explicit correlation is to avoid solving higher order equations respect to Z-factor, which often leads to multiple solutions and increases computation efforts. Gas viscosity is an intrinsic property of gas itself that can be is affected by its compositions, pressure and temperature. Lee-Gonzalez-Eakin model (1966) is used to estimate the gas viscosity: The porous medium is assumed to be perfectly elastic so that no plastic deformation occurs. The constitutive equation can be expressed in terms of effective stress ( 𝜎 𝑖𝑗 , ), strain ( 𝜀 𝑖𝑗 ), and pore pressure (P):

𝜇 𝑔 = 10 -4 𝐾𝑒
𝜎 𝑖𝑗 = 2𝐺𝜀 𝑖𝑗 + 2𝐺𝜈 1 -2𝜈 𝜀 𝑘𝑘 𝛿 𝑖,𝑗 -𝛼𝑃𝛿 𝑖,𝑗 ( 23
)
where 𝐺 is the shear modulus, 𝜈 is the Poisson's ratio, 𝜀 𝑘𝑘 represents the volumetric strain, 𝛿 𝑖,𝑗 is the Kronecker delta defined as 1 for 𝑖 = 𝑗 and 0 for 𝑖 ≠ 𝑗, and 𝛼 is the Biot's effective stress coefficient. The strain-displacement relationship and equation of equilibrium are defined as:

𝜀 𝑖𝑗 = 1 2 (𝑢 𝑖,𝑗 + 𝑢 𝑗,𝑖 ) (24) 𝜎 𝑖𝑗,𝑗 + 𝐹 𝑖 = 0 ( 25
)
where 𝑢 𝑖 and 𝐹 𝑖 are the components of displacement and net body force in the i-direction. Combining Eq. ( 23)-( 25), we have a modified Navier equation in terms of displacement under a combination of applied stress and pore pressure variations:

𝐺∇𝑢 𝑖 + 𝐺 1 -2𝜈 𝑢 𝑗,𝑖𝑖 -𝛼𝑃𝛿 𝑖,𝑗 + 𝐹 𝑖 = 0 (26)
Darcy's law cannot describe the actual gas behavior and transport phenomena in nano-porous media. In such nanopore structure, fluid flow departs from the well-understood continuum regime, in favor of other mechanisms such as slip, transition, and free molecular conditions. The Knudsen number [START_REF] Knudsen | Die Gesetze der Molekularstromung und der inneren Reibungsstromung der Gase durch Roḧren (The laws of molecular and viscous flow of gases through tubes)[END_REF] is a dimensionless parameter that can be used to differentiate flow regimes in conduits at micro and nanoscale. For conduit with radius r, it can be estimated by

𝐾 𝑛 = 𝜇 𝑔 𝑍 𝑃𝑟 √ 𝜋𝑅𝑇 2𝑀 (27)
Table 1 shows how these different flow regimes, which correspond to specific flow mechanisms, can be classified by different ranges of 𝐾 𝑛 .

𝐾 𝑛 0 -10 -3 10 -3 -10 -1 10 -1 -10 The apparent permeability of shale matrix can be represented by the following general form [START_REF] Wang | Impact of Shale Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas-Slippage, Desorption and Geomechanics[END_REF]:

𝑘 𝑎 = 𝑘 ∞ 𝑓(𝐾 𝑛 ) (28)
where 𝑘 ∞ is the intrinsic permeability of the porous medium, which is defined as the permeability for a viscous, non-reacting ideal liquid, and it is only determined by the nano-pore structure of porous medium itself. 𝑓(𝐾 𝑛 ) is the correlation term that relates the matrix apparent permeability and intrinsic permeability using Knudsen number. Sakhaee-Pour and Bryant (2012) developed lab experiments based correlations for 𝑓(𝐾 𝑛 ). They proposed a first-order permeability model in the slip regime and a polynomial form for the permeability enhancement in the transition regime using regression method:

𝑓(𝐾 𝑛 ) = { 1 + 𝛼 1 𝐾 𝑛 Slip Regime 0.8453 + 5.4576𝐾 𝑛 + 0.1633𝐾 𝑛 2 Transition Regime 0.1 < 𝐾 𝑛 < 0.8 (29)
where 𝛼 1 is permeability enhancement coefficient in slip regime. To ensure the approximation of continuity of 𝑓(𝐾 𝑛 ) at the boundary region of slip and transition regime, where no existing model available, 𝛼 1 is set to be 4. Intuitively, Eqs. ( 27) and ( 29) predict net increases of 𝐾 𝑛 and 𝑓(𝐾 𝑛 ) with decreasing pore pressure. However, decreasing pore pressure can also lead to reduction in pore radii and, in turn, reduce the intrinsic permeability 𝑘 ∞ . It can be seen from Eq. ( 28) that the evolution of matrix apparent permeability 𝑘 𝑎 is determined by two competing mechanisms (i.e., the reduction in 𝑘 ∞ and increase in 𝑓(𝐾 𝑛 ) during production).

Wang and Marongiu-Porcu (2015) conducted a comprehensive literature review on gas flow behavior in shale nano-pore space and proposed a unified matrix apparent permeability model, which bridges the effects of geomechanics, non-Darcy flow and gas adsorption layer into a single mode, by considering the microstructure changes in nano-pore space. In a general porous media, the loss in cross-section area is equivalent to the loss of porosity,

𝜙 𝜙 0 = 𝑟 2 𝑟 0 2 (30)
where variables with a subscript 0 correspond to their value at the reference state, which can be laboratory or initial reservoir conditions. Laboratory measurements by [START_REF] Dong | Stress-Dependence of the Permeability and Porosity of Sandstone and Shale from TCDP Hole-A[END_REF] show that the relationship between porosity and stress follows a power law relationship, and can be expressed using the concept of mean effective stress 𝜎 𝑚 :

𝜙 = 𝜙 0 ( 𝜎 𝑚 𝜎 𝑚0 ) -𝐶 ∅ ( 31
)
𝐶 ∅ is a dimensionless material-specific compaction constant that can be determined by lab experiments. For shale samples, the values of 𝐶 ∅ range from 0.014 to 0.056 [START_REF] Dong | Stress-Dependence of the Permeability and Porosity of Sandstone and Shale from TCDP Hole-A[END_REF]. Combining Eq. ( 30) and Eq. ( 31) leads to the relationship between pore radius and local stress

𝑟 = 𝑟 0 ( 𝜎 𝑚 𝜎 𝑚0 ) -0.5𝐶 ∅ (32)
When the adsorption layer is considered, the thickness of the gas adsorption layer, δ, can be interpolated based on a Langmuir type functional relationship:

𝛿 = 𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 ( 33
)
where 𝑑 𝑚 is the average diameter of gas molecules residing on the pore surface. And the effective pore radius (32) can be modified as:

𝑟 = 𝑟 0 ( 𝜎 𝑚 𝜎 𝑚0 ) -0.5𝐶 ∅ -𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 (34)
The relationship between intrinsic permeability and pore radius can have the following relationship

𝑘 ∞ 𝑘 ∞0 = ( 𝑟 𝑟 0 ) 𝛽 ( 35
)
where 𝛽 is the coefficient that define the sensitivity of permeability to the changes of pore radius. Different shale formations may have different nano-pore structure typology, which leads to different values of 𝛽. In this study, 𝛽 equals 2 by assuming the overall all intrinsic permeability resembles fluid flow in a capillary tube (Beskok and Karniadakis 1999). Combining Eq. ( 34), Eq. ( 35) and Eq. ( 28), we have the final expression of matrix apparent permeability:

𝑘 𝑎 = 𝑘 ∞0 (𝑟 0 ( 𝜎 𝑚 𝜎 𝑚0 ) -0.5𝐶 ∅ -𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 ) 𝛽 𝑟 0 𝛽 𝑓(𝐾 𝑛 ) (36)
Based on the formulations introduced above, a fully coupled shale gas reservoir simulator with embedded discrete fracture networks is constructed. Newton-Raphson method [START_REF] Ypma | Historical development of the Newton-Raphson method[END_REF] and finite element analysis [START_REF] Zienkiewicz | The Finite Element Method[END_REF] are used to solve all the coupled equations numerically.

In order to investigate and differentiate the impact of various mechanisms on rate transient behavior and flow regime identification, we first start with a base case scenario where conventional assumptions for RTA are valid. This includes uniform hydraulic fracture spacing, homogenous reservoir properties, and does not account for gas desorption, matrix permeability evolution and pressure/stress dependent fracture conductivity. In addition, non-stimulated reservoir volume (Non-SRV) does not contribute to production. Fig. 4 shows the idealized conceptual model for the base case simulation. Even though it is possible to simulate a complete section of a horizontal well with multiple transverse fractures, it is more efficient to simulate a unit SRV/drainage volume and apply symmetric conditions along the boundaries to account for late time boundary dominated flow, and then the total gas production from a horizontal well can be determined by adding the contributions from each SRV or drainage volume unit. So in this article, the simulated production refers to the production from a single SRV or drainage volume unit, unless otherwise specified. Table 2 shows all the input parameters. Constant production pressure is applied to the wellbore and no-flux conditions are applied at SRV boundaries. 2. Input Parameters for Base Case Simulation Fig. 5 shows the pressure distribution after 50 and 500 days of production. As expected, the pressure gradient only varies in the direction that perpendicular to the hydraulic fracture, which drives linear flow from formation to hydraulic fracture. It can be also noticed that after 50 days of production, transient linear flow front edge has not reached the drainage boundaries, while after 500 days of production, boundary dominated flow already emerged, and the pressure at the drainage boundaries has dropped to 28.1 MPa from the initial value of 35 MPa. Fig. 6 shows the production data on log-log and the square root of time plots. It can be clearly observed that transient linear flow dominates early-time production with signatures of a half-slope of RNP′ on the log-log plot and a linear trend of RNP on the square root time plot. A closer observation of the log-log plot, one can realize that the early-time transient flow (slope of RNP ′ is 0.5) terminates around 400 days of material balance time, followed by a transitional flow period (slope of RNP ′ is between 0.5 and 1) with mixed transient linear and boundary dominated flow, fully boundary dominated flow (slope of RNP′ is 1) developed at around 650 days of material balance time and dictates the rest of production life. Once the flow regimes and its duration are identified, and knowing the slope of the linear portion on the square root of time plot, one can estimate matrix permeability, fracture surface area, drainage area, etc., using methods that proposed by numerous authors as discussed in the previous section, provided that the conceptual model of Fig. 4 is valid. It should be noted that the plots of RNP can be replaced with the plots of Pressure Normalized Rate (PNR), which is the reciprocal of RNP. Then in such cases, the magnitude of the slope should remain the same, but with opposite sign.

Factors and Mechanisms Impact Rate Transient Behavior

Impact of Gas Desorption

In the above base case simulation, adsorption gas is not included in the model. In this section, we set Langmuir volume, 𝑉 𝐿 , and Langmuir pressure, 𝑃 𝐿 , to be 0.0013 m 3 /kg and 10 MPa, respectively, such that around 40% of total gas in place exists as adsorption gas. Fig. 7 shows that more gas can be produced when adsorption gas is included in the model, this is because more total gas initial in place if adsorption gas is considered. What factors control shale gas production and its decline trend will not be extensively discussed in this article, readers can refer to [START_REF] Wang | Impact of Shale Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas-Slippage, Desorption and Geomechanics[END_REF], [START_REF] Wang | What Factors Control Shale Gas Production and Production Decline Trend in Fractured Systems: A Comprehensive Analysis and Investigation[END_REF] for more detailed analysis on shale gas production mechanisms. Transient Linear Flow Slope between ½ to 1: Transitional Flow gas, it seems that adsorption gas does not impact early time flow regime identification because linear transient flow is still apparent with half-slope in the first a few hundred days of production, however, the slope of RNP′ with adsorption gas converges to a unit in much later time, following an extremely long transitional flow period.

Fig. 7 Cumulative production with and without adsorption gas

Fig. 8 Log-log plot with and without adsorption gas

Fig. 9 shows RNP on the square root of time plot with and without the inclusion of adsorption gas. It can be observed that the adsorption gas does make a noticeable difference on this diagnostic plot, which leads to smaller slope and longer duration of early-time transient linear flow on RNP data. If Eq.( 6) is used without correction for adsorption gas, 𝐴 𝑓 √𝑘 can be overestimated. In addition, because the termination of transient linear flow occurs later, the distance of investigation and drainage volume can also be overestimated if the impact of adsorption gas is overlooked.

Fig. 9 Square root of time plot with and without adsorption gas

Impact of Matrix Apparent Permeability

The apparent permeability in the matrix of shale gas formations can be altered by non-Darcy flow, gas slippage, the release of adsorption layer and rock compaction during production (as illustrated in Fig. 3). In this section, we examine the impact of matrix permeability evolution on log-log and the square root of time plots. Here we assume the average matrix pore radius, 𝑟 0 , is 3 nm and the dimensionless compaction constant, 𝐶 ∅ , is 0.35. Fig. 10 and Fig. 11 show the values of Knudsen Number and the corresponding matrix apparent permeability in the drainage volume after 50 and 500 days of production. The simulation results indicate that Knudsen Number increases as pore pressure declines, and the value of Knudsen Number reveals that slip and transition microscopic flow regimes (refer to Table 1) dominates the entire simulated domain. It can be also observed that the matrix apparent permeability is highest (around 183 nd) in the low-pressure zone region that adjacent to the fracture, while remaining the same level as initial matrix permeability (50 nd) in the region that has not been tapped by pressure depletion. The local matrix apparent permeability will evolve as pressure transient front edge continues to propagate: starting from areas that close to the hydraulic fracture to the rest of the drainage volume. 12 shows the cumulative production with and without matrix permeability evolution. The results indicate that the evolution of matrix permeability during production can indeed make a difference in gas production, due to the enhanced local matrix apparent permeability. However, it should be mentioned that how much benefits can be gained from matrix permeability evolution depends on many factors, such as the average pore radius, the existence of natural fractures, etc [START_REF] Wang | Impact of Shale Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas-Slippage, Desorption and Geomechanics[END_REF]. So the matrix permeability evolution may not be as influential as it is in this case if average pore radius is larger or well-connected conductive fracture networks exist.

Fig. 12 Cumulative production with and without matrix permeability evolution

Fig. 13 shows RNP′on log-log plot with and without matrix permeability evolution. The early-time transient flow regime can be clearly identified for both cases with the half-slope. However, compared to the case of static matrix permeability, the slope of RNP′ with matrix permeability evolution converges to a unit much earlier with a shorter transitional period. This is because enhanced local matrix permeability during production advances the occurrence of boundary dominated flow. Fig. 14 shows RNP on the square root of time plot with and without matrix permeability evolution. It can be observed the slope of the linear portion of RNP data is smaller and the early-time transient linear flow ends sooner if local matrix permeability enhancement is considered during production. This is expected because the pressure depletion front edges reach the SRV boundary sooner when matrix permeability is higher. If Eq.( 6) is used to estimate 𝐴 𝑓 √𝑘, then 𝑘 is not the original matrix permeability, it represents the characteristic apparent permeability inside the drainage volume during the period of early-time transient linear flow.

Fig. 14 Square root of time plot with and without matrix permeability evolution

Impact of Fracture Networks

In this section, the role of natural fractures, their distribution and pressure/stress dependent conductivity on gas production and rate transient behavior will be examined. Because natural fractures are formed within geological time scale and the formation itself may have gone through multiple tectonic events, the distribution and orientation of the natural fractures do not necessarily relate to the current in-situ stresses. For the purpose of simplicity and easy interpret results, the distribution and orientation of natural fractures are randomly placed inside the SRV in this study. Both hydraulic fracture and natural fracture are modeled with DFN using Eq.( 11) to solve pressure along fracture directions, the only difference is they have different values of conductivity and sensitivity to pressure/stress. The advantage of using DFN is that fractures can be models with arbitrary length, orientations and conductivity. And more importantly, it overcomes the limitations of dual-porosity models [START_REF] Barenblatt | Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata][END_REF][START_REF] Warren | The behavior of naturally fractured reservoirs[END_REF] that cannot capture pressure and rate transient behavior of most fractured reservoirs with horizontal wells (Kuchuk and Biryukov 2014). To better discretize complex fracture geometry, unconstructed mesh is applied. In the first example, we assume the permeability of natural fracture, hydraulic fracture, and shale matrix remain constant during production. And the average width and permeability of natural fracture are 0.1 mm and 200 md, respectively. All the other input parameters are the same as the base. Fig. 15 shows the pressure distribution in the simulated SRV unit after 50 and 500 days of production. Compared to Fig. 5, the existence of well-connected, conductive fracture networks dramatically changed the depletion pattern, and much more reservoir volumes can be depleted (but not evenly depleted) during the same period. So there is no surprise that more gas can be produced with the benefits of larger productive fracture surface area from fracture networks, as shown in Fig. 16. Fig. 17 shows RNP′on the log-log plot with and without fracture networks. Compared to the base case scenario where only the hydraulic fracture provides productive fracture surface, it seems the early-time transient linear flow regime disappears when complex, conductive fracture networks exist, and transitional flow dominated early time production. However, a closer observation of Fig. 18, we can still identify a linear portion RPN data, during which the linear flow period is too short to be detectable on a log-log scale plot. This is reasonable, because the existence of abundant natural fractures not only disturbs the linear flow from matrix to the hydraulic fracture, but also compartmentalizes the drainage volume into smaller depletion unit with irregular typology. So once the linear flow from matrix to the hydraulic fracture/natural fractures starts to interfere with each other, the signature of early-time transient linear flow vanishes on diagnostic plots. In this case, the termination time of transient linear flow should not be used to infer fracture spacing or the size of SRV. Next, we examine the impact of matrix apparent permeability evolution on log-log and the square root of time plots, as we did in section 3.2, but now with the incorporation of fracture networks inside the SRV. Fig. 19 shows the local matrix apparent permeability inside the SRV after 50 and 500 days of production. Compared to Fig. 11, it can be observed that the enhanced matrix permeability region covers more reservoir volume because of more effective pressure depletion with the aid of fracture networks. However, if we compare Fig. 20 with Fig. 12, one can notice that the impact of matrix permeability evolution during production has much less impact on cumulative gas production. This is due to the fact that with abundant, conductive fracture networks, the general the role of matrix permeability itself is diminished, so does its evolution. Finally, we examine the impact of dynamic fracture permeability evolution during production. The conception of pressuredependent-permeability for fractures, can sometimes be referred to as pressure-dependent-conductivity, stress-dependentpermeability or stress-dependent-conductivity in literature. For consistency, the term "pressure-dependent-permeability (PDP)", will be used throughout this article, to represent the changes in fracture conductivity as local effective stress increases and fluid pressure declines. Because fracture conductivity is the product of fracture permeability and fracture width, shift fracture permeability has the same effect as shift fracture conductivity. In addition, the fracture storage effects have a negligible impact on gas production, because the driving forces of primary shale gas recovery come from the gas volume expansion and gas desorption with declining pressure. Fig. 23 shows the relationships between fluid pressure inside fracture and fracture permeability that will be implemented as pressure-dependent-permeability for fractures in the following cases. It assumes the fact that the conductivity of natural fracture is much smaller than that of hydraulic fracture and it is more sensitive to pressure/stress changes, due to lack of the support of proppants. Fig. 24 shows the cumulative production with and without pressure-dependent fracture permeability. As expected, less gas is produced if fractures exhibit pressure-dependent permeability, because the conductivity of both hydraulic fracture and natural fracture is impaired as effective stress increases during production.

Fig.24 Cumulative production with and without pressure-dependent fracture permeability

Fig. 25 shows RNP′on the log-log plot with and without pressure dependent fracture permeability. It can be observed that the pressure-dependent fracture permeability is more likely to impact the slope of RNP′ during early-time linear flow and transitional flow period. Compared to Fig. 17, we can realize that it is difficult to set apart the individual impact of pressuredependent fracture permeability, because it is the fracture networks that still dominate the slope of RNP′ data on a log-log plot. In summary, if abundant, conductive and well-connected fracture networks exist inside the SRV, then the period of early-time transient linear flow regime can be too short to be identified on a log-log plot (but it can be identified on the square root of time plot), and the transitional flow period (the slope of RNP′ is between 0.5 and 1) tends to last for months or even years. Field cases (Jorge 2016; 2017) also supported that the transitional flow period can be long lasting and is strongly affected by the heterogeneity of complex fracture networks. So very short early-time linear flow that ensued by a long period of transitional flow can be an indication that the production is dominated by the productive surface area of complex fracture networks. Even though the dynamic changes of local matrix apparent permeability and fracture conductivity can occur during production, their effects on rate transient behavior can be overshadowed by that of fracture networks, so it is extremely difficult, if not impossible, to quantify their individual influence from production and pressure data alone.

Impact of Non-Simulated Reservoir Volume (non-SRV)

Sub-optimal completion design (stress interference from multiple simultaneously propagating fractures), reservoir heterogeneity (alternating brittle and ductile rocks, non-uniform distributed weak planes and faults) and screen out at some perforations/clusters can lead to unevenly stimulated reservoir volume along the horizontal wellbore, as shown in Fig. 27. The reservoir volumes that are not overlapped by SRV and remain its original flow capacity (natural fracture systems are not stimulated), will be classified and investigated as non-SRV in the following case. Fracture Networks With Dynamic Fracture Permeability Fracture Networks With Constant Fracture Permeability Fig. 28 shows the pressure profiles in the simulated domain after 500 and 1500 days of production. It can be observed that after 500 days of production, most area within the SRV has been well depleted. But the pressure in the Non-SRV region still remains relatively high. After 1500 days of production, pressure depletion inside the non-SRV becomes the main drive mechanism to sustain the production.

Fig. 28 Pressure distribution after 500 and 1500 days of production inside the simulated drainage volume Fig. 29 shows cumulative production and RNP′on the log-log plot with and without the existence of non-SRV. Within the first 200 days of material balance time, the cumulative production and the RNP′ data are identical for both cases, because pressure depletion only occurs within the SRV during this period and the non-SRV has not "feel" the pressure disturbance yet. However, after 300 days of material balance time and once non-SRV starts to contribute to production, the slope of RNP′ data of both cases diverges and more gas can be produced is the drainage volume is not bounded. Most interestingly, when non-SRV and fracture networks coexist, even though early-time transient linear flow is not identifiable on a log-log scale, a period of half-slope temporarily occurs in the late times (500-1000 days of material balance time), which indicates pressure depletion is dominated by the macroscopic linear flow from the non-SRV region to the SRV. Finally, if production time is long enough (10000 days of material balance time), radial flow emerges (slope equals zero) with large size of non-SRV. Normally, radial flow is rarely observed in shale gas production wells because of long horizontal wellbore and extreme low matrix permeability. Even if radial flow develops at the very late period of a well's lifetime, the only extra information we can obtain is the permeability outside the SRV, similar to "composite reservoir radial flow [START_REF] Chu | A new model for a fractured well in a radial, composite reservoir[END_REF][START_REF] Spath | Well test solutions for vertically fractured injection wells[END_REF])", and it does not reflect the reservoir properties inside the SRV and completion efficiency. In addition, shale gas wells decline rapidly and most production comes from the first 5 years of production [START_REF] Baihly | Shale Gas Production Decline Trend Comparison Over Time and Basins Revisited[END_REF][START_REF] Wang | What Factors Control Shale Gas Production and Production Decline Trend in Fractured Systems: A Comprehensive Analysis and Investigation[END_REF]) and the production from late-time radial flow period contributes little to the cumulative production of a well. In other words, only linear flow and boundary dominated flow are useful as regard to production prediction and reservoir characterization (i.e., estimate A f √k, DOI and the size of SRV).

Impact of heterogeneous completion

As discussed in Fig. 27, hydraulic fractures can be of unequal length and spacing in reality, as opposed to ideally, conceptualized homogenous completion with fractures of equal length and spacing (shown in Fig. 4). Fig. 31 shows the pressure distribution after 50 and 500 days of production in the case of a heterogeneous completion, with unevenly distributed hydraulic fractures along a horizontal wellbore. From the results, we can observe that in the first 50 days of production, each hydraulic fracture mostly drains from its own adjacent region. However, after 500 days of production, it is clear that the drainage area of each hydraulic fracture has overlapped and started to interfere with each other. Unlike numerical modeling and simulation, field measurement at the surface can't distinguish how much gas is produced from which fracture or stage, so often the total production rate from an entire horizontal well is used to analyze the rate transient behavior. Fig. 33 shows the log-log plot with the simulation results from unevenly distributed hydraulic fractures (shown in Fig. 31). We can observe that the early-time transient linear flow period lasts about 300 days of material balance time, and then followed by a long-term transitional flow period. Fully boundary dominated flow only starts to emerge after 5 years of material balance time. The departure from early-time transient linear flow can also be identified on the square root of time plot, as shown in Fig. 34. Under this scenario, the duration of linear flow period can only be interpreted to infer fracture spacing between the closest hydraulic fractures (hydraulic fractures #2 and #3 in this case). Compared to Fig. 6, we can conclude that the existence of fracture networks or a heterogeneous completion with unevenly distributed hydraulic fractures shortens earlytime transient linear flow and stretches transitional flow period.

Model Validation with Field Example

In the previous section, sensitivity analysis of different mechanisms on gas production and rate transient behavior has been presented and discussed using our presented fully coupled model. In this section, we present a production history match for the Marcellus Shale gas well using our fully coupled model. This dataset has been previously studied by [START_REF] Eshkalak | An Integrated Reservoir Model for Unconventional Resources, Coupling Pressure DependentPhenomena[END_REF]. Similar to our presented model, [START_REF] Eshkalak | An Integrated Reservoir Model for Unconventional Resources, Coupling Pressure DependentPhenomena[END_REF] modeled the non-Darcy flow effects in the matrix with an exponential decline for the change in propped-hydraulic-fracture and induced-fracture permeability as a function of confining pressure. The choosing of matching parameters and fracture permeability estimation are discussed by [START_REF] Mirani | Production-Pressure-Drawdown Management for Fractured Horizontal Wells in Shale-Gas Formations[END_REF]. Fig. 35 shows the production history matches obtained with two different models. The result shows that our model yields a reasonable match of the production data. Note that that the early-time match of the model with the actual field data is relatively poor in the first 6 months. This is primarily related to the flow-back methodology used on the well, that involves high drawdown effort to unload fracturing fluids and resulting in two-phase flow. For example, 3D seismic data can be used to infer the density of natural fractures [START_REF] Jenkins | Quantifying and predicting naturally fractured reservoir behavior with continuous fracture models[END_REF]), but it needs to be heavily calibrated with a sufficient sample statistic of control data, such as core analysis, image logs, shear wave polarization, etc. [START_REF] Hunt | On calibrating curvature data to fracture density: causes[END_REF]. Microseismic data can be used to calibrate fracture networks and constrain fracture dimensions [START_REF] Clarkson | integration of rate-transient and microseismic analysis for unconventional gas reservoirs: Where reservoir engineering meets geophysics[END_REF][START_REF] Zhou | Integration of microseismic and well production data for fracture network calibration with an L-system and rate transient analysis[END_REF], individual stage fracture effectiveness can be assessed with production logging, distributed acoustic sensing (DAS) and distributed temperature sensing (DTS) data [START_REF] Pang | Production Performance Modeling of Shale Gas Wells with Non-Uniform Fractures Based on Production Logging[END_REF][START_REF] Wheaton | A Case Study of Completion Effectiveness in the Eagle Ford Shale Using DAS/DTS Observations and Hydraulic Fracture Modeling[END_REF]). The use of proppant tracer followed by spectral gamma ray logging can help identify the unstimulated area in a horizontal wellbore [START_REF] Leonard | Evaluating and Optimizing Refracs-What the Diagnostics Are Telling Us[END_REF]. Hydraulic fracture geometry and SRV permeability can also be correlated to offset well pressure data through poroelastic response [START_REF] Roussel | Introduction to Poroelastic Response Monitoring -Quantifying Hydraulic Fracture Geometry and SRV Permeability from Offset-Well Pressure Data[END_REF]. Pressure data from diagnostic fracture injection test (DFIT) can be used to determine representative matrix permeability, minimum in-situ stress, characterize the roughness of fracture walls and estimate the pressure/stress dependent conductivity of un-propped fractures [START_REF] Wang | A New Variable Compliance Method for Estimating In-Situ Stress and Leak-off from DFIT data[END_REF]2018). Flow back data can also be useful to infer effective fracture pore volume if fracture compliance can be estimated from DFIT [START_REF] Fu | Estimating Effective Fracture Pore Volume From Flowback Data and Evaluating Its Relationship to Design Parameters of Multistage-Fracture Completion[END_REF]. The data of stretched transitional flow period can possibly yield information on fracture complexity and heterogeneity (Jorge 2016; 2017). Despite significant efforts have been made trying to gauge the dimensions, typology of stimulated fractures and the size the stimulated reservoir volume based on the shift and progression of flow regimes using rate transient analysis, how to address the issue of extreme non-unique interpretation and calibrate rate transient analysis with other measurements and modeling still remain a challenge. This requires a good understanding of the accuracy, limits of various information sources and the linkage between them. Nevertheless, the early-time transient linear flow regime is least influenced by pressure interference, reservoir heterogeneity and completion efficiency, and thus, provides us the most reliable information to assess the total productive fracture surface area and matrix flow capacity.

Conclusions

Although advances in drilling and completion technology have enabled commercial production from ultra-low permeability reservoirs, such as shale and tight gas reservoirs, there remain significant challenges to the optimization of field development.

Provided sufficient data quality, rate transient analysis can be used to provide information about reservoir properties and completion efficiency, such as the size of stimulated reservoir volume, original gas in place, effective fracture surface area, matrix permeability, etc. However, current common practice for analyzing rate transient behavior still largely depends on the idealized conceptual model, with the assumptions of Darcy's flow, homogeneous reservoir, static fracture conductivity, matrix permeability, and uniform fracture spacing and length. But field and downhole measurements constantly indicate that this may not be the case in reality, as production contributions vary from stage to stage, and in many cases, a significant portion of perforation clusters or fracturing stages do not contribute to production at all. Previous analytical models that based on the idealized conceptual model only yield a homogeneous deterministic estimation and fail to capture the effects of heterogeneity. So it is imperative to have a better understanding of rate transient behavior in heterogeneous fractured reservoir and its implications when to evaluate completion efficiency and estimate reservoir properties. In this study, we presented a unified shale gas reservoir model, which incorporates real gas transport, nano-flow mechanisms and geomechanics into fractured shale systems. The solutions from this model enable us to investigate how previously overlooked mechanisms/scenarios could affect our interpretation of rate transient behavior. Conclusions reached from the analysis of this paper include the following:

1. Shale gas production and rate transient behavior are significantly impacted by fully coupled mechanisms, reservoir heterogeneity, fracture networks and completion efficiency. 2. The existence of adsorption gas elongates the duration of transitional flow period, where the slope of RNP′ is between 0.5 and 1 on a log-log plot. In addition, 𝐴 𝑓 √𝑘 and the size SRV can be overestimated if the impact of adsorption gas is neglected. 3. Matrix apparent permeability evolution during production advances the occurrence of boundary dominated flow and shortens the duration of transitional flow. If early-time transient linear flow period is used to estimate 𝐴 𝑓 √𝑘, then k should be interpreted as the characteristic apparent permeability inside the SRV during that period, other than taken as the original matrix permeability. 4. If abundant, conductive and well-connected fracture networks exist inside the SRV, then the period of early-time transient linear flow regime can be too short to be detectable on a log-log plot, and the ensued transitional period tends to be severely stretched. In this case, the termination time of transient linear flow should not be used to infer fracture spacing or the size of SRV, because of depletion interference between fracture networks. 5. When hydraulic and natural fractures exhibit pressure/stress dependent conductivity, it reduces the effective fracture surface area and leads to longer termination time for transient linear flow regime. 6. When SRV is surrounded by non-SRV, the half-slope of RNP ′ data can occur in the late times on a log-log plot, following a long transitional flow period. This late linear flow regime indicates pressure depletion is dominated by the macroscopic linear flow from the non-SRV region to the SRV. 7. Similar to the impact of fracture networks, a heterogeneous completion with unevenly distributed hydraulic fractures with varying length can shorten early-time transient linear flow and stretches transitional flow period. In such case, the duration of linear flow period can only be interpreted to infer fracture spacing between the closest hydraulic fractures. 8. Short early-time linear flow with long transitional flow period is an indication of either existence of abundant complex fracture network or heterogeneous completion with unevenly distributed hydraulic fractures that are productive. 9. Even though it is still a challenge to estimate fracture dimensions and the size of drainage volume with enough confidence based on the shift and progression of flow regimes due to the extreme non-unique interpretations, but early-time transient linear flow do offer us valuable information to infer the total productive fracture surface area and rock matrix flow capacity, because it is least influenced by reservoir heterogeneity and completion efficiency. = Matrix porosity at reference conditions

Nomenclature

Fig. 3

 3 Fig. 3 Mechanisms that alter shale matrix apparent permeability during production (Wang and Marongiu-Porcu 2015)

Fig. 4

 4 Fig. 4 Illustration of base case (idealized conceptual model) simulation domain with spatial discretization

Fig. 5

 5 Fig. 5 Pressure distribution after 50 and 500 days of production inside the drainage volume

Fig. 8 shows

 8 Fig.8shows RNP′on a log-log plot with and without the inclusion of adsorption gas. Compared to the case without adsorption

  Fig. 10 Knudsen Number in the drainage volume after 50 and 500 days of production

  Fig. 13 Log-log plot with and without matrix permeability evolution

  Fig. 15 Pressure distribution after 50 and 500 days of production inside SRV with fracture networks

FigFig. 18

 18 Fig. 17 Log-log plot with and without fracture networks

Fig. 19 FigFig. 21 and

 1921 Fig. 19 Matrix apparent permeability in the SRV after 50 and 500 days of production with fracture networks

  Fig. 23 Pressure-dependent fracture permeability

Fig. 26

 26 Fig.26shows RNP on the square root of time plot with and without pressure dependent fracture permeability. It indicates that the pressure dependent fracture permeability reduces the effective fracture surface area (smaller slope during early-time transient linear flow) and leads to longer termination time for linear flow regime.

Fig. 27

 27 Fig.27 Illustration of unevenly distributed, asymmetric hydraulic fractures that contribute to production

Fig. 29

 29 Fig. 29 Cumulative production and log-log plot with and without non-SRVFig.30shows RNP on the square root of time plot with and without the existence of non-SRV. It seems that the slope of RNP data is a straight line during the whole simulated production life when non-SRV is included. This is because after a transitional flow period, which ensued from a short early-time transient linear flow inside the SRV, the non-SRV regions start to

Fig. 31

 31 Fig.31 Pressure distribution after 50 and 500 days of production from unevenly distributed hydraulic fractures with varying lengthFig.32shows production rate and cumulative production for each hydraulic fracture. We can observe that even though each hydraulic fracture has different initial production rate, they all have very steep production decline trend, especially during the first year of production. Hydraulic fractures #2 and #4 have higher cumulative production because of their longer fracture length. Hydraulic fracture #4 has the highest cumulative production, despite it has the same fracture length as hydraulic fracture # 2. This is due to the fact that hydraulic fractures #2 suffers more severe production interference from adjacent fractures.

Fig. 33

 33 Fig. 33 Log-log plot with unevenly distributed hydraulic fractures

Fig. 35

 35 Fig. 35 Comparison of history matches with different models

Gas Adsorption Capacity (mg/g) Pressure (KPa)

  Their work indicates that despite rock matrix compaction, the apparent permeability in shale matrix continues to increase during production.

	1.4
	1.2
	1
	0.8
	0.6
	0.4		Experiments at 45	°C
	0.2		Experiments at 65 Experiments at 85	°C °C
	0		Experiments at 105 °C
	0	10000	20000	30000

Table 1 . Microscopic Fluid Flow Regimes Defined by Ranges of 𝑲 𝒏 (Roy et al. 2003)

 1

	1	> 10 1

  𝐴 𝑓 = Total productive fracture surface area (infinite conductivity), 𝑚 2 𝑏 = Hyperbolic exponent for Arps equation, 𝑚 2 𝐵 = Parameter for fracture pressure-dependent permeability, 1/𝑃𝑎 𝑐 𝑡 = Total formation compressibility, 1/𝑃𝑎 𝐶 ∅ = Material constant for pressure dependent porosity 𝑑 𝑓 = Fracture width, 𝑚 𝑑 𝑚 = Diameter of absorbed gas molecules, 𝑚 𝐷 𝑖 = Initial nominal decline rate for Arps equation, 1/𝑚𝑜𝑛𝑡ℎ 𝐸 = Young's modulus, 𝑃𝑎 𝑓(𝐾 𝑛 ) = Non-Darcy flow correction term 𝐹 𝑖 = Net body force along i-direction, 𝑃𝑎 𝐺 = Shear modulus, 𝑃𝑎 𝑘 = Matrix permeability, 𝑚 2 𝑘 𝑎 = Apparent permeability, 𝑚 2 𝒌 𝒂 = Apparent permeability tensor, 𝑚 2 𝑘 𝑓 = Fracture permeability, 𝑚 2 𝒌 𝒇 = Fracture permeability tensor, 𝑚 2 𝑘 𝑓,𝑖 = Fracture permeability at initial reservoir conditions, 𝑚 2 𝑘 ∞ = Matrix intrinsic permeability, 𝑚 2 𝑘 ∞0 = Matrix intrinsic permeability at reference conditions, 𝑚 2 𝐾 𝑛 = Knudsen number 𝑚 = Total gas content, 𝑘𝑔/𝑚 3 𝑚 𝑎𝑑 = Gas adsorption mass per unit volume, 𝑘𝑔/𝑚 3 𝑚(𝑃) = Pseudo-pressure, 𝑃𝑎/𝑠 𝑀 = Molecular weight, 𝑘𝑔/mol 𝒏 Gas production rate at generic time t, 𝑚 3 /𝑠 𝒒 𝒇 = Flow rate vector in the fracture per unit height, 𝑚 3 /𝑠/𝑚 𝒒 𝒈 = Velocity vector of gas phase, 𝑚/s 𝑄(𝑡) = Cumulative production at generic time t, 𝑚 3 𝑄 𝑓 = Mass source term in fracture, 𝑘𝑔/𝑚 3 /𝑠 𝑄 𝑚 = Mass source term in matrix, 𝑘𝑔/𝑚 3 /𝑠 𝑟 = Effective pore radius, 𝑚 𝑟 0 = Effective pore radius at reference conditions, 𝑚 𝑅 = Universal gas constant, 8.3145 𝐽/𝑚𝑜𝑙/𝐾 𝑅𝑁𝑃 = Rate normalized pressure 𝑅𝑁𝑃′ = Derivative of rate normalized pressure

		= Normal vector to fracture surface
	𝑃	= Reservoir pressure, 𝑃𝑎
	𝑃 𝑖	= Initial reservoir pressure, 𝑃𝑎
	𝑃 𝑤𝑓	= Wellbore flowing pressure, 𝑃𝑎
	𝑃 𝐿	= Langmuir pressure, 𝑃𝑎
	𝑃 𝑝𝑟	= Pseudo-reduced pressure
	𝑃 𝑟𝑒𝑓	= Reference pressure, 𝑃𝑎
	𝑞(𝑡) = 𝑡 𝑒 = Material balance time, 𝑠
	𝑇	= Reservoir temperature, 𝐾
	𝑇 𝑝𝑟	= Pseudo-reduced temperature
	𝑢 𝑖,𝑗	= Component of displacement, m
	𝑉 𝐿	= Langmuir volume, 𝑚 3 /𝑘𝑔
	𝑍	= Gas deviation factor
	𝛼	= Biot's coefficient
	𝛼 1	= permeability enhancement coefficient
	𝛽	= Coefficient for pore radius dependent intrinsic permeability
	𝛿	= Thickness of gas adsorption layer, 𝑚
	𝛿 𝑖,𝑗	= Kronecker delta
	𝜀 𝑖𝑗	= Elastic strain
	𝜀 𝑘𝑘	= Volumetric strain
	𝜇 𝑔	= Gas viscosity, Pa• s
	𝜈	= Poisson's ratio
	𝜌 𝑔	= Gas density, 𝑘𝑔/𝑚 3
	𝜌 𝑔𝑠𝑡	= Gas density at standard condition, 𝑘𝑔/𝑚 3
	𝜌 𝑚	= Matrix density, 𝑘𝑔/𝑚 3
	𝜎 𝑖𝑗	= Effective stress, 𝑃𝑎
	𝜎 𝑚	= Mean effective stress, 𝑃𝑎
	𝜎 𝑚0	= Mean effective stress at reference conditions, 𝑃𝑎
	𝜙 𝑓	= Fracture porosity
	𝜙 𝑚	= Matrix in situ porosity
	𝜙 𝑚0