Numerical investigation of Non-Darcy flow regime transitions in shale gas production

Shale gas reservoirs are organic rich formations with often ultra-low permeability. Gas is stored in free and adsorbed form. Conventional Darcy flow cannot fully describe the gas transport in such porous media. It is thus crucial to study the shale gas production considering different flow regimes and time dependent permeability, which can improve well-induced fracture design and ultimate gas recovery. In particular, this paper will focus on the transition in non-Darcy flow regimes near fracture-matrix interfaces using mathematical modelling. Especially, we investigate conditions at which these effects vanish, and Darcy flow assumptions become reasonable. The model describes a representative well-induced high permeability fracture surrounded by shale matrix. Investigated Non-Darcy mechanisms include apparent permeability, Knudsen diffusion, gas desorption and Forchheimer flow. Pressure depletion is the main driving force for single phase gas flow from the matrix to the fracture and from the fracture to the well. Pressure dependent gas desorption is defined by Langmuir isotherm and is a key production mechanism. This model is implemented in Matlab using Marcellus shale data. Scaling the model shows that recovery of gas depends on two dimensionless number that incorporates geometry relations, time scales of flow, intrinsic parameters of the porous media, non-Darcy constants, adsorption and boundary conditions. The dimensionless numbers define respectively if 1) the fracture or matrix limit the gas production rate 2) if non-Darcy flow is significant in the fracture or matrix. When one of the media limit production, the non-Darcy flow in the other medium has reduced importance and can be excluded from the model. Non-Darcy flow is important if it limits flow in the medium limiting the production. By checking the magnitude of the selected dimensionless numbers, the modelling approach can be determined in advance and significant computational time can be saved. The proposed model provides a tool for interpretation of complex shale gas production systems. It can be used for screening of flow regimes at different operational configurations and hence appropriate modelling approaches. The model can be used to optimise fracture network design and potentially in identifying stimulation operations that may significantly improve production rates and ultimate recovery from unconventional gas reservoirs.


Introduction
Technological advances in hydraulic stimulation of shale reservoirs have caused a fundamental shift to the exploration-and-production industry. These unconventional reservoirs typically have extremely low matrix permeability (10-100 nD, Cipolla et al., 2010) and exhibit gas stored both in free and adsorbed form. Gas flows from the nanopores in the matrix to the hydraulic fractures and then to the horizontal wells. This transport of gas comprises several flow mechanisms as investigated by a large number of scientists and engineers over many years (Bird, 2002;Javadpour et al., 2007;Javadpour, 2009;Civan, 2010;Civan et al., 2011;Beskok and Karniadakis, 1999;Blasingame, 2008;Moridis et al., 2010;Klewiah et al., 2019).
One of the key mechanisms is the non-Darcy flow; the traditional linear equation for flow in porous media based on Darcy's law is not sufficient for accurately describing high-rate flows. Non-Darcy flow occurs when inertial forces may no longer be neglected compared with viscous forces (Hagoort, 2004). That is very common near gas production wells or in the near-wellbore region, especially in fractures where local velocities can be very high. Bybee (2006) suggested that in hydraulic fracture stimulation, non-Darcy flow can have a major effect on reduction of a propped half-length to a considerably shorter "effective" half-length, thus lowering the productive capability of the well and overall reserves recovery. Moreover, flow-capacity can reduce by 5%-30% in low-rate wells due to non-Darcy effects (Bybee, 2006). To account for this nonlinear behavior, an inertial term called the Forchheimer term is added to Darcy's equation. Forchheimer (1901) gave the empirical Forchheimer equation to model gas flow more accurately at high flow rates (Mustapha et al., 2015;Li and Engler, 2001;Belhaj et al., 2003;Jones, 1987;Ling et al., 2013;Barree and Conway, 2005;Zeng and Zhao, 2008). Al-Rbeawi (2018) showed that non-Darcy flow has a significant effect on the pressure profile of unconventional gas reservoirs, especially at early production time. Luo and Tang (2015) through semianalytical modelling concluded that non-Darcy flow in the fracture mainly reduces the effective conductivity. This varying conductivity and non-Darcy flow in the fracture make the pressure curves deviate from the type curves. Several efforts have been made over the past 10 years to identify the effects of non-Darcy flow on overall gas production from shale reservoirs (Wang and Marongiu-Porcu, 2015;Fan et al., 2019;Al-Rbeawi, 2019;Luo and Tang, 2015;Pang et al., 2018;Wang et al., 2017). However, there appears to lack a clear understanding on exactly where the transition from Darcy to non-Darcy flow occurs, quantifying this transition and assessing how its importance can be estimated a priori. The objective of this work is to address these issues. As an analysis tool, we consider a 1Dþ1D combined fracture/matrix model that allows systematic evaluation of the role and magnitude of the different mechanisms. This extends work presented in Berawala et al. (2019) focusing on flow regime characterization to now also consider and focus on non-Darcy flow mechanisms. Similar to the works by Mainguy and Ulm (2001) and Andersen et al. (2014Andersen et al. ( , 2015 the model consists of a high-permeable fracture (length L y ) with width 2b. This depicts a typical hydraulic fracture in a real-field scenario. The fracture can have non-uniform width and is symmetrically surrounded by shale matrix of fixed length L x and low permeability as shown in Fig. 1. We implicitly assume equally spaced perforation intervals by assuming fixed matrix length. The gas is stored densely in the matrix by adsorption (modelled by a Langmuir isotherm), in addition to free gas in the pores. Apparent permeability is used to account for gas slippage effects, effective stress, adsorption and flow regimes relevant due to the nano-pore structure of the shale matrix. The pressure gradient towards the fracture and the well causes free gas from the matrix nanopores to flow. With pressure depletion, gas adsorbed onto the kerogen material desorbs into the pore space where it can flow as free gas to the fracture. A transfer term takes care of the communication between the fracture and the matrix. The system consists of a pressure-diffusion equation for the fracture which is coupled with a pressure-diffusion equation in the matrix. The model is scaled to derive dimensionless numbers that characterize the model. The Forchheimer term is incorporated into the flow equations using a correction factor f denoting flux reduction compared to Darcy flow. The scaled system gives rise to non-linear partial differential equations which are solved numerically using an operator splitting approach. Using the scaled model we then address the following questions of practical importance: -How can we quantify the flow transition from Darcy to non-Darcy? -What are the conditions under which non-Darcy effects in the matrix become significant for gas recovery? -How does non-Darcy flow affect flow regimes in shale gas production?
Marcellus field and literature data are used to parameterize the model. Sensitivity analysis is performed to see the effect on gas recovery with time and 2D distributions of scaled pressure, and the transition factor fðzÞ.

System geometry
Assume a hydraulic fracture extending perpendicularly from a well, along which we define the positive y-axis, starting from the well perforation. The fracture has width 2bðyÞ and it is assumed the width can vary from the perforation until its length L y . The fracture width is assumed to vary linearly and recover gas from the matrix only in its perpendicular direction (x-direction). The matrix is assumed to behave identically on both sides of the fracture. This results in a net no-flow boundary at x ¼ L x (the matrix half length). In the following we study the matrix and fracture only on the right-hand side of the symmetric system ð À b < x < L x Þ.

Mass conservation equation
Consider a domain with volume V containing gas in free and adsorbed form. The mass of gas in the porous media volume changes due to flow in and out of the interface ∂Ω with area A as expressed by the mass balance equation (LeVeque, 2002): where ϕ is porosity, ρ g ðp g Þ gas density, a g adsorbed gas (mass per solid volume), u Darcy mass flux vector, n is the unit normal vector pointing out of Ω; and p g is the gas pressure.

Fracture domain
In the fracture, gas adsorption is negligible, i.e., a f g ¼ 0: The fracture width, denoted 2bðyÞ, can vary with distance from the well. Considering a volume dV ¼ 2bðyÞ h dy→0 we get from (1): Since the fracture is surrounded by matrix symmetrically, the two source terms contribute identically: À

Fig. 1.
System geometry: the near well reservoir is seen from above where a fracture with variable width extends from a well perforation with length L y . Shale matrix surrounds the fracture on both sides with total length 2L x (typical perforation interval) . and we obtain:

Matrix domain
In the matrix, it is assumed that all flow is directed in the x-direction (towards the fracture), while flow in the y-direction is ignored. Considering a volume dV ¼ dx h dy→0 we get from (1):

Non-Darcy flow
Forchheimer's equation is defined by (Forchheimer, 1901): where γ is Forchheimer's constant. When γ ¼ 0 the formula reduces to This coefficient is usually obtained from experimental data. However, there are several correlations available in the literature to evaluate the Forchheimer's constant. In this paper, we use the correlation given by Tek et al. (1962): where C β is non-Darcy flow constant and k is the apparent permeability. It is useful that: where sð ⋅Þ denotes the sign function with value �1: We further have that: To get velocity in a form comparable with Darcy's law ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi The latter term is in the form: where which has a limit f→1 when z→0 and f→0 when z→ ∞: We name fðzÞ the transition factor as it denotes the transition of flow from Darcy to non-Darcy. In other words, the velocity is always less or equal to the velocity from Darcy's law. The value of fðzÞ ¼ 1 would indicate Darcy flow and fðzÞ < 1 would mean non-Darcy flow.

The volume factor and density relation
Introduce the inverse volume factor b g , using the real gas law: which implies: Assuming that the gas is ideal Z ¼ 1, we get: where b ' g is the inverse volume factor differentiated with respect to pressure. We note that b ' g is constant and has unit of inverse pressure.

Shale gas adsorption
From (14-16) we have: where we have defined: b a g ðpÞ ¼ ð1 À φÞ φρ gs b ' g a g ðpÞ: The pressure dependency is related through a Langmuir adsorption relation: where b a max is the max capacity of the shale (in units of pressure) to store gas on the surface and p L is the pressure at which half this capacity has been obtained.

Apparent permeability correction
Due to the nano-pore structure of the shale matrix, Darcy's law cannot describe the actual gas behavior and transport phenomena. Fluid flow departs from the continuum flow regime, in favour of other flow mechanisms such as slip flow, transition flow and free molecular conditions. The Knudsen number (Knudsen, 1909) which is a dimensionless parameter is used to differentiate between these flow regimes, for conduit with effective radius r e , it is defined as: (20) where T is absolute temperature, Z is gas compressibility factor, r e is effective radius of flow path, R is the universal gas constant and M is gas molecular weight. The apparent permeability of shale matrix can be represented by the following general form that relies only on Knudsen number K n ; and the effective intrinsic permeability k ∞e (Karniadakis and Beşk€ ok, 2001): k ¼ k ∞e f ðK n Þ (21) Florence et al. (2007) extended this derivation to characterize the non-Darcy gas flow in shale formations: where α K is the rarefraction parameter: Considering the effect of matrix compaction and adsorbed layer on the nanopore geometry, the effective intrinsic permeability is given as (Jiang and Yang, 2018): where r e is the effective radius of flow path and τ is the tortuosity of rock.
Huang and Ghassemi (2015) and Cao et al. (2016) gave generalized formulation that incorporates the overall contribution from effective stress, adsorption and flow regimes for the apparent gas permeability:

Summary of model
Substituting (6-25) into (4) and (5), we summarize the system for the flow of gas in the fracture-matrix system: These flow equations must be solved together with the initial and boundary conditions.

Initial conditions
Initially, the fracture and matrix have the same reservoir pressure p init . The adsorbed gas content in the matrix is defined from the isotherm at this pressure:

Boundary conditions
The well is perforated at y ¼ 0 with a known pressure: There is pressure and mass flux continuity across the fracture-matrix boundary. The fracture is closed (or has negligible production) from the matrix in y-direction. Similarly, the matrix has, due to symmetry, no flow at its outer boundary:

Scaling of the model
We now scale the system by introducing the following dimensionless variables: 2b 0 is the average width of the fracture. D is the gas diffusion coefficient resulting from the absolute pressure, scaled by the average pressure. The scaled coordinates obey 0 � X; Y � 1. The scaled pressure P ranges between 0 and 1 (corresponding to the well pressure in absolute terms). The apparent permeability is scaled using a reference permeability defined as: Accordingly, we scale the transition factor fðzÞ using the k ref and Δp: where: Applying these dimensionless variables to (32-37), we can define two time scales (extending those defined in Berawala et al. (2019) to include non-Darcy effects): � τ f representing diffusion of gas from the fracture to the well, � τ m representing diffusion of free and adsorbed gas from the matrix to the fracture, given as follows: Note that we have introduced the parameter ðG ' Þ ref , motivated as follows: Let the quantity of gas in free and adsorbed form be represented by G ¼ P þÂ g . A typical retardation factor is then: ðG ' Þ ref � 1 denotes both the factor of increased time to produce gas from the matrix due to adsorption, but also the increased quantity of gas released from adsorption during the pressure depletion. Further, we scale time with respect to the fracture diffusion time scale: After scaling, the coupled transport system (26) and (27) can be expressed in the following form: All the constant terms are collected into the two following dimensionless variables α and β: Physically, α represents the ratio of the time scales involved in gas diffusion from the fracture and gas diffusion from the matrix (including desorption), respectively. β denotes the pore volume ratio of the matrix relative to the fracture and consists of the additional amount of gas produced from the matrix due to adsorption, indicated by ðG ' Þ ref .
Comparing to Berawala et al. (2019), the updated model represented by (26-29) now consists of non-Darcy flow velocity u represented in a form comparable to Darcy flow through a transition factor f. It also accounts for permeability variation in the matrix due to gas slippage effects, adsorption and effective stresses. The scaling of model also leads to new definitions of dimensionless variables. However, it is important to note that diffusion time scale of matrix τ m comprises of reference transition factor f m ref whose value changes with matrix non-Darcy flow constant from case to case.

Simulation results
In this section, we study the behavior of the model (Eqs. (39)-(41)) by considering Marcellus shale Langmuir isotherm parameters defined in Table 2. We also perform sensitivity analysis to various input parameters to identify the conditions under which the non-Darcy effect becomes significant. In particular, we plot overall gas recovery vs time and show distributions of scaled pressure, transition factor and the relative amount of total mass in the system at 15% production of the mass initially in the reservoir, denoted by RF ob .
The system is solved by an operator-splitting approach, similar to that described by Andersen et al. (2014Andersen et al. ( , 2015, Berawala et al. (2019) and Berawala and Andersen, (2020). In this approach, we alternatively solve for flow in the y-direction (fracture diffusion) and flow in the x-direction (fracture/matrix diffusion and desorption). We refer to Appendix A for a detailed solution procedure. In total 600 cells were used with 20 equal cells in y-direction and 30 equal cells in positive x-direction. The numerical solution was validated by Berawala et al. (2019) by comparing the model with the established industry software Eclipse (GeoQuest, 2009) for a case with negligible amount of adsorbed gas in the matrix, a fracture with uniform width and gas transportation driven by Darcy flow. The full system (both sides of matrix surrounded by fracture) was then modelled using 80⋅20⋅1 ¼ 1600 blocks and solved fully implicit. Using gas recovery and pressure vs time profiles at different locations relative to the fracture and well, the numerical solution was found consistent with that of Eclipse.
Mean pore radius of shale matrix is usually in the range of 1-100 nm (Javadpour et al., 2007;Loucks et al., 2009;Zou et al., 2012;Yao et al., 2013), we use a representative value 14 nm for simulation. Using the input parameters defined in Table 2, we get Knudsen number (K n Þ in the range of 0.01-0.1 which indicates slip-flow or transition flow regime (refer Table 1). This implies the mean-free path of gas molecule is less or of the same magnitude as pore size of the matrix. In this regime, the gas transport is mainly governed by Knudsen diffusion and the conventional Darcy's law equation with no-slip boundary conditions cannot be applied. The permeability described by Eq. (25) takes into account this slippage effect. The corresponding apparent permeability for the reference case parameters is shown in Fig. 2(b).
Eq. (12) described in 2.3 denotes the transition factor f which is a function of non-Darcy flow velocity z. The transition factor will always be less or equal to 1. fðzÞ ¼ 1 indicates Darcy flow and a value less than 1 would indicate how significant the non-Darcy effect is. We plot the transition factor f against z as shown in Fig. 3. We see that the gas transport is governed by non-Darcy effects when z > 10 À 1 : To establish the conditions under which this can happen, we perform sensitivity analyses (3.1) on various input parameters such as non-Darcy flow constant (C β Þ, pore size (r m Þ, fracture permeability (k f Þ, and shape (b 0 Þ and size. Further, we also compare how each of these parameters affect gas recovery with/without non-Darcy flow in the matrix and fracture.

Reference case demonstration
Using the reference case parameters listed in Table 2, we present scaled average gas pressure and recovery profiles against time for four systematically varied non-Darcy flow constants C β , see Fig. 4. The same values are applied to fracture and matrix. The C β ¼ 0 case indicates that Darcy flow is considered both in matrix and fracture. As seen from Fig. 4, the recovery process goes much faster for Darcy flow C β ¼ 0 compared to the non-Darcy cases C β ¼ 1e À 6; 3e À 6 and 9e À 6m À 2.5 . This indicates that with increasing the magnitude of non-Darcy flow constant, the gas is produced at a much slower rate. To report how significant the non-Darcy effect for individual cases is, we report the reference transi-  Table 3. The reference values are calculated using the pressure gradient between initial reservoir and well pressure divided over the entire length of the matrix (refer to Eq. (35)). For C β ¼ 1e À 6 m À 2.5 , we get f f ref ¼ 0:11 and f m ref ¼ 0:87, which indicates that the flow is reduced by 89% and 13% compared to Darcy flow for the same pressure gradient in fracture and matrix respectively. Also, from Eq. (11), we see that for the same pressure gradient, non-Darcy flow will give lower velocities than Darcy. In general, both Darcy and non-Darcy Forchheimer models predict the same behavior at low velocity. But at high velocities like in fracture, non-Darcy models results in reduced velocities limiting the overall gas recovery.
To further understand as to why non-Darcy effects limits the production, we plot in Fig. 5 distributions of scaled pressure, total mass and transition factor after RF ob ¼ 15% for three values C β ¼ 0; 1e À 6 and 3e À 6 m À 2.5 . Scaled total mass is defined as the relative amount of gas currently in place to initial mass of gas in matrix and fracture (refer to Appendix C, Eq. (68)). If C β ¼ 0 (Darcy flow), we see that the gas entering the fracture is instantly produced to the well. A zero-scaled pressure in the fracture for this case and uniform pressure and total mass distribution surrounding the fracture is observed. The flow here is mainly governed by matrix and we are in a matrix-controlled regime. In the case with C β ¼ 1e À 6; 3e À 6 m À 2:5 (non-Darcy flow), the gas requires some time to leave the fracture and a significant pressure is observed in the fracture. This high-pressure gradient reduces the production rate from the matrix and leads to more non-uniform production around the fracture. We can thus say that non-Darcy effects leads to significant residence time in the fracture and controls the rate of recovery. The flow becomes more fracture-dominated or fracturecontrolled regime. Inspired by works on spontaneous imbibition (Rangel-German and Kovscek, 2002;Andersen et al., 2014); Berawala et al. (2019) showed that also the production of shale gas can be classified into matrix-or fracture-controlled. A similar approach is applied here with respect to the role of non-Darcy flow.

Table 1
Classification of gas-flow regimes according to Knudsen number, K n (Roy et al., 2003).

Role of individual fracture and matrix non-Darcy flow constants
In order to understand the importance of non-Darcy flow in the fracture and matrix domains, we consider two cases where we turn off the non-Darcy effect alternately in each domain and plot scaled pressure and gas recovery versus time. These two cases are compared against Darcy flow denoted by dashed line in Fig. 6 and against the reference case where we used same C β ¼ 1e À 6 m À 2.5 for both fracture and matrix.
For the case where non-Darcy flow is considered only in the matrix (C f β ¼ 0; C m β ¼ 1e À 6 m À 2.5 , green line), we obtain similar recovery as in the case of Darcy flow. This shows that for the input parameters mentioned in Table 2, non-Darcy effects in the matrix do not play a significant role. When C f β ¼ 1e À 6; C m β ¼ 0 m À 2.5 (orange line), we get much lower recovery compared to the Darcy flow. The gas transport from matrix to the well is then fracture dominated, i.e. the time scale of transporting gas through the fracture limits the gas production compared to producing the gas from the matrix. This also follows from the mathematical formulation of diffusion time scale for fracture (Eq. (38)). High C β in fracture and its intrinsic properties give very low f f ref resulting in higher diffusion time for gas in fracture compared to matrix.

Effect of pore size
In this section, we investigate the role of matrix pore size r m on gas recovery. We consider four pore radii r m ¼ 10; 20; 40; 80 nm and plot Knudsen number, apparent permeability in Fig. 7(a) & (b). As seen from the figure, increasing pore radii in the matrix, increases the Knudsen number. For higher pore radii, pore size becomes comparatively larger than the mean-free path of gas molecules and gas is mainly driven by viscous forces. Apparent permeability is proportionally linked to the Knudsen number (refer Eq. (25)). Increase in Knudsen number with increasing pore size gives higher apparent permeability as shown in Fig. 7(b). Fig. 7(c) shows the recovery profile when only Darcy flow is considered in the system. It can be seen that at higher pore radii, gas travels from matrix to the fracture at much faster rate due to high apparent permeability, which is then produced from the fracture instantaneously. However, when the same cases are plotted in presence of non-Darcy effects both in matrix and fracture (Fig. 7(d)), we see delay in production. As seen in 3.1.1, non-Darcy effect also here cause a shift from matrix-controlled flow regime towards fracture-controlled flow regime. A high-pressure gradient is created in the matrix due to which gas is not uniformly produced from fracture surroundings. The fracture then limits the flow of gas and thus; the gas is not instantaneously produced to the well. Interestingly, for cases with sufficiently low pore radii (r m < 40 nm), gas recovery seems to be less sensitive in the  Table 2.  presence of non-Darcy effects as compared to Darcy flow in the system. This is because non-Darcy effects incurs residence time in the fracture in addition to increase travel time for gas to flow from matrix to the fracture.

Effect of fracture permeability
In the following, we consider the role of fracture permeability k f by plotting gas recovery for the system with/without non-Darcy effects. The fracture permeability is varied between 1 mD and 1000 mD. As seen in Fig. 8(a), fracture permeability barely has any influence on gas recovery for Darcy flow. The effect is dominant only when k f is very low (~1mD). Low fracture permeability gives lower matrix-to-fracture pressure gradient, resulting in lower recovery rate compared to the  case with much higher recovery (1000 mD). This was also observed by Berawala et al. (2019). For non-Darcy flow ( Fig. 8(b)), we see that fracture permeability becomes very important and the production is limited by the fracture until k f e100 mD, a fracture controlled regime. However, when non-Darcy effects are turned off in the matrix, the recovery becomes more sensitive with increasing fracture permeability. The effect of non-Darcy flow in matrix becomes more important at higher fracture permeability.

Effect of fracture shape and size
We compare three cases with average fracture width ðb 0 Þ 0.05 m, 0.02 m (reference), and 0.09 m to evaluate the effect of fracture size. The fracture width in all the cases above was assumed to have constant width. Here, we also consider cases where the fracture shape varies. The fracture width decreases linearly with distance from the well and is defined by three parameters; the length L y , the average half-width b o and the max-to-min width ratio b max =b min : For a uniform fracture width, i.e. b max =b min ¼ 1; we get bðyÞ ¼ 2b 0 . For each of the fracture widths considered, we use two subcases with bmax b min ratio of 1 and 10 to evaluate the effect of fracture shape. The bmax bmin ratio of 10 indicates that the fracture is 10 times narrower at the tip of the fracture (y ¼ L x Þ compared to at the well (y ¼ 0Þ, as illustrated in Fig. 1.  Fig. 9 shows the simulated gas recovery for all cases with Darcy flow (a) and with non-  We see that for Darcy flow, gas recovery is very weakly sensitive to fracture shape and size, i.e. very similar profiles are seen. This indicates that the flow is mainly dominated by the matrix. However, at lower fracture width (b 0 ¼ 0:009 m), we observe that the shape of the fracture becomes more important.
For the gas transport with non-Darcy effects in matrix and fracture ( Fig. 9(b)), gas recovery becomes strongly sensitive to fracture shape. The flow of gas is delayed with decreasing fracture size and is controlled by fracture properties. This is reasonable as we have high k f compared to the matrix. Moreover, for each fracture width, we note that gas is recovered at a faster rate with high bmax bmin ratio. Berawala et al. (2019) explained that the narrower the fracture is at the end; the less space does the gas have to diffuse towards the well resulting in a local pressure buildup. In such case, more gas is produced from the matrix in the near well region compared to th regions at the far end of the fracture. For the cases with widest fractures (b 0 ¼ 0:05 m), the difference is not significant and a close to uniform production is seen along the fracture.

Interpretation using dimensionless numbers
In the above discussion, we have demonstrated how non-Darcy flow constants, matrix pore radius and fracture properties affect the production of shale gas with/without non-Darcy effects. In particular, we observed that non-Darcy effects typically shifts the flow towards

Table 4
Input parameters for simulation cases shown in Fig. 10, selected such that ω ¼ αβ is nearly constant for 4 values: ω i ¼ 10 À 3 ; 10 À 2 ; 10 À 1 and 10 0 . Other unspecified parameters are given by reference case values in Table 2. fracture-controlled regime as compared to matrix-controlled regime with Darcy flow. In this section, we demonstrate the cases where matrix properties become dominant as compared to fracture properties in the presence of non-Darcy effects. We do this by interpreting the results with the help of two dimensionless numbers α and β defined in Eq. (43). Berawala et al. (2019) showed that the production of shale gas can be classified into matrix-controlled or fracture-controlled based on the magnitude of the product of α and β (this was first demonstrated for advection-spontaneous imbibition flow in Andersen et al. (2014)). They coined the parameter ω given by: From (42), we get: If ω≪1, the flow of gas is completely controlled by the time scale of diffusion from the matrix. For larger ω, the residence time in the fracture is significant and further delays the process. Thus, when gas recovery is plotted for the cases with low ω and high non-Darcy flow constant for matrix and fracture, we expect to see matrix-controlled flow regime. To perform this test, we present eight simulation cases of gas recovery to give approximate values of the product ω ¼ αβ ¼ 1e À 3; 1e À 2; 1e À 1; 1e0. Both linear and non-linear parameters are varied as explained in Table 4. Other parameters are kept constant as listed in Table 2 unless otherwise is specified. All the simulation cases are presented in terms of gas recovery vs time in Fig. 10.
In particular, we compare two scenarios:  Table 4 are kept constant so that the recovery is only affected by the matrix non-Darcy effect.
Logarithmic time axis is used in Fig. 10 since the simulations span over a wide range of time scale.
We observe that for cases with low value of ω (1e-1), gas recovery is strongly sensitive to non-Darcy effects in the matrix. Only at high ω ¼ 1 (grey lines), we see negligible effect of non-Darcy flow in matrix. This is because at high ω, flow is dominated by fracture properties and residence time for the flow of gas in fracture plays a key role. On the other hand, for high ω, flow is controlled by the residence time in matrix and the recovery only depends on matrix properties.
We repeat the simulation cases described in Table 4 and plot them against the time scaled against τ m . However, this time for cases with no non-Darcy effect in matrix, we vary the other input parameters and adjust them in such a way that we get same ω as their corresponding pair with non-Darcy effect in both fracture and matrix. These changes are reflected in Table 5 The extended 1Dþ1D model is a useful tool to evaluate sensitivity of input parameters, to understand the role of non-Darcy effects in matrix and fracture and to qualitatively study the shale gas production system. However, the model does not consider changes in effective stresses during production. The resulting changes in fracture or matrix porosity and permeability due to geomechanical effects might alter some of the results discussed in this paper. Moreover, the model considers flow of gas only from stimulated reservoir volume, i.e. the domain affected by the hydraulic fracture. However, flow of gas from beyond the tip of fracture and cross-flow could also contribute to overall recovery. These effects should be evaluated before extending the model to field scale application.

Conclusions
In this paper, we have presented a model for production of shale gas by incorporating non-Darcy and gas slippage effects. The system consists of a single fractured extended vertically from a well perforation and is surrounded symmetrically by shale matrix. The model presented is derived in such a way that it helps us to investigate the transition from Darcy to non-Darcy flow. With the help of Forchheimer's equation, the role of non-Darcy effects were for different conditions of geometry and intrinsic properties of the fracture-matrix system. The results were further interpreted by dimensionless numbers. From the numerical investigations presented, we draw the following conclusions: � Non-Darcy flow shifts the flow regime towards fracture dominated.
The non-Darcy effects are more pronounced in fracture than matrix and cause greater increase in fracture diffusion time than matrix diffusion time. � Theoretical and numerical results indicated that the model cases could be classified according to matrix-dominated for ω≪1 (the matrix then limits the gas production) and fracture-dominated where ω � 1 (the fracture limits the gas production). � Non-Darcy flow in the matrix is significant when the flow regime is matrix dominated and non-Darcy flow in matrix significantly reduces matrix flow rate. When any of these conditions is not met, non-Darcy flow in the matrix is not relevant. This is mathematically equivalent to the dimensionless number ω≪1 and f m ref ≪1 as expressed in our model. Fig. 10. Absolute gas recovery vs scaled time. Comparative test where ω ¼ αβ is approximately for 4 values: ω i ¼ 1e À 3; 1e À 2; 1e À 1 and 1e0: Input parameters are varied in 8 tests as described in Table 4. ω seems to characterize the flow regime of the fracture-matrix system. Unspecified parameters are given by reference case values in Table 2.
� The gas recovery is more sensitive to the size and shape of the fracture in presence of non-Darcy flow compared to Darcy flow. � At sufficiently high fracture permeability (k f > 100 mD for our base case), diffusion time in the fracture reduces. Recovery profiles then become very similar and unchanging with permeability both for Darcy and non-Darcy flow. The impact of non-Darcy flow in the matrix becomes more sensitive at high fracture permeability.
� The magnitude of the f-factor helps to quantify the transition of Darcy flow to non-Darcy flow. At high non-Darcy flow constant C β , non-Darcy effects in the matrix show greater sensitivity due to lower f m ref value and should be considered for evaluation of shale gas production.
� Because of the parameters appearing in ω, we can conclude that matrix properties will control the production if f m ref is low, matrix is much less permeable than fracture, the fracture volume is low and the fracture spacing is large.

Declaration of interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  Fig. 11. Scaled gas recovery vs scaled time. Comparative test where ω ¼ αβ is constant for 4 values: ω i ¼ 1e À 3; 1e À 2; 1e À 1 and 1e0: Input parameters are varied in 8 tests as described in Table 5. ω seems to characterize the flow regime of the fracture-matrix system. Unspecified parameters are given by reference case values in Table 2.

CRediT authorship contribution statement
System 1-Fracture diffusion in y-direction.
We set ∂ X P ¼ 0 and ∂ T b A g ¼ 0 indicating no adsorption. This gives: System 2-Fracture/Matrix diffusion and desorption and flow in the x-direction.
Here, we assume no flow in y-direction and set ∂ y P ¼ 0, we get: The system is further split into two subsystems: 1) We only consider diffusion of free gas with adsorbed mass held fixed and 2) equilibrate adsorbed gas with free gas in the matrix.

Fracture matrix diffusion.
With no desorption, we set ∂ T b A g ¼ 0 and solve the diffusion system: Desorption. No flow is considered.
The conserved property here is P þ b A g . From the definition of the Langmuir isotherm and mass conservation we obtain: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where p eq is the equilibrium pressure adjusted for b a g after fracture-matrix diffusion of free gas. G denotes conserved property with units of pressure, G ¼ p g þ b a g ðp g Þ: The equilibrium pressure p eq is then scaled and returned from the adsorption-correction. Refer to Berawala et al. (2019) for detailed procedure.

B. Discretization
Let the y-axis be discretized into j ¼ 1 : N y cells and the matrix into i ¼ 1 : N x cells. System 1-Fracture diffusion The scaled (half) width BðyÞ is constant for a given cell j denoted by B j . The conserved property is P, which is integrated over the grid cell gives: The flux is further discretized as: At the fracture boundaries we set: At the fracture-matrix interface we have: