Chemical Potential-Based Modeling of Shale Gas Transport

the Shale gas plays an increasingly important role in the current energy industry. Modeling of gas ﬂ ow in shale media has become a crucial and useful tool to estimate shale gas production accurately. The second law of thermodynamics provides a theoretical criterion to justify any promising model, but it has been never fully considered in the existing models of shale gas. In this paper, a new mathematical model of gas ﬂ ow in shale formations is proposed, which uses gas density instead of pressure as the primary variable. A distinctive feature of the model is to employ chemical potential gradient rather than pressure gradient as the primary driving force. This allows to prove that the proposed model obeys an energy dissipation law, and thus, the second law of thermodynamics is satis ﬁ ed. Moreover, on the basis of energy factorization approach for the Helmholtz free energy density, an e ﬃ cient, linear, energy stable semi-implicit numerical scheme is proposed for the proposed model. Numerical experiments are also performed to validate the model and numerical method.


Introduction
Shale gas has become a significant energy resource over the last decade. Shale gas refers to natural gas composed of primarily methane, which is trapped within the pores of fine-grained sedimentary rocks with rich micropores and relatively low permeability. The shale gas reservoirs differ from conventional natural gas reservoirs that apart from free gas in the pores and fractures, a certain amount of gas is adsorbed onto the solid surfaces, and as a result, it can not only store shale gas but also generate gas [1][2][3][4]. Experimental investigations have indicated that adsorbed gas storage capacity is primarily affected by shale reservoir conditions, such as temperature, pressure, and shale matrix pore structures [2,3,5,6]. A number of adsorption models have been developed to describe the methane adsorption in shale gas reservoirs. The Langmuir model [7] and Dubinin-Radushkevich (D-R) model [8] are the most popularly employed models to describe the gas adsorption in micropore-rich materials. Methane under shale formation conditions usually stays at the supercritical state, and consequently, the classical models that require a saturation pressure cannot be applied to describe the gas adsorption under supercritical conditions [3,9]. By the use of gas density rather than gas pressure, modified formulations of the Langmuir and D-R models have been developed for supercritical conditions [3,9]. In this paper, it is shown that the use of gas density will also be beneficial for ensuring thermodynamical consistency of the models.
Modeling of gas flow in shale media plays a crucial role in predicting shale gas production [10][11][12][13][14][15][16][17][18][19]. For gas flow in tight porous media, the most remarkable phenomenon is the so-called Klinkenberg effect [20], which results from slip flow of gas molecules through very small pores. This effect leads to the apparent permeability that is generally greater than the absolute permeability of a porous medium [14,15]. By using the apparent permeability, the shale gas flow equation can be simply formulated as the form of Darcy's law, which states that gas velocity is proportional to the pressure gradient [12][13][14][15][16][17][18][19]21].
As general principles, the fundamental laws of thermodynamics play a significant role in modeling of various physical problems [22][23][24]. Specially, the second law of thermodynamics states that any spontaneous process in an isolated system will always lead to an escalation in the entropy of this system. In terms of the second law of thermodynamics, for an isothermal system, any dynamical process should obey an energy dissipation law [22], and a promising model should preserve this property. Nevertheless, energy dissipation laws for modeling of shale gas transport have been scarcely studied so far. In this paper, a new mathematical model is proposed, which uses gas density instead of pressure as the primary variable and introduces chemical potential gradient instead of pressure gradient as the primary driving force. It is rigorously proved that the proposed model obeys an energy dissipation law, and thus, it is thermodynamically consistent (that is, it obeys the second law of thermodynamics).
Numerical algorithms that preserve the energy dissipation law at the discrete level, the so-called energy stable methods, are preferred as well [25][26][27]. In general, it is a quite challenging issue to construct such methods since the commonly used explicit or implicit scheme could not inherit a discrete energy dissipation law. The energy factorization (EF) approach for the Helmholtz free energy density [26] is a newly developed approach to design efficient energy stable numerical methods. An appealing feature of this approach is that it leads to linear, easy-to-implement, and energy stable numerical schemes, and this advantage is more notable for numerical simulation of realistic fluids. Due to its excellent features, this approach has been successfully extended to phase-field models [28,29]. In this work, using the EF approach, an efficient, linear, energy stable semiimplicit numerical scheme is constructed for the proposed model of shale gas transport.
The rest of this paper is organized as follows. In Section 2, a chemical potential-based model of gas flow in shale media is proposed, which is proved to obey an energy dissipation law. An efficient, linear, energy stable semi-implicit numerical method is proposed in Section 3. In Section 4, numerical experiments are performed to validate the proposed model and numerical scheme. Finally, some concluding remarks are provided in Section 5.

Model Equations
In this section, a chemical potential-based model of shale gas transport is proposed. The chemical potential is defined as the derivative of Helmholtz free energy density. In Appendix A, we elaborate on the Helmholtz free energy density determined by Peng-Robinson equation of state [30].
Molar density (mol m -3 ) of methane is denoted by c. For specific temperature TðKÞ, the Helmholtz free energy density, denoted by f ðcÞ, is a function of molar density c as described in Appendix A. The chemical potential is defined as the derivative of Helmholtz free energy density function f ðcÞ with respect to molar density where μ is the chemical potential (Pa·mol -1 m 3 ) and f ðcÞ is the Helmholtz free energy density (Pa). For specific temperature and pressure, a cubic equation is solved to obtain molar density of methane, and moreover, the solutions of c are not unique in general. In contrast, for specific temperature and molar density, the pressure can be uniquely and explicitly calculated from the Peng-Robinson equation of state [30], which is formulated in Appendix B. Consequently, under the constant temperature, molar density is preferred to pressure as the primary variable in numerical simulation.
For specific temperature, the pressure becomes a function of molar density c only. Moreover, the pressure relates to the Helmholtz free energy density and chemical potential as [22] where p is the pressure (Pa). The chain rule gives the relationship between the pressure gradient and chemical potential gradient.
Let ϕ denote the porosity. In tight reservoirs with abundant micropores, the Klinkenberg effect cannot be neglected, and thus, the apparent permeability [14,15] is expressed as where b is the Klinkenberg slippage factor (Pa) and K 0 is the intrinsic permeability (md). The slippage factor accounts for the slippage effect on permeability of gas in reservoirs. Various formulations for the slippage factor have been developed in the literature, and the following formulation [13] is used in this work: where η is the gas viscosity (Pa·s), σ is the comprehensive coefficient, R is the universal gas constant (m 3 ·Pa·mol -1 ·K -1 ), and M w is the molar weight of gas (g·mol -1 ). On the basis of the apparent permeability, the velocity can be described by the form of Darcy's law as where u is the Darcy velocity (m·s -1 ) and K app is the apparent permeability (md). Applying relation (3), the following chemical potential-based velocity formulation is obtained: In the shale gas reservoir, there is a large amount of gas adsorbed onto micropores in addition to free gas. The adsorbed gas has no mobility unless it converts to free gas, and consequently, it is assumed to have no contribution to the free energy. In the context of shale gas, c stands for the free gas density, and meanwhile, c ads represents the adsorbed gas molar density. In order to describe the adsorption of gas, there have been several models proposed in the literature, which can be classified into two classes according to their primary variables, namely, the pressure-based models and density-based models. One famous pressure-based model is the classical Langmuir isotherm adsorption model [7]. The pressurebased models are not suitable to describe the gas adsorption under supercritical conditions [3,9]. The density-based models have been developed and extensively applied to characterize the adsorption that occurs in shale media. Two modified adsorption models using molar density as the primary variable are described in Appendix C, which are employed in this paper due to their consistency to the chemical potential-based model. The mass accumulation of both free and adsorbed gas is given by c tot = ϕc + ð1 − ϕÞc ads , where c ads is the adsorbed gas density (mol·m -3 ). The mass balance equation can be expressed as Furthermore, a density-dependent function is defined as and then, (8) can be rewritten as For an isothermal dynamical system, the second law of thermodynamics leads to a certain energy dissipation law [22]. A distinctive feature of the proposed model is that it obeys an energy dissipation law, and as a consequence, it is thermodynamically consistent. Let Ω be a connected and smooth space domain with the boundary ∂Ω, and further, let n denote the normal unit outward vector to ∂Ω. Multiplying (10) by μðcÞ and then integrating it over the domain, it follows that The definition of chemical potential and integration by parts lead to ð Substituting (7) into (12) The right-hand side term of (13) represents the energydissipation rate. In order to demonstrate this property clearly, the adsorption effect is disregarded and the no-flow boundary condition is applied on the entire boundary of the domain, i.e., u · n = 0. In this situation, total free energy within the domain is defined as and from (13), the following energy dissipation law can be derived: which implies that total free energy would be dissipated over time until an equilibrium state is reached.

Numerical Method
In this section, a linear, efficient numerical scheme is proposed, and it is able to ensure the energy dissipation law at the discrete level. Denote by τ the time step size and also denote by t n = nτ the time. In the time discrete scheme, c n stands for the approximation of molar density c at time t n .
3.1. Energy Factorization Approach. The proposed scheme employs the EF approach for the Helmholtz free energy density [26], which leads to the linear discrete chemical potential inheriting the energy dissipation property. The Helmholtz free energy density can be expressed as a sum of three contributions where Here, R is the universal gas constant and the parameters α and β are described in Appendix A. The energy factorization approach gives the discrete chemical potential as where In (22), GðcÞ is an intermediate energy function, which is defined as

Geofluids
where λ is the dimensionless stabilization constant. Apparently, the discrete chemical potential μ n+1 is a linear function of c n+1 only. From the physical point of view, the boundedness of molar density c is assumed as where ρ is the upper bound of molar density. Let ε 0 = ϱβ. If the stabilization constant λ in (24) is chosen such that It has been proved in [26] that the discrete chemical potential (20) satisfies the following energy inequality: 3.2. Discrete Scheme. On the basis of the discrete chemical potential given in (20)-(23), the proposed semi-implicit time discrete scheme reads as follows:

Geofluids
The apparent permeability is explicitly calculated as where the pressure can be directly obtained by (2) using c n .

Combining (28) and (29) yields a single linear equation of c n+1 as
which is easy to be solved due to the fact that μ n+1 is a linear function of c n+1 .       Figure 9: Initial adsorption density in example 2.

Geofluids
As depicted in Appendix C, it is assumed that φðcÞ > 0. It is ready to prove that the proposed numerical scheme preserves a discrete energy dissipation law. Multiplying Equations (28) and (29) by μ n+1 and c n ∇μ n+1 , respectively, and then integrating them over the domain Ω, we obtain It follows from (32) and (33) that          10 Geofluids On the other hand, it follows from (17) that As a consequence of (34) and (35), the following discrete energy dissipation law is obtained Equation (31) can be further spatially discretized using the cell-centered finite difference (CCFD) method, which is a popular discrete method in various fluid flow fields [25,26].

Numerical Examples
In this section, the performance of the proposed model and numerical method is demonstrated through some numerical tests. The spatial domain in all numerical tests is a square domain ½0, 10 × ½0, 10, and the length unit is meter. A uniform rectangular mesh with 3600 grid cells is used to divide the domain. The gas is methane and its physical properties are listed in Table 1. The modified Dubinin's adsorption model is used, and the initial velocity is taken to be zero for all tests.
The Hagen-Poiseuille equation is used to calculate the intrinsic permeability of gas flow in cylinders [13].
where r refers to the averaged pore diameter (nm).

Example 1: Verification of the Energy Dissipation
Property. In this example, the aim is to verify the energy dissipation property of the model and numerical scheme. For  11 Geofluids this purpose, the adsorption of methane is disregarded and the no-flow boundary condition is applied on the entire boundary of the domain, i.e., u · n = 0. In this situation, total free energy is expressed as From (36), the following energy dissipation law at the discrete level is obtained which implies that total free energy would be decreasing with time steps until an equilibrium state is reached. The physical parameters used in this example are listed in Table 2.
The stabilization constant in (24) is taken as λ = 0:0375. The time step size is chosen as τ = 20 hours, and 50 time steps are simulated. The initial methane density is generated by a random way.
The proposed model and numerical scheme are employed to simulate this problem. Dynamics of molar density and pressure is illustrated in Figures 1 and 2. For comparison, the classical pressure-based model is also applied to simulate the same problem, and the corresponding molar density and pressure are shown in Figures 3 and   4. The classical model are described in detail in Appendix D. Although two models produce slightly different results, it is clearly observed that due to spatial inhomogeneity of molar density, the fluid system tends to reach an equilibrium state where molar density becomes uniform in space. Figure 5 depicts total free energy profiles of the proposed model and numerical scheme. Total energy is always decreasing monotonously over time, and thus, theoretical result is validated. For comparison, total free energy computed by the classical model is also plotted in Figure 6. Despite the decrease of total energy occurring in the simulation of the classical model, total energy is obviously oscillating and not monotonously decreasing over time as shown in the zoom-in plot.
The performance of two models in preserving the property of total mole conservation is also compared. The relative error of total moles is defined as Figure 7 depicts the relative errors of total moles computed by two models. It can be observed that the proposed model can guarantee this key property very well within the range of roundoff errors, whereas the classical pressurebased model can cause a nonnegligible mass loss even in this simple problem. This means that the use of molar density as      13 Geofluids the primary variable has a great advantage in guaranteeing mass conservation.

Example 2.
A gas reservoir with multiple highpermeable layers is considered in this example. The intrinsic permeability of the reservoir is illustrated in Figure 8. The layers with high permeability are denoted by where L x = L y = 10 m. The initial gas density is uniform and equal to 200 mol/m 3 , while the initial adsorption density is illustrated in Figure 9. The gas is flowing out from the production boundary Γ, which locates on the right end of the domain The boundary condition on the production end is specified as where c b is the gas production density. The rest boundaries are closed, and no mass transfer across these boundaries takes place. The physical parameters used in this example are listed in Table 3. The stabilization constant in (24) is taken as λ = 0:0148, and the time step size is chosen as τ = 0:8 hour.
Tests with c b = 20,40,60 mol/m 3 are performed. The gas density distributions at various times computed with c b = 20,40,60 mol/m 3 are illustrated in Figures 10-12, while the corresponding pressure contours are shown in Figures 13-15. The cumulative gas production profiles under different values of c b are also plotted in Figure 16. These results indi-cate that the production density has a great effect on the gas production; in fact, the increase of c b can largely reduce the cumulative gas production. The adsorption densities at different times computed with c b = 20,40,60 mol/m 3 are shown in Figures 17-19. The total adsorption amounts under different values of the production density are illustrated in Figure 20. It is clearly observed that the remaining total adsorption amounts are against the increase of production density. This means that a large production density may not be preferred for the sake of achieving optimal gas production.

Conclusions
A new mathematical model of gas flow in shale media has been proposed. Different from the existing models, the proposed model uses gas density instead of pressure as the primary variable, and it employs chemical potential gradient rather than pressure gradient as the primary driving force. This distinctive feature brings up with thermodynamical consistency of the proposed model; that is, the model obeys an energy dissipation law, which implies the satisfaction of the second law of thermodynamics. Physically, energy stable numerical methods that preserve a discrete energy dissipation law are strongly demanded for numerical approximation of the model. For this purpose, on the basis of the energy factorization approach, an efficient, linear, energy stable semi-implicit numerical scheme is proposed for the proposed model, which inherits an energy dissipation law at the discrete level. Numerical experiments are performed, and numerical results show that the model and numerical method can not only ensure energy stability but also provide physically reasonable results.