On a boundary layer problem related to the gas flow in shales

The development of gas deposits in shales has become a significant energy resource. Despite the already active exploitation of such deposits, a mathematical model for gas flow in shales does not exist. Such a model is crucial for optimizing the technology of gas recovery. In the present article, a boundary layer problem is formulated and investigated with respect to gas recovery from porous low-permeability inclusions in shales, which are the basic source of gas. Milton Van Dyke was a great master in the field of boundary layer problems. Dedicating this work to his memory, we want to express our belief that Van Dyke’s profound ideas and fundamental book Perturbation Methods in Fluid Mechanics (Parabolic Press, 1975) will live on—also in fields very far from the subjects for which they were originally invented.


Introduction
A wealth of information about the geological structure of shales covering length scales from the pore size to deposit size can be found in the article of Silin and Kneafsey [1] (see also [2]).The structure of a typical microelement can be represented as follows (Fig. 1): there exists a matrix, which is a porous and perhaps even fissurized medium.Its permeability, although sometimes small, nevertheless is sufficient for consideration by the classical approach of subterranean fluid mechanics [3][4][5].Within the matrix, inclusions are distributed, having moderate, and sometimes even large, porosity, but very low permeability, on the order of a nanodarcy (∼10 −21 m 2 ) or less; the standard permeability of oil and gas deposits is 100 millidarcys (∼10 −13 m 2 ).The low permeability of the inclusions is x Fig. 1 A low-permeability inclusion (gray) is embedded within a porous matrix shown by the white, dotted region, which may be fissurized.When exploitation begins and the pressure in the matrix is lowered, a boundary layer of gas flow (dark gray) develops and grows in width over time.The flow in the boundary layer is determined as a function of time and the coordinate x, which is normal to the inclusion boundary and directed inside the inclusion related to the smallness of pores or channels whose diameter, in contrast to ordinary rocks forming oil and gas deposits (1-100 µm), is on the order of nanometers.
In typical shales, it is thought that a large majority of the gas is contained within the inclusions.Due to the extremely low permeability of the inclusions, classical subterranean fluid mechanics would predict a negligible amount of flow on a realistic time scale.However, experimental data point to substantial amounts of flow for reasons that are not fully understood.Here we develop a mathematical model of gas flow in shales to explain this observation, based upon the hypothesis that the inclusions undergo a nanostructural transformation when subjected to a very large pressure gradient, causing a substantial increase in permeability.The mathematical model makes predictions about the expected gas flow rate from a shale deposit.

Mechanism of gas flow from the inclusions
When the exploitation of a gas deposit in shales begins, the gas pressure in the matrix is reduced relatively rapidly [3,5].Due to the extremely low permeability of the inclusions, there appears rather quickly a pressure difference at the boundaries between the inclusions and the matrix.The permeability of the inclusions is small, but nevertheless nonzero, and therefore in the narrow region at the boundaries between the inclusions and the matrix there appears a very large pressure gradient.
Our basic hypothesis is formulated as follows.Due to the large pressure gradient, there appears a boundary layer in the inclusions (shown in dark gray in Fig. 1) where the shale undergoes a nanostructural transformation, corresponding to breakage and/or deformation of the pores and channels in the material.This causes the boundary layer to become significantly more permeable for gas.Since the pores are narrow, the influence of the Knudsen number Kn, defined as the ratio of the molecular mean free path to a representative length scale, should also be assumed.This leads to an expression for permeability of the form k = F(∇ p, Kn). ( Two items should be mentioned.Firstly, due to the small thickness of the boundary layer, the tangential components of the pressure gradient are negligibly small, so the gas flow is one-dimensional, and in fact the permeability in the boundary layer depends on the component ∂ x p, where x (Fig. 1) is the coordinate along the normal to the inclusion boundary and directed inside the inclusion.Secondly, in a more general case, the permeability could also depend on the pressure p, although this is not considered in detail here (dependence on pressure, but not the dependence on the pressure gradient, was used by Silin et al. [2]).
It is important that the pores in the boundary layer have no characteristic, selected length scale, and therefore in the intermediate asymptotic interval there is no characteristic quantity having the dimension of the pressure gradient.Here is the characteristic length size of the inclusion and λ is the fundamental nanometer length scale, 1 Å, a pixel in the description of nanoscale phenomena [6].From this assumption it follows that the permeability in the boundary layer is the power function of the pressure gradient where, in a more general case, A and m can depend on the Knudsen number and A can also depend on the pressure.Indeed, if ∂ x p is denoted by z, then the ratio of permeabilities at two different values of z, z = z 1 and z = z 2 , is F(z 1 )/F(z 2 ).If we pass to different units of pressure or length, then this ratio will be equal to F(αz 1 )/F(αz 2 ), where α is a constant factor.If there exists no other characteristic argument of the dimension of pressure gradient, then the ratio of permeabilities preserves the same value, so that We have the right to select α = z −1 2 such that this equation assumes the form Differentiating by z 1 and setting z 1 = z 2 = z yields whereupon integration gives F(z) = Az m .

Model of gas flow from the inclusions
Using the basic equation of mass conservation and the Darcy law where φ is the porosity and μ is the viscosity of gas, we obtain the basic equation This equation would remain valid in the more general case where φ, μ, and A could have a dependence on p, although in the present article we do not consider this and assume that they are constants.We assume further that the gas is thermodynamically perfect and that the flow is isothermal, so that ρ( p) = C p, where C is a temperature-dependent constant.Using these assumptions, we reduce Eq. ( 8) to the final equation where is a constant that can depend on the Knudsen number.In the classical case where m = 0 and A = k, the Boussinesq-Leibenson equation is obtained, where the constant this time is equal to κ = k/(2φμ).

Boundary conditions and self-similar solution to the boundary layer problem
Equation (8) for the problem is complemented by the two boundary conditions and the initial condition where t 0 is the time of the beginning of the exploitation.P is the gas pressure in the pristine state, and p 0 is the pressure in the matrix, established rapidly-we can assume instantaneously-after the beginning of the exploitation, and after that remaining constant.The derivation of the second boundary condition in Eq. ( 12) follows the usual steps of classical boundary layer theory in a viscous fluid.It is important to note that the relative depression (P − p 0 )/P is small-this is a well-known fact in the practice of the development of gas deposits.We can now solve the mixed Cauchy-boundary value problem given by Eqs.(9), (12), and (13).The solution depends on the quantities The dimensions of these quantities in the LT [ p] class are Standard application of dimensional analysis gives so that the solution is self-similar.We denote the independent variable ξ as so that Substituting Eqs. ( 17) and (18) into Eq.( 9) with conditions given by Eqs. ( 13) and ( 12), we obtain Equation ( 19) can be solved numerically using the Runge-Kutta method.To do this, it is convenient to rewrite the equation as a first-order system The shooting method is used, whereby the system given by Eqs. ( 21) and ( 22) is treated as an initial-value problem.An integration step of 4 × 10 −6 is used, which is small enough to ensure numerical errors will be negligible.We now come to the results of numerical integration.Figure 2a, b shows graphs of f (ξ ) and d f /dξ for the case where σ = 0.95 for various values of m.The most significant feature suggested by the computations is that for m > 0 the solutions have a compact support: pressure values are different from 1 only over a finite interval 0 ≤ ξ < ξ 0 (m). Figure 2c, d shows the same quantities for the case where σ = 0.975, indicating that for a given value of m the corresponding value of ξ 0 (m) is reduced.In Fig. 3 the values of d f /dξ at ξ = 0 that are determined using the shooting method, and the values of ξ 0 (m), are plotted as a function of m for the two values of σ that are considered.
Figure 4a shows plots of the pressure gradient for σ = 0.95 when rescaled by the values from Fig. 3.In Fig. 4b, plots of rescaled permeability are given, by taking a rescaled pressure gradient to the mth power.It is instructive that the result is universal (m-independent); the universality was found by numerical integration.For σ close to one, the universality can be proved analytically.Indeed, in such cases u = d f /dξ is small, f is close to one, and u.Therefore, the third term in Eq. ( 22) can be neglected, and the equation takes the form Integration of this equation under the condition u = 0 at ξ = ξ 0 gives At ξ = ξ 0 , u = d f /dξ = 0, so that Figure 5 demonstrates the smallness of the difference between the rescaled permeability obtained by numerical integration and the parabola given by Eq. ( 26).The prediction of the gas flow rate for a given deposit is of crucial importance.According to the previous analysis, we obtain for the inflow rate per unit volume Q at time t Here S is the specific surface of the inclusions (surface per unit volume).Equation ( 27) clearly indicates that the inflow rate decreases over time according to a power law.The exponent in the decay rate varies from − 1 2 for m = 0 down to −1 as m → ∞.Classical subterranean fluid mechanics would predict a decay of the gas inflow rate according to (t − t 0 ) −1/2 , and thus the model we present predicts that higher rates of decay will be observed.

Conclusion
Gas recovery from shales is considered today to be one of the most promising directions in natural gas production.In the present article a mathematical model of the gas flow in shales is presented based on the assumption of a nanostructural transformation that significantly increases the permeability of low-permeability inclusions containing the (basic) amount of gas in shales.A boundary layer on the surface of the inclusions is formed during the process of exploitation.It is demonstrated that the flow in this boundary layer is self-similar, and a complete solution is constructed.The ultimate formula for the gas inflow rate from the inclusions to the porous matrix is obtained, and it predicts a range of decay rates over time that are higher than for classical subterranean fluid mechanics.

Fig. 2
Fig. 2 Plots of a f and b d f /dξ for the case of σ = 0.95, and plots of c f and d d f /dξ for the case of σ = 0.975

Fig. 3 Fig. 4
Fig. 3 Plots of values of a d f /dξ that are found using the shooting method as a function of m, and b ξ 0 as a function of m for the two values of σ that have been considered (a) (b)

Fig. 5
Fig.5 Plots of difference between rescaled permeability and the approximate parabolic solution, for several values of m, for the case where σ = 0.95