A coupled compressible flow and geomechanics model for dynamic fracture aperture during carbon sequestration in coal

This paper presents the development of a discrete fracture model of fully coupled compressible fluid flow, adsorption and geomechanics to investigate the dynamic behaviour of fractures in coal. The model is applied in the study of geological carbon dioxide sequestration and differs from the dual porosity model developed in our previous work, with fractures now represented explicitly using lower‐dimensional interface elements. The model consists of the fracture‐matrix fluid transport model, the matrix deformation model and the stress‐strain model for fracture deformation. A sequential implicit numerical method based on Galerkin finite element is employed to numerically solve the coupled governing equations, and verification is completed using published solutions as benchmarks. To explore the dynamic behaviour of fractures for understanding the process of carbon sequestration in coal, the model is used to investigate the effects of gas injection pressure and composition, adsorption and matrix permeability on the dynamic behaviour of fractures. The numerical results indicate that injecting nonadsorbing gas causes a monotonic increase in fracture aperture; however, the evolution of fracture aperture due to gas adsorption is complex due to the swelling‐induced transition from local swelling to macro swelling. The change of fracture aperture is mainly controlled by the normal stress acting on the fracture surface. The fracture aperture initially increases for smaller matrix permeability and then declines after reaching a maximum value. When the local swelling becomes global, fracture aperture starts to rebound. However, when the matrix permeability is larger, the fracture aperture decreases before recovering to a higher value and remaining constant. Gas mixtures containing more carbon dioxide lead to larger closure of fracture aperture compared with those containing more nitrogen.


| INTRODUCTION
Fluid flow and storage in fractured geologic formations are complex, and a comprehensive understanding of the coupled phenomena involved is of great interest in a wide range of geotechnical engineering applications, such as deep geothermal energy exploitation and carbon dioxide (CO 2 ) sequestration. [1][2][3] However, many challenges remain unsolved. For example, as a promising technique for reducing CO 2 emissions to the atmosphere, sequestration of CO 2 in deep coal seams presents a problem of loss of injectivity, especially at the early stage of injection. 4,5 In fractured reservoirs, the fractures provide the most likely conduits for fluid flow; hence, the hydraulic behaviour in fractured rock depends largely on the structure of the fracture network. On the other hand, the fluid flow is also able to induce the dynamic redistribution of mechanical stress in the reservoir, leading to changes in the hydraulic properties of fractures. It can therefore be said that the fluid flow and geomechanical behaviour of rock formations play an important role in many fields of geoscience.
The characteristics of flow and transport processes in fractured porous media have been studied extensively. Three typical modelling approaches have been proposed: the continuum model, the dual continuum model and the discrete fracture model. 6 Because of its realistic description, the discrete fracture model has received considerable attention over the last few years. [7][8][9] In a discrete fracture model, the fractures can be represented within the domain in a spatially explicitly manner, and the effects of individual fractures on fluid flow and transport can be accounted for explicitly. 10 Models have been introduced for single-phase flow 7,11 and two-phase flow in a fractured media. 8,12 However, most of the works mentioned above do not take the geomechanical effects into account and focus on the hydraulic problem or the energy production rate.
Generally, geologic materials contain pores and other cavities filled with fluids under saturated or unsaturated conditions. The pore pressure of fluid existing in the void space of a geological formation and deformation is closely related. The coupled flow-deformation problem has been studied in depth in recent decades based on the classical poroelastic theory of Biot. 13 The effective stress principle is one of the most fundamental contributions in this theory, which allows the application of the quasi-static force balance law of continuum solid mechanics to coupled fluid flow and geomechanics problems. Many historical studies existing in the literature analyse the stress-sensitivity in conventional porous media by means of quasi-static poroelasticity. 14,15 A number of works also use quasi-static poroelasticity to investigate the coupled flow-deformation problem in fractured porous media with a dual porosity model, in which the fracture deformation is described implicitly with a dynamic permeability defined for the fracture continuum. 16,17 Coupled slightly compressible fluid flow and mechanical behaviour in fractured porous media have been studied extensively. To represent the fracture deformation explicitly, the discrete fracture model has been more widely used recently in coupled fluid flow and geomechanics problems. A fracture is defined as two surfaces in contact in the discrete fracture model presented by Garipov et al,18 in which a mechanical model for the fractures is derived to describe the changes in the stress and the displacement fields through the surfaces representing the fractures. Norbeck et al 19 present the framework for a numerical model that is capable of calculating the coupled interaction of mass transfer between fractures and the surrounding rock matrix, fracture deformation and fracture propagation. Instead of solving the coupled system of equations, some works use in situ stress conditions to estimate the stress field within the reservoir, and utilize an experimental relation to compute the fracture aperture deformation due to the stress field. [20][21][22][23] However, the study of compressible fluid transport and fracture deformation induced by extra physical or chemical processes, other than the change of mechanical stress, is limited. For example, CO 2 can adsorb onto the internal walls of coal causing it to swell. Most previous works treat fracture properties as material parameters in the simulation models of naturally fractured coal seams. 24 Huang and Ghassemi 20 present a poroelastic model to investigate coupled desorption/diffusion-rock deformation phenomena including fracture deformation (both normal and shear deformation with dilation) and its impact on gas flow in naturally fractured gas shale reservoirs. Simulation results show that the time evolution of fracture aperture is strongly affected by gas desorption during production, especially in the nearwellbore region. In general, coal seams contain a large number of fractures, and these fractures provide the major conduits for fluid flow. Therefore, an explicit and realistic description of the dynamic behaviour of preexisting fractures, such as fracture opening or closure, is important and necessary to understand the flow characteristics during CO 2 injection.
To understand the dynamic behaviour of fractures induced by variation of pressure, in situ stress and adsorptioninduced coal swelling during CO 2 injection, and to overcome an inherent limitation of dual porosity model presented in our previous work, 25 this paper presents a discrete fracture model of fully coupled compressible fluid flow, adsorption and coal matrix/fracture deformation. The effect of fractures is represented explicitly with lower-dimensional interface element. The fluid flow in two interactive domains (matrix and fractures) is solved assuming flow continuity between the two domains and using the principle of superposition. The dynamic correlation between the properties of the equivalent porous medium and fracture aperture is established using matrix-fracture stress-strain models that relate deformation of fracture to effective normal stress and shear stress changes. The numerical model is implemented in the coupled thermal, hydraulic, chemical and mechanical (THCM) model, COMPASS, developed at the Geoenvironmental Research Centre by  The accuracy of the new developments to the numerical model is tested against results obtained from the literature. The developed model is then applied to investigate the temporal evolution of fracture aperture during CO 2 injection into coal seams and to quantify the influence of gas injection pressure and composition, adsorption and matrix permeability on the dynamic behaviour of the fracture aperture.

| GOVERNING EQUATIONS FOR DISCRETE FRACTURE MODEL
Compared with conventional porous media, models of fractured porous media must consider the influence of discontinuities on the flow and deformation behaviour. The discrete fracture model is employed in this work, in which the computational domain is segregated into two regions: discrete fractures as hydraulic conduits and the rock matrix as the bulk porous medium. Discrete fractures can be idealized as lower-dimensional geometric objects, for example, a 1D entity is employed to represent fractures in 2D problems, as shown in Figure 1. The whole domain Ω is in which Ω represents the entire domain, Ω m represents the domain occupied by matrix, Ω f,d is the domain occupied by the d-th fracture fracture and w d is the aperture of the d-th fracture. In this section, the theory and governing equations for fluid flow and mechanical deformation are presented for each subdomain. By modelling the fracture system and porous matrix as distinct regions within the computational domain, each has its own flow variable, namely, the gas concentrations in the fractures (C f ) and porous matrix (C m ). The subscript f denotes the matrix, and the subscript f denotes the fractures. The displacements vector (u) is defined as the primary variable for mechanical behaviour. The model assumptions are as follows: (1) the gas is free phase in the natural fractures; (2) the gas may be free or adsorbed phase in the rock matrix; (3) the rock matrix is homogenous and isotropic and the deformation is linearly elastic; (4) no phase change occurs.

| Matrix flow model
The rock matrix without fractures is represented as a continuum and flow through the matrix is governed by mass conversation. Under isothermal conditions, the general form of the Eulerian continuity equation can be expressed as 28 where n m is the matrix porosity, c i m is the concentration of the i-th gas component, c i s is the adsorption mass term for the i-th component per unit volume of medium, v m is the flow velocity of fluid, D eff is the effective diffusion coefficient, F I G U R E 1 Schematic representation of the discrete fracture model q i sm are the sink terms, q i fm represent the leakoff flow between fracture and matrix and n c is the number of gas components.
By ignoring the effect of gravity, the velocity v m can be given using Darcy's law as where k m is the intrinsic permeability of matrix and μ g is the gas mixture viscosity. u m is the gas pressure, which can be expressed in terms of the sum of the concentrations of the gas components in the gas phase using the real gas law, given by where R is the universal gas constant, T is the temperature and Z m is the compressibility factor of fluid in matrix, which in this work is calculated using the Peng-Robinson equation of state. 29 Multicomponent gas adsorption can generally be described by the extended Langmuir. Taking this approach, the amount of each gas per unit volume of medium is described as where ρ c is coal density, V i is the adsorbed concentration, given as where V i L is the Langmuir volume constant of ith gas component and b i L is the Langmuir pressure constant of ith gas component.
The effective diffusion coefficient, D i eff , in Equation 1 considers the effect of the porous medium on fluid diffusion, which can be derived from the free fluid coefficient D i m as 28 : where τ is the tortuosity factor, related to the porosity, which is calculated using widely used relation provided by Millington and Quirk. 30 Substitution of Equations 3-7 into the mass balance Equation 2 yields: where δ ij is Kronecker's delta tensor with δ ij = 1 for i = j, δ ij = 0 for i 6 ¼ j.

| Fracture flow model
During gas injection into coal seams, coupled gas flow and deformation lead to changes in fracture aperture and thus fracture permeability, with the stress-dependence of fracture properties influencing the fluid flow. Thus, an independent gas mixture flow model is required for the fractures, the mass balance equation for which can be written as 31 where n f is porosity of the fracture system, identified as unity in this work, w is the fracture aperture, c i f is the concentration of the i-th gas component in the fractures and l is the dimension along the fractures. v f is the average gas velocity through the fracture, q sf is the sink term and q i mf is the leakage flux from the sides of the fracture surface to the surrounding porous media.
Expanding the term on the left hand side of Equation 9 yields Due to the sequential implicit numerical method used in this work (details given in Section 3), the fracture aperture is updated in each time step, and the temporal derivative of fracture aperture in the second term on the left hand of equation 10 has not be considered; thus, the mass balance equation for flow in fractures is reduced to It is assumed that the fluid flow along the fracture obeys Poiseuile's law, so that the flow rate obeys the well-known cubic law and, hence, v f is given as 32 where u f is the gas pressure, expressed as where Z f is the compressibility factor of fluid in fractures.

| Mechanical behaviour in the porous matrix
The mechanical response of a poroelastic medium with fluid flow is described by Biot's theory. Under the assumption of small deformation, the stress equilibrium equation is where σ ij is the total stress tensor and F i is the component of the body force vector. Based on the principle of effective stress, the relationship between the total stress and the effective stress can be written as where σ 0 ij is the effective stress tensor, α = 1 − K/K s is the Biot coefficient, K = E/3(1 − 2v) and K s are the bulk moduli of matrix and solid skeleton respectively, E is Young's modulus of coal and v is Poisson's ratio of coal. In keeping with the usual convention in geotechnics, the stresses (total and effective) are taken as tensile positive, whereas fluid pressure is taken as positive in compression.
The deformation of the porous matrix is assumed here to be segregated into two components: (1) due to mechanical stress and (2) due to sorption in the presence of adsorptive fluids. Therefore, strain can be expressed as The stress-strain constitutive relation is defined as where G = E/2(1+v) is shear modulus, λ = Ev/(1+v)(1 − 2v) is the Lamé constant and ε e v = ε e ii is the elastic volumetric strain of the bulk porous medium.
The strain-displacement relation is expressed as where ε ij is the strain tensor and u i is the solid displacement vector. ε ad is the sorption-induced volumetric strain, which is generally described by Langmuir-type equation as where ε i L is the Langmuir adsorption strain constant of i-th component. From Equations 16 and 17, the incremental volumetric strain can be expressed as where σ = σ ii =3 is the mean normal stress or confining pressure. Using the definition of porosity, the porosity change of a deforming rock matrix can be expressed as 33 where K p is the bulk modulus of pores. The volumetric variation of the porous medium satisfies the Betti-Maxwell reciprocal theorem: Transformation of Equation 20 yields Substituting Equations 22 and 23 into Equation 21 produces the change of matrix porosity: Hence. the porosity can be expressed as where Δε eff = Δε v + Δu m K s − Δε ad is the total effective volumetric strain. Because the total effective volumetric strain Δε eff is very small, generally Δε eff ( 1 and α n m0 ) 1, Equation 25 can be simplified to give 34 where n m0 is the initial porosity. The temporal evolution of porosity is then given by The permeability varies with porosity, which can be described by the widely used cubic relationship between permeability and porosity, 33 given as where k m0 is the initial permeability of the matrix system.

| Mechanical behaviour in discrete fractures
The fracture system properties depend largely on the reservoir stress field. Two distinct mechanisms are generally applied to describe the variation of fracture aperture, the first being the variation of normal effective stress acting on the surface of the fracture and resulting in relative displacements between the surfaces of the fracture, directly altering the fracture aperture as shown in Figure 2A. The second is the shear stress acting along the fracture, which is likely to make the fracture walls slide relative to each other and may cause a change in the fracture aperture through normal dilation, as shown in Figure 2B. 20,22 Fracture initiation as well as the propagation and mechanisms related to rock failure are not addressed in this work. The fracture opening or closure in the normal direction is controlled by the normal effective stress and normal fracture stiffness. A hyperbolic model developed by Bandis et al 35 and Barton et al 36 is widely used to describe the relation between normal effective stress and the response of fracture aperture in normal closure. Figure 3 shows the relationship between normal stress and fracture aperture, formulated as where w nmax is the maximum closure of the fracture aperture, K n0 is the initial normal fracture stiffness and σ 0 n is the effective normal stress acting on the fracture. Normal fracture stiffness is a measure of the fracture's sensitivity to normal stress and can be calculated as 35 The effective normal stress acting at the fracture surface is defined as 31 The shear dilation induced normal displacement of fractures has been observed in experiments 36 ; however, shear dilation does not always occur when shear stress is applied to the fracture. Before reaching the shear strength, τ c , only tangential displacement caused by shear stress occurs and no normal displacement is observed. Shear dilation starts to trigger normal dilation when the shear stress reaches the shear strength, τ c , as shown in Figure 4. When the shear stress reaches the peak shear stress, τ p , the plastic deformation of fractured rock is generated. The rock failure is not considered in this study, and the details about shear dilation after failure can be found in literature. 22,23 The shear dilation, d s , induced by an associated shear displacement, u s,a , is given as 23 where φ is the dilation angle.  Figure 4, the tangential displacement that generates shear dilation is a function of shear stress and shear stiffness of the fractures, as follows: where K t is the shear stiffness of fractures. The final aperture of the fractures is then calculated by The fracture is usually represented using a local coordinate system, as illustrated in Figure 5. In order to update the fracture aperture, the stress field calculated in the global coordinate system should be converted into the local coordinate system along the fracture, as F I G U R E 5 Transformation between the x 0 -y 0 local coordinate system (dash line) and x-y global coordinate system (solid line) σ n = σ xx sin 2 θ − 2σ xy sinθcosθ + σ yy cos 2 θ, where θ is the orientation of the fracture.

| Coupled field equations
Substituting the partial derivatives of matrix porosity n m with respect to time from Equation 27 into the mass balance Equation 8 and substituting Equations 12-13 into Equation 9 yields the final gas flow equations: where P represents the strain matrix and the mapping vector m T = (1,1,1,0). For convenience, the governing equation for mechanical behaviour is rewritten in incremental form: Therefore, Equations 33-(35) define a model for fully coupled adsorptive gas flow and mechanical deformation. It can be seen that the mechanical and mass transfer coupling terms naturally exist in the equilibrium and flow continuity equations.

| NUMERICAL APPROACH AND SOLUTION PROCEDURE
In this work, the Galerkin variation method with finite element discretization is used to solve the above model of all partial differential equations with certain initial and boundary conditions. The unknown variables, in this case, the matrix gas concentration, C m , fracture gas concentration, C f , and displacements, u, are approximated by interpolating the variables at nodes in each element through the following functions: where N I m and N I f are the standard finite element shape functions of node I m and I f for the discretized porous domain and fracture domain, respectively. N m and N f are the set of all nodes in the discretized porous domain and fracture domain, respectively, and the symbol^denotes the approximate value of the primary variable.
To reduce the dimension of fractures, we assume the fluid concentration should be continuous at the matrixfracture interface, that is,ĉ j I m =ĉ j I f = c j , and the discrete fracture elements must be located on the edges of porous matrix elements, sharing the same nodes as shown in Figure 6. The coupling between the two flow systems is achieved by using the principle of superposition, which has been applied in many works. 7,9,37 By using the principle of superposition, the mass exchange term between the two systems can balance off, and no explicit calculation for mass exchange is required. Before superimposition of two discretized flow equations, the local coordinate system should be transformed into the global coordinate system. The transformation relation between the two sets of coordinate systems is achieved by rotation matrix: where θ is the angle between the positive x-axis in the global coordinate and the positive x 0 -axis local coordinate, as shown in Figure 5.
If the unknown variables in the above field equations (34)(35)(36) are replaced with the approximate values of Equation 37, the superimposed discretized flow and deformation equations can be obtained after multiplying the test shape function and integration by parts. The global system equation is given as where C m c j c j = Σ where t is the external traction forces on the boundary Γ t , and q i m and q fj are the fluid flux normal to the boundary surface of the porous matrix and at the fracture inflow or outflow points, respectively. In this discrete fracture model, the assembly process of the matrix and fracture equations into the stiffness matrix is shown in Figure 7.
The aforementioned numerical formulation has been incorporated into the numerical code, COMPASS, which is a finite element-based simulator of coupled thermo-hydro-chemo-mechanical behaviour in porous media. For temporal discretization, an implicit midinterval backward difference time-stepping algorithm is employed. This has been found as a suitable solution for the highly non-linear class of equations such as the current application problem. 26,38 Figure 8 presents the numerical implementation procedure for this coupling model of fluid flow and deformation. A sequential implicit numerical approach is used to couple the fluid flow and geomechanics problem in fractured media. After specifying the initial information, such as initial and boundary conditions for the primary variables, material parameters and fracture-related information (fracture location, initial aperture and orientation), both systems of equations are simultaneously solved by iteration until convergence is achieved. Time-dependent matrix and fracture properties are updated in each time step based on the numerical results of the previous time step. The process continues until the specified simulation time is reached.

| VERIFICATION
To present the fluid-solid coupling process and assess the performance of the numerical model, a verification test is presented comparing the numerical results with those published in the literature. 22,23 A horizontal well-pair pattern with one producer and one injector is simulated. Four separate fractures are located between the two wells, as follows: (a) fracture length is 200 m, interconnecting the two well bores; (b) fracture length is 100 m, connected with the injection wellbore; (c) fracture length is 100 m, connected with the production wellbore; and (d) fracture length is 100 m without a connection to either wellbore, as shown in Figure 9. The displacements at the left and bottom sides are constrained in the both horizontal directions, respectively. Isotropic in situ stresses of 20 MPa are applied on the top and the right sides, the initial pressure is 10 MPa, and the injection and production wells operate at constant pressure of 20 and 10 MPa, respectively. It is assumed that the four fractures have the same properties initially. The no-load aperture of fractures is set to 0.25 mm, and the maximum normal fracture closure is 0.24 mm. A single component fluid with a compressibility of 5.39e-9 Pa −1 is selected. The input data used for the verification are listed in Table 1, which are chosen from work of Moradi et al. 22 The numerical results at steady state flow are compared with those obtained from Gu et al 23 and Moradi et al. 22 A comparison of fracture apertures at the steady state condition obtained in the present work and in the works by Gu et al 23 Figure 10. The variation in each fracture aperture along the width of the domain is almost linear due to the fact that the pressure change between the injection and production wells is approximately linear. This test verifies the interaction between the pressure field and the fracture deformation field. It can be seen that there is a good agreement between the two results, demonstrating that the theoretical model of fracture deformation developed in this work has been accurately implemented in the numerical model.

| Simulation conditions
Coal is a fractured rock with a clearly defined natural cleat system. The interconnected fracture network provides the major conduits for the flow of pore fluids because of its higher permeability; therefore, the fracture opening and closure has a significant effect on the fluid flow behaviour. However, the cleats run perpendicular to the bedding plane, and the permeability of coal parallel to the bedding plane is significantly larger than in the direction perpendicular to the bedding plane. 39 To investigate fracture opening and closure in the direction parallel to bedding plane, simulation is restricted to two dimensional problems, and plain strain condition is considered here. In this section, a set of simulation examples of the developed model aim to illustrate the effects of coupled gas flow, adsorption and deformation in fractured coal. Figure 11 presents a schematic of the model geometry used for the simulations. The reservoir dimension is considered to be 20 m× 20 m, with the injection well located at the centre with a diameter of 0.1 m, as shown in Figure 11A. Using symmetry to reduce the computational expense and complexity, only a quarter of the reservoir is simulated, as shown in Figure 11B. The fractures near the wellbore have a direct and significant effect on the well injectivity; four fractures near the wellbore are taken into account, as shown in Figure 11C, two of which run in x F I G U R E 9 Domain of the 2D problem used for model direction with the other two being y direction. It is assumed that the coal seam is located at an average depth of 1250 m and the average overburden rock density is 2300 kg/m 3 . The initial isotropic in situ horizontal stresses can be calculated as F y = F x = − 12MPa, which are assigned to the outer boundaries. The displacement constraints are assigned normal to the inner symmetric boundaries. The initial reservoir pressure is approximately 0.1 MPa. In the simulation, this pressure is assumed to be induced by a small quantity of the gas that will be injected, that is, CO 2 , N 2 or CO 2 /N 2 mixture. For gas flow, zero-flux boundaries are applied at the outer boundaries. Two measuring points are set to monitor the variations of fracture aperture, that is, point A (0.40, 0.25) and point B (0.65, 0.50). The parameters for the simulations are listed in Table 2. Many of these parameters were chosen from the literature. [40][41][42] A constant loading speed is assumed, and the characteristic time for reaching the specified injection pressure is set 0.5 h. A series of simulations under different injection conditions are performed, as listed in Table 3, to investigate the dynamic behaviour of fractures in coalbeds. The numerical simulation results are presented in terms of (1) the impact of injection pressure and adsorption-induced coal swelling, (2) the impact of different CO 2 /N 2 mixture injection and (3) the impact of matrix permeability.

| Impact of injection pressure and adsorption induced swelling
To investigate the impact of injection pressure and adsorption-induced swelling on the evolution of the fracture aperture, four simulation tests are conducted, at different injection pressures of 4, 6 and 8 MPa, the pure CO 2 is selected as injected gas in these four simulations. In order to understand the effects of adsorption-induced swelling, a simulation scenario without adsorption is presented. In each of these simulations, a single gas component is selected, and the initial matrix permeability is set to k m0 =1.0e-17 m 2 . Figure 12 shows the variation of fracture aperture at detection points A and B at the different injection pressures with adsorption induced swelling effect considered. In order to analyse the influence of changing stress field on the variation of fracture aperture, the time evolutions of the stress field at  Figure 15.
It can be seen from Figure 12 that, although different injection pressures are used, the variation in fracture aperture at both measuring points follows a similar trend. Compared with the variations of fluid pressure and stress at measuring points A and B, shown in Figures 13 and 14, it can be observed that the change of fracture aperture depends mainly on the normal stress acting on the fracture surface. At the early stage, there is an increase in fracture aperture. The fracture pressure increases rapidly to reach the same as the injection pressure due to the higher permeability of fractures, leading to the decrease of effective stress normal to the fracture surfaces in fracture domain. At this stage, although the gas pressure remains very low and there is little swelling occurring in coal matrix, the effective stress normal to the fracture surfaces has a dominant role in changing the fracture aperture. The fracture aperture starts to decline rapidly after reaching a maximum value. This occurs as gas flows rapidly into the rock matrix from the fractures, with matrix swelling being localized to the vicinity of the fractured zone, especially for the higher injection pressure, because there is higher pressure gradient formed in the area closed to fractures. Gas pressure in the matrix beyond this swelling area still remains very low because no gas has reached that area, and the external boundary stays unmoved. This swelling stage is referred to as local swelling. 43,44 During local swelling, the stress field in the vicinity of the fractures behaves volumetrically, and local stress concentration results in the compression of fractures. The evolution of fracture aperture is similar to the case of a constant volume boundary. When gas diffuses into the matrix far away from the fractures and the gas pressure in the matrix continues to increase, the swelling area induced by gas adsorption extends. As the swelling area propagates to a certain size, the swelling becomes macro-swelling, the external boundary starts to move outwards, and the stress field around the fractures is controlled by the external boundary condition and behaves nonvolumetrically. 44 At this stage, the fracture aperture stops declining and starts to rebound, as shown in Figure 12.
However, there are some differences apparent in Figure 12 that can be highlighted for discussion. Higher injection pressures induce a more rapid increase of fracture aperture at the early stage, as well as leading to a larger opening of the fractures. The recovery of fracture aperture for higher injection pressure takes place earlier due to faster diffusion into the rock matrix, causing the transition of local swelling to macro swelling to occur earlier. Although the change in fracture aperture in different orientations is different, the redounding process for fractures is almost same, which implies that the stress variation is not uniform at the early period; however, it is the same when the swelling transits from local swelling to macro swelling.
As described above, fracture closure or opening is strongly dependent on the stress normal to the fracture surfaces. However, the stress change in the direction normal and parallel to the fracture is different, as shown in Figures 13 and  14. The x direction and y direction are normal and parallel to the fracture where measurement point A is located,

T A B L E 3
Simulation cases for the investigation of fracture aperture responses to injection under different conditions Case 1: Impact of injection pressure and adsorption induced swelling u in = 6MPa, ε ad = 0; u in = 8MPa, ε ad 6 ¼ 0 u in = 6MPa, ε ad 6 ¼ 0; u in = 4MPa, ε ad 6 ¼ 0 respectively, whereas x direction and y direction are parallel and normal to the fracture containing measurement point B. The effect of CO 2 adsorption-induced swelling on stress change is considered. At the early stage, the compressive stress parallel to fracture increases. Fractures are quick to be filled with the injected fluid due to their higher permeability, which leads to a larger swelling area in the direction parallel to the fractures compared with the direction normal to the fractures; the local stress concentration in the direction parallel to fractures takes place earlier. It is interesting to note that when the local stress concentration in the direction normal to fractures occurs and normal compressive stress increases, the compressive stress parallel to the fracture starts to decrease. During local swelling, the shear stress also shows a graduate increase, even exceeding the shear strength. However, due to larger magnitude of shear stiffness, the shear dilation makes little contribution to normal displacement of fracture. The maximum shear stress at detection point A is about 7 MPa, at the same time, the normal stress reaches about 16 MPa. If Mohr-Coulomb criterion is selected as failure determination, 23 even though cohesion strength is assumed to be zero, the calculated peak shear stress is still larger than the maximum. It can be concluded that the coal deformation remains elastic in these simulations. After the transition from local swelling to macro swelling, the compressive stress in both directions reaches the same value, and there is almost no shear stress at both measuring points. This implies that the stress field after reaching macro swelling varies uniformly. The impact of adsorption-induced swelling on the variation in fracture aperture is illustrated in Figure 15. For the simulation without adsorption, the fracture aperture increases with time as expected with effective stress dependency. The evolutions of fracture aperture at both detection points show similar patterns. Conversely, the variation of fracture aperture involving adsorption is relatively complex, as described above.

| Impact of different CO 2 /N 2 mixture injection
To investigate the role of the injected gas component on the variation of fracture aperture, a selection of CO 2 /N 2 gas mixture is made with different compositions, namely, pure CO 2 , 25%CO 2 :75%N 2 , 50%CO 2 :50%N 2 , 75%CO 2 :25%N 2 and pure N 2 . The injection pressure is kept constant at 6 MPa for each simulation scenario and the initial matrix permeability is set to be k m0 =1.0e-17 m 2 . Simulation results showing the evolution of fracture aperture at measuring points A and B are shown in Figures 16 and 17, respectively.
It can be seen that the fracture aperture at both locations shows a similar variation for each CO 2 /N 2 mixture injected. The effects of adsorption on fracture aperture are again illustrated. Where the gas mixture contains more CO 2 , a smaller increase of fracture aperture occurs at the early stage because its growth is supressed by the larger swelling induced by CO 2 than by N 2 , resulting in a larger local stress concentration. The change in fracture aperture caused by pure N 2 injection is close to that without adsorption, shown in Figure 15. The gas mixture containing more N 2 can also lead to a small decrease in fracture aperture, although the minimum value of aperture surpasses the initial aperture, which is different to the injection of pure CO 2 .

| Impact of matrix permeability
Three simulations with pure CO 2 injection at 6 MPa are presented that investigate the impacts of transport behaviour in the rock matrix on the evolution of fracture aperture. Three different matrix permeabilities are used, that is, k m0 =1.0e-17 m 2 , k m0 = 1.0e-15 m 2 , and k m0 = 1.0e-14 m 2 , and the simulation results are shown in Figure 18. It can be seen that the transport behaviour in the porous matrix has a significant effect on the evolution of the fracture aperture. Compared with a lower matrix permeability, there is no opening of fracture occurring for the simulation with a higher matrix permeability at the early stage. Conversely, the fracture aperture decreases at the early stage and then starts to recover to a higher value and keep constant. The flow rate in the higherpermeability matrix is larger, and the swelling area enlarges quickly; the local swelling plays a controlled role in the aperture variation over the normal effective stress from the beginning in the early period. In addition, the recovery process of fracture aperture for a higher matrix permeability occurs earlier than that for lower matrix permeability. Furthermore, the time taken to reach the final steady state is less for a higher matrix permeability. This is also because the area of matrix swelling induced by gas adsorption expands quickly, reducing the time required for transition from local swelling to macro swelling. In addition, the fracture aperture at two detection points displays no difference for high permeability medium. It is worth pointing out that although the both simulation cases with matrix permeability of 1.0e-15 and 1.0e-14 m 2 have reached the equilibrium state at the late stage, the eventual apertures of fractures at measurement point are slightly different whereas this difference is acceptable, which can be attributed to the calculation error, especially the superposition of swelling strain increment in each loop. Different time steps in each loop were used for these two simulation cases to achieve convergence of numerical solution.

| CONCLUSIONS
Naturally occurring fractures are commonly important features in candidate rock formations for carbon dioxide sequestration and unconventional gas exploration. To understand the effect of adsorptive gas injection on fracture deformation, a discrete fracture model of fully coupled compressible fluid flow, adsorption and geomechanics is presented in this work, which consists of a matrix-fracture fluid transport model, a matrix deformation model and a stress-strain model for fracture deformation. A sequential implicit discrete fracture model is presented and implemented into an existing coupled numerical modelling platform based on the Galerkin finite element method. The effect of fractures is represented explicitly through the use of lower-dimensional interface elements. The model is verified using previously obtained solutions as benchmarks. Based on the model developed in this work, the effects of injection pressure, adsorption, different CO 2 /N 2 mixture injection and rock matrix permeability are analysed. The major findings are summarized as follows: 1 The change of fracture aperture is mainly controlled by the normal stress acting on the fracture surface. Although the injection pressure may be different, the variation in fracture aperture follows a similar trend. The fracture aperture increases at the early stage due to the decrease of effective stress, after which the adsorption-induced local swelling can cause a local stress concentration, resulting in the compression of fractures. As the gas diffuses into the region far away from the fracture, the swelling area extends and the local swelling becomes macro swelling, with the fracture aperture starting to rebound. A higher injection pressure can induce a more rapid increase of fracture aperture at the early stage, as well as leading to a larger opening of the fractures. The recovery of fracture aperture at higher injection pressures also occurs earlier. 2 Compared with the complex evolution of fracture aperture involving adsorption, the impact of injecting an adsorption-free gas on the fracture aperture is monotonously increasing. The fracture aperture increases with time as expected with effective stress dependency. 3 The fracture aperture shows a similar variation at both measuring points for each CO 2 /N 2 mixture considered for injection. The gas mixture containing more CO 2 causes a smaller increase of fracture aperture at the early stage due to the fact that CO 2 adsorption can induce larger swelling than that by N 2 . The gas mixture containing more N 2 can also lead to a small decrease in fracture aperture, although the minimum value of aperture surpasses the initial aperture. 4 The transport behaviour in the porous matrix has a significant effect on the evolution of fracture aperture. Compared with a lower matrix permeability, no opening of fractures occurs for the simulation with a higher matrix permeability at the early stage. Conversely, the fracture aperture decreases and then starts to recover to a higher value and keep constant. The recovery process of the fracture aperture for a higher matrix permeability occurs earlier than that for a lower matrix permeability. In addition, the time taken to reach the final steady state is less for a higher matrix permeability.