A Mixed Element Method for the Desorption-Diffusion-Seepage Model of Gas Flow in Deformable Coalbed Methane Reservoirs

We present a desorption-diffusion-seepage model for the gas flow problem in deformable coalbed methane reservoirs. Effects of fracture systems deformation on permeability have been considered in the proposed model. A mixed finite element method is introduced to solve the gas flow model, in which the coalbed gas pressure and velocity can be approximated simultaneously. Numerical experiments using the lowest order Raviart-Thomas (RT 0 ) mixed element are carried out to describe the dynamic characteristics of gas pressure, velocity, and concentration. Error estimate results indicate that approximation solutions could achieve first-order convergence rates.


Introduction
Coalbed methane (CBM) is an abundant, low cost energy source, which has become a viable alternative source to conventional fuel.CBM reservoirs are dual porosity systems consisting of coal matrix and fracture network [1].In CBM reservoirs, most of coalbed gas is stored in the coal matrix as adsorbed gas, and only a small amount is stored in fracture systems as free gas.Comparing with conventional natural gas reservoirs, the effect of fracture systems deformation on permeability is significant in CBM reservoirs [2][3][4].Meanwhile, the gas flow in fracture systems is driven not only by pressure field but also by gas concentration field [5,6], so we need to consider Darcy seepage and Fick diffusion simultaneously in multifield.As coalbed gas is released from CBM reservoirs, fracture system pressure reduces, resulting in gas diffusion in the coal matrix and gas desorption from coal matrix surface to fracture systems.Therefore, the gas flow process in CBM reservoirs includes gas desorptiondiffusion in coal matrix and gas seepage in fracture systems [7][8][9][10].Usually, the fracture systems are considered as the gas flow passage and the coal matrix is treated as the source for fracture systems.
Based on the above understandings about coalbed gas flow mechanism, a series of conventional mathematical and numerical models [5,6,11,12] have been developed, obtaining some useful computational and simulation results.Young [11] presented a computer model for CBM reservoirs to determine key data and to describe the variability in CBM reservoir properties.Unsal et al. [12] proposed a numerical model for multiphase flow in CBM reservoirs using a fracture-only model.Thararoop et al. [5,6] developed a compositional dual porosity, dual permeability CBM numerical model to describe the pressure and concentration dynamic behaviors.Nevertheless, the existing CBM numerical simulation methods are based on finite difference scheme and mainly focus on gas pressure and concentration, ignoring gas flow velocity.
Mixed finite element method was initially introduced by engineers [13] in the 1960s and has been applied to many areas such as solid and fluid mechanics, which could approximate both vector variable (e.g., the fluid velocity) and scalar variable (e.g., the pressure) simultaneously and give a high order approximation of both variables.Comparing with the standard finite element method only employing a single finite element space, mixed finite element method need construct two different finite element spaces.Raviart and Thomas [14] introduced the first family of mixed finite element spaces for second-order elliptic problems.Nedelec [15] extended these spaces to three-dimensional problems.Then, Brezzi et al. [16,17] and Chen and Douglas [18] presented many mixed finite element spaces.

Mathematical Problems in Engineering
In addition, the study about numerical method and numerical analysis for fluid flow problems in porous media has been a research hotspot for the past decades [19][20][21][22][23][24], due to its wide applications in various engineering areas.Douglas et al. [19][20][21] have done a lot of useful work for the fluid flow problem in porous media.For the incompressible miscible displacement problem in porous media, Wang [22] introduced ELLAM-MFEM to solve equation of the problem.A component-based Eulerian-Lagrangian method [23] was used to treat the multicomponent, multiphase flow problem in porous media.Li and Sun [24] presented an unconditional convergence Galerkin-mixed element approximation for the incompressible miscible flow problem.However, until now, we have not found any paper considering mixed element approximation scheme for the gas flow problem in deformable CBM reservoirs.
The aim of this paper is to study the gas flow problem in deformable CBM reservoirs using mixed element method.For this purpose, we introduce a mixed finite element method to approximate the coalbed gas pressure and velocity simultaneously.Another objective of this paper is to carry out some numerical experiments using the lowest order Raviart-Thomas mixed element, so as to show the convergence rate of the mixed element approximation method, and to analyze coalbed gas flow characteristics.
This paper is organized as follows.In Section 2, considering the effects of fracture system deformation, we present a desorption-diffusion-seepage model of gas flow in deformable CBM reservoirs for the two-dimensional problem.In Section 3 we introduce a mixed finite element approximation scheme for the coalbed gas flow model.In Section 4 numerical experiments using the lowest order Raviart-Thomas (RT 0 ) mixed element are carried out.

Coalbed Gas Flow Problem in Deformable CBM Reservoirs
In this section, considering the pressure dependence of permeability and porosity, we derive a mathematical model for the gas flow problem in deformable CBM reservoirs.

Basic Assumptions and Gas State Descriptions
(i) CBM reservoirs are treated as dual porosity systems consisting of coal matrix and fracture network.
(ii) CBM reservoirs are isothermal.When pressure variation is not significant, we can assume the gas viscosity is constant under isothermal conditions [25].
(iii) The gas diffusion process in the coal matrix is pseudostatic and described by Fick's first diffusion law.
(iv) Gas absorption/adsorption is described by Langmuir adsorption isotherm.
The equation of real gas state is where  is the gas pressure,  is the gas volume,  is the amount of substance,  is the gas temperature,  is the gas mass,  is the gas molar mass,  is the gas deviation factor, and  is the universal gas constant.From (1) the gas density  can be written as According to the Langmuir isotherm and the state equation of real gas, the adsorbed gas concentration   in coal matrix and free gas concentration   in fracture systems can be written, respectively, as where   is the gas pressure in coal matrix,   is the maximum adsorption concentration of the coal matrix,   is the Langmuir pressure,   is the gas pressure in fracture systems, and   is the porosity of fracture systems.

Gas Flow Model.
Considering the effects of fracture system deformation, fracture system permeability tensor K  (,   ) can be written as where  0 is the initial gas pressure in fracture systems, K 0 is the fracture system permeability tensor under  0 , and  is the permeability modulus with respect to gas pressure.Coalbed gas flow in fracture systems is driven by pressure field and concentration field simultaneously.According to Darcy's law and Fick's diffusion law, the motion equation of gas flow in fracture systems can be written as where u is the total volume velocity of gas flow driven by multifield,   is the gas viscosity, and D  is the gas diffusion coefficient in fracture systems.Then, according to the mass conservation law, the continuity equation of gas flow in fracture systems is where   is the gas volume of interporosity flow from coal matrix to fracture systems.Using the state equation of real gas, the divergence term of (6) could be unfolded as the following form: Similarly, the capacity term could be unfolded as the following form: where   is the compressibility factor of   and   is the compressibility factor of , defined as For the expression brevity, define and substitute motion equation ( 5) into ( 7); the continuity equation ( 6) can be rewritten as follows: Since the seepage velocity u is very low, the quadratic term u 2 usually can be omitted in engineering application.Therefore, the more common form of continuity equation is as follows: In addition, there are some other studies [26,27] considering the influence of quadratic term u 2 .For example, Ranjbar et al. [26] have considered the gas density variation in space and solved the nonlinear equations using semianalytical methods.However, they did not consider the effect of gas desorption and deformable media.But their method may form a basis for application of semianalytical methods for problem like the one studied here and it will be useful to the readers.
According to Fick's first diffusion law, the pseudosteady state equation of gas diffusion in coal matrix can be written as where   is the average gas concentration of a coal matrix block,   (  ) is the adsorbed concentration on the coal matrix surface and is defined as is the gas diffusion coefficient in coal matrix block, and  is the geometrical factor.Combining ( 11) and ( 13), with corresponding initial and boundary conditions, we have the coalbed gas flow model as follows: where Ω is the objective region for the two-dimensional problem with boundary Ω, [0, ] is a time interval with  ∈ (0, ∞], and   is the region of a coal matrix block.

Notations and Mixed Element Approximation Scheme
In this section, we will present the mixed element approximation scheme for the coalbed gas flow model and give some notations.

Some Notations and Basic Approximation Results
. For the convenience of subsequent analysis, we first give notations and basic approximation results.Define the inner product on Ω and its norm: Recall that the Sobolev space   (Ω) is the closure of  ∞ (Ω) in the norm Define the function spaces ,  and their norms as follows: Let  ℎ ×  ℎ ⊂  ×  be the mixed element function spaces such as RT  with index  and discretization parameter ℎ.Now we define the RT projection Π ℎ :  →  ℎ , satisfying Define standard  2 projection  ℎ :  →  ℎ , satisfying 3.2.Mixed Element Approximation Scheme.We consider the case that the seepage-diffusion tensor K(,   ) is a diagonal matrix and  , > 0 ( = 1, 2) for the two-dimensional problem.That is to say, there exist positive constants   and   such that, for all  ∈  2 , For the practical situation, there exist   and   such that the scalar coefficient (,   ) satisfies Also, we can easily verify that the coefficients  , (  ), (  ), and   (  ) satisfy the Lipschitz continuity.Since K(,   ) is a diagonal and invertible matrix, from (5) we have that Using the integration by parts and applying the boundary condition, we can define the mixed weak formation, which is to find   ∈ , u ∈  such that where n is the unit exterior normal vector to the Ω.
Replacing the original pressure   and velocity u by their approximations, we get the semidiscrete mixed element approximate problem, which is to find  ℎ ∈  ℎ , u ℎ ∈  ℎ such that Let Δ > 0,  = /Δ, an integer, and   = Δ,  = 0, 1, . . ., .We can define the full discrete mixed element approximate scheme with backward Euler time-discretization as follows: Since gas diffusion equation ( 13) is an ordinary differential equation, when   ℎ is known, we can get   (  ) as follows: Then, the gas volume of interporosity flow from coal matrix to fracture systems could be calculated as follows: In practical calculation, alternative and iterative method is used to solve ( 25)-( 27), and specific calculation procedure is listed as follows.
Let  ℎ be a quasiregular triangulation for the twodimensional rectangular region Ω with mesh size ℎ.For a triangular element ( 1  2  3 ), let x  = (  ,   ) be the coordinate of vertex   and   the edge opposite to vertex   .Let n  be the outward unit vector on the edge   and ℎ  the length of the perpendicular dropped from the vertex   onto the edge   ,  = 1, 2, 3. We denote by ϝ ℎ the set of all edges of the triangulation  ℎ .
The pressure function space  ℎ is the space of piecewise constant functions: and the dimension of the space  ℎ is equal to  1 -the total number of elements   in the triangulation  ℎ .Each function  ℎ ∈  ℎ can be represented as a linear combination: where   (x) is the basis function associated with the element   , satisfying The velocity function space  ℎ is chosen to be the RT 0 (Ω,  ℎ ), defined as and the dimension of the space  ℎ is equal to  2 -the total number of edges   ∈ ϝ ℎ .Each function k ℎ ∈  ℎ can be represented as a linear combination where V  = k ℎ ⋅ n  and   is the basis function associated with the edge   ∈ ϝ, satisfying To be specific, for a triangle element   , the basis functions   ( = 1, 2, 3) are determined by the following formulas: where x  = (  ,   ) are the coordinates of the vertices   ( = 1, 2, 3) and || is the area of triangle element   .Now we use the basis functions   (x) ( = 1 ⋅ ⋅ ⋅  1 ) and   (x) ( = 1 ⋅ ⋅ ⋅  2 ) to represent   ℎ and u  ℎ : Substituting ( 35) into ( 25) and setting  ℎ =   ( = 1 ⋅ ⋅ ⋅  1 ) and k ℎ =   (x) ( = 1 ⋅ ⋅ ⋅  2 ), full discrete mixed element approximation scheme (25) can be written as follows: We introduce matrices and vectors as follows: where Therefore, (36) can be written in a matrix form as follows: where B  is the transpose of B. Solving linear algebra problem (39), we can get approximation solutions P  and U  at th time step, simultaneously.

Numerical Examples
In this section, we carry out numerical experiments using the lowest order Raviart-Thomas mixed finite elements RT 0 for the gas flow problem in deformable CBM reservoirs.
Example 1.In order to verify the convergence, the test region is selected as unit square; that is,  1, 2, and 3, respectively.We can find that the approximation solutions  ℎ , u ℎ , and  ,ℎ could achieve first-order convergence rate (see Figure 1).
Example 2. For a simple application, the scale of objective region Ω is 100 m × 100 m and the mesh size ℎ = 2.5 m.There is a production well at the center of region Ω with the production rate   .Simulation parameters are listed in Table 4.
Substituting the parameter value into calculation procedure, we can get the numerical solutions, so as to describe the dynamic characteristic of gas flow in deformable CBM reservoirs.
Since average pressure could reflect the decrement of reservoir driving energy during production process, Figures 2 and 3 illustrate the gas average pressure in fracture systems varying with time  and show the effects of   and  on average pressure, respectively.In Figure 2, the model has the initial reservoir pressure of 10 MPa, and diffusion coefficient of coal matrix   is taken as 1 × 10 −7 , 5 × 10 −8 , 1 × 10 −8 , and 5 × 10 −9 m 3 /s, respectively, so as to analyze the effects of   on average pressure in fracture systems.
In Figure 2, we have found that pressure variation curves could be divided into 3 stages.Firstly, at the initial stage, since coalbed gas could not desorb and diffuse from matrix into fracture systems instantly, average pressure in fracture systems declines rapidly.At this stage, we have found when   is increasing, the initial stage became shorter.It is because that larger   could make fracture systems get pressure compensation from coal matrix more quickly.Secondly, at the stable stage, with gas desorption increasing from the coal matrix into fracture systems, fracture systems could get more pressure compensation; thus, the average pressure drop speed slows down.We have also found, with a bigger   , fracture systems could maintain a higher average pressure, while the stable stage would be shorter.Lastly, at the last stage, due to the decline of gas concentration in matrix gas desorption rate decreases, resulting in the pressure drop speed becoming rapid once again.Meanwhile, when   is larger, the last stage appears at an earlier time.The double porosity characteristics of CBM reservoirs have been just shown through the pressure variation process presented above.In Figure 3, the permeability modulus  is set as 1 × 10 −1 , 1 × 10 −3 , and 1 × 10 −5 , respectively, so as to analyze the effects of  on average pressure.We have found that pressure drop curves with different  are almost coincident Relative error at early time, which indicates that fracture deformation has no evident effects on permeability and pressure at this stage.Essentially, with the fracture pressure dropping, the effects of fracture deformation on permeability become significant, and, conversely, the permeability could also affect fracture pressure.That is to say, the effects of permeability and pore pressure are mutual.Subsequently, we have found pressure drop curves become apart in Figure 3. Also, we have found that the larger the  is, the lower the average pressure will be.From the definition of , we know that  reflects  the effects of fracture system deformation on the permeability K  .Thus, when  is larger, K  decreases more significantly with   dropping, such that CBM reservoirs can only get less pressure compensation from outer boundary, resulting in a lower average pressure.
Figure 4 shows the bottom hole pressure varying with time , and the diffusion coefficient   of coal matrix is taken as 1 × 10 −7 , 5 × 10 −8 , 1 × 10 −8 , and 5 × 10 −9 m 3 /s, respectively, so as to analyze the effects of   on bottom hole pressure.Since the bottom hole pressure is the direct reflection of production effects, the bottom hole pressure drops sharply from 10 MPa to about 4.25 MPa after starting production, and then there is a pressure recovery due to the gas desorption from coal matrix to fracture systems.Thus, this causes an early fluctuation before reaching a stable state.In the subsequent stable stage, we have found that a bigger diffusion coefficient   could make more coalbed gas desorb into fracture systems and result in a higher bottom hole pressure.Then, at the last stage, due to the gas desorption decreasing, the effects of   become weak and bottom hole pressure curves tend to be close.Put simply, the effects of   are evident at stable stage, while they are not so evident at the early and last stage.Coalbed gas flow in fracture systems is driven by pressure field and concentration field simultaneously; thus, the effects of gas diffusion in fracture systems need to be investigated.In Figure 5, we set fracture diffusion coefficient   as 0, 1×10 −2 , 1 × 10 −3 , and 1 × 10 −5 m 3 /s, respectively, so as to analyze the effects of gas diffusion in fracture systems.We have found that when   is larger, the bottom hole pressure is higher.Specifically, the pressure difference between   = 1 × 10 −2 and   = 0 is about 0.35 MPa all through the stable stage in Figure 5.That is to say, existence of gas diffusion in fracture systems could make fracture systems get more pressure compensation from outer boundary, resulting in a higher bottom hole pressure.Thus, it is very meaningful to consider gas diffusion in fracture system.
In order to illustrate the gas desorption dynamic process from coal matrix into fracture systems, we show the gas average desorption rate varying with time  in Figures 6 and 7.In Figure 6, the diffusion coefficient   is taken as 1 × 10 −7 , 5 × 10 −8 , 1 × 10 −8 , and 5 × 10 −9 m 3 /s, respectively, so as to investigate the effects of matrix diffusion coefficient on gas desorption process.According to equilibrium desorption model, we know that gas desorption rate is affected by gas average concentration   in matrix, fracture pressure   , and diffusion coefficient   .During production process, as fracture pressure   and gas average concentration   decline, there are two stages in gas desorption process.At first half stage, due to quick drawdown of fracture pressure, gas concentration   (  ) on matrix surface declines sharply, and then the gas average desorption rate increased to the maximum.Then, at the second half stage, the desorption rate declines due to the reduction of gas average concentration   .In Figure 6, we have also found when   is larger, average desorption rate will be higher and peak desorption rate appears at an earlier moment.In Figure 7, the fracture diffusion coefficient   is taken as 0, 1 × 10 −2 , 1 × 10 −3 , and 1 × 10 −5 m 3 /s, respectively.Since the effects of   on desorption rate are indirect and reflected by means of fracture pressure   , the effects of   are not evident at the starting and ending time.However, we have still found that the larger the   is, the higher the peak desorption rate will be.
In order to reflect the remaining reserves of coalbed gas in coal matrix, we show the gas average concentration in coal matrix varying with time  in Figures 8 and 9.In Figure 8,   is taken as 1 × 10 −7 , 5 × 10 −8 , 1 × 10 −8 , and 5×10 −9 m 3 /s, respectively.We have found that the gas average concentration declines slowly at the early stage on account of the existence of desorption delay phenomenon.Then, due to gas desorption increasing from coal matrix into fracture systems, the average concentration in coal matrix declines quickly at the later stage.At last stage, due to gas desorption rate decreasing, the concentration drop speed becomes slow once again.Also, in Figure 8, we have found concentration curves are almost coincident at early time.At subsequent stages, when   is larger, the average concentration is lower and the concentration drop speed is quicker.It is because that larger   could accelerate gas desorption rate.In Figure 9, we have found that the effects of initial permeability  0 are not evident at the early stage, due to the existence of desorption delay phenomenon.And then, at the later half stage, the average concentration increases with  0 increasing.It indicates that the CBM reservoir with higher permeability could get more pressure compensation from outer boundary, resulting in higher pore pressure and matrix concentration.

Conclusions
(1) In this paper, we present a desorption-diffusionseepage model of gas flow in deformable CBM reservoirs, in which the effects of fracture system deformation have been considered.
(2) We introduce a mixed finite element method to approximate the coalbed gas pressure and velocity, simultaneously, and establish error estimates for the numerical solutions.
(    could achieve first-order convergence rate.Numerical results in Example 2 have shown that the existence of fracture system deformation phenomenon will result in a lower gas average pressure.

9 Figure 2 : 5
Figure 2: Gas average pressure of fracture systems with different   .

Figure 3 :
Figure 3: Gas average pressure of fracture systems with different .

Figure 4 :Figure 5 :
Figure 4: Bottom hole pressure of CBM reservoirs with different   .

Figure 6 : 2 D
Figure 6: Gas average desorption rate with different   .

Figure 7 :
Figure 7: Gas average desorption rate with different   .

Figure 8 :
Figure 8: Gas average concentration in coal matrix with different   .

Figure 9 :
Figure 9: Gas average concentration in coal matrix with different  0 .

Table 1 :
Relative errors and convergence rates for  ℎ .

Table 2 :
Relative errors and convergence rates for u ℎ .

Table 3 :
Relative errors and convergence rates for  ,ℎ .
) We carry out numerical experiments using the lowest order Raviart-Thomas mixed finite elements RT 0 .Numerical results in Example 1 indicate that the approximation solutions  ℎ , u ℎ , and  ,ℎ