A Mathematical Model of Gas and Water Flow in a Swelling Geomaterial—Part 1. Veriﬁcation with Analytical Solution

: Gas generation and migration are important processes that must be considered in a safety case for a deep geological repository (DGR) for the long-term containment of radioactive waste. Expansive soils, such as bentonite-based materials, are widely considered as sealing materials. Understanding their long-term performance as barriers to mitigate gas migration is vital in the design and long-term safety assessment of a DGR. Development and the application of numerical models are key to understanding the processes involved in gas migration. This study builds upon the authors’ previous work for developing a hydro-mechanical mathematical model for migration of gas through a low-permeable geomaterial based on the theoretical framework of poromechanics through the contribution of model veriﬁcation. The study ﬁrst derives analytical solutions for a 1D steady-state gas ﬂow and 1D transient gas ﬂow problem. Using the ﬁnite element method, the model is used to simulate 1D ﬂow through a conﬁned cylindrical sample of near-saturated low-permeable soil under a constant volume boundary stress condition. Veriﬁcation of the numerical model is performed by comparing the pore-gas pressure evolution and stress evolution to that of the results of the analytical solution. The results of the numerical model closely matched those of the analytical solutions. Future studies will attempt to improve upon the model complexity and investigate processes and material characteristics that can enhance gas migration in a nearly saturated swelling geomaterial.


Introduction
The primary purpose of a deep geological repository (DGR) for the long-term management of radioactive waste is to contain and isolate wastes to minimize impact to the environment and radiological exposure to people. In developing a safety case for a DGR, which provides the supportive arguments that the long-term solution for the management of radioactive waste will be protective of human health and the environment over the long term, relevant features, events, and processes (FEPs), must be evaluated [1,2]. One such process with the potential means for radiological exposure to the biosphere is the generation of radioactive gas which may migrate to the surface [3].
Gas could be generated through a number of processes including the degradation of organic matter, radioactive decay of the waste, corrosion of metals producing hydrogen gas (H 2 ), and the radiolysis of water producing H 2 [3][4][5]. If production exceeds the containment capacity of the engineered barriers or host rock, these gases could migrate through the engineered barriers and/or the host rock [6,7]. The preferential migration pathway of these radioactive gases, to potentially expose people and the environment to radioactivity, might be through the access and ventilation shafts, as these components are typically part of the repository design.
In recent years, a number of international projects have focused on the topics of gas generation and migration, with a focus on the impact of gas build-up and migration through an engineered barrier system (EBS) [3,8]. Expansive or swelling soils, such as bentonite-based materials, are currently the preferred choice of seal materials used for an EBS. Understanding the long-term performance of these seals as barriers against gas migration is an important factor in the design and long-term safety assessment of a DGR.
A wealth of laboratory and field-scale experimental studies have investigated gas transport processes through natural (host rock) and engineered barriers. These studies provide a wealth of evidence that suggest that at gas pressures above a critical level, there is formation of pressure-induced preferential flow pathways and dilation of the clay, resulting in increased gas flow. In addition, a number of mathematical models have been developed to simulate the gas transport processes observed through these laboratory and field-scale studies with some success [9][10][11][12]. However, no studies to date have been able to determine the exact mechanisms which control gas entry, flow, and pathway sealing [4,6,8,[13][14][15][16][17]. Development and use of numerical models are of key importance in understanding the processes involved and their use in long-term safety assessments.
Dagher et al. [9] developed a fully coupled hydro-mechanical (HM) linear-elastic mathematical model for advective-diffusive visco-capillary controlled two-phase flow through geomaterials in order to model the first two transport mechanisms proposed by Marschall et al. [18]. Results from a constant volume 1-D flow experiment was used to validate the model. A number of parametric studies were investigated to assess the contribution of advection of poregas, diffusion of dissolved gas in porewater, advection of dissolved gas in porewater, and inclusion of mechanical deformation (linear elasticity) on flow behaviour with increasing gas pressures over time. Additionally, sensitivity analyses were conducted to gain an understanding of the influence of a number of soil properties on flow behaviour, such as the effect of modifying the air-entry value (AEV), intrinsic permeability, and initial porosity of the soil specimen. Finally, the study investigated the use of a linear elastic damage model to better represent the experimental results. Although the model results reproduced some of the general features noted in the experimental results, the model was unable to simulate dilatancy-controlled gas flow.
In this paper, the authors build upon previous published work by first performing a verification study of the linear poro-elastic hydro-mechanical model proposed in Dagher et al. [9]. Analytical solutions for a simplified 1D steady-state and 1D transient single-phase flow through a low permeability geomaterial problem are presented. Verification of the numerical model presented in Dagher et al. [9] is then performed by comparing the results obtained by the numerical model to the results of the analytical solutions. The model is then used to investigate the effects on gas flow of different processes and phenomena, including: (i) heterogeneity within the geomaterial, (ii) the Klinkenberg "slip flow" effect on gas permeability, and (iii) swelling and dessication of the geomaterial. The results of the simulation of the above factors are presented in the companion paper to this one.
This research continues to be, in part, a contribution to Task A of the current project phase of the international working group for the DEvelopment of COupled models and their VALidation against EXperiments (DECOVALEX-2019). This task, led by the British Geological Survey (BGS), further attempts to identify the physical HM mechanisms required to adequately model dilatancy-controlled gas migration.

Study Overview
This study expands upon the work performed by Dagher et al. [9] on the development of a mathematical model for gas migration (two-phase flow) through a low-permeable swelling geomaterial. Using the theoretical framework of poromechanics, a verification study is performed through the following: The mathematical model follows the general formulation by Dagher et al. [9]. The applicable constitutive relations and governing equations for conservation of momentum, water mass and gas mass are presented below.

Assumptions for Poromechanical Behaviour
This paper adopts the following assumptions for the mechanical behaviour of the porous medium as presented in Dagher et al. [9]: Bishop's modified effective stress, with a χ parameter generalized from the work of Khalili and Khabbaz [19], linear poro-elasticity following the general framework of poromechanics.

Bishop's Modified Effective Stress Principle
The general form for the equation of poro-mechanics can be expressed by, where σ ij is total normal stress tensor acting on the soil element (Pa) σ ij is the effective stress tensor (Pa) δ ij is the Kronecker delta (identity tensor) (adimensional) and p is the porefluid pressure (Pa) Many effective stress equations have been proposed to characterize the stress-state of an unsaturated soil or porous media [20]. This paper proposes the use of Bishop's effective stress principle, which is dependent on both net normal stress and matric suction, and may be more suitable for expansive clays, and is described by Equation (2).
where p g is the poregas pressure (Pa) p w is the porewater pressure (Pa) σ ij − p g is the net normal stress (Pa) p g − p w is the matric suction (Pa) χ is a parameter related to the degree of saturation of the soil (unitless) Expanding and rearranging for σ, The porefluid pressure can be defined as, Khalili and Khabbaz [19], proposed the following unique relationship for the determination of χ based on the ratio of suction over the air entry value, also termed the suction ratio, where p g − p w b is the air-entry value (AEV) of the soil and only applies when the matric suction >AEV [19].

Poro-Elastic Stress-Strain Relation
The increment of the effective stress tensor, dσ ij is related to the increment of the total strain by the constitutive relation: dσ ij = C ijkl dε kl (6) where C ijkl is the stiffness tensor (Pa). ε kl is the total strain tensor (adimensional). Assuming small deformations, the total strain is related to the components of the displacement tensor as: where ε kk = u k,k = trε ij = ε vol is the volumetric strain (adimensional).
is the spatial derivative of the displacement vector (m).

Constitutive Relations for the Hydraulic Behaviour
This paper adopts the following constitutive relations for the hydraulic behaviour as presented in Dagher et al. [9]: • van Genuchten equation for the soil-water characteristic curve [9]. • Huang et al. [21] equation depicting the relationship between the AEV and void ratio. • Darcy's law for two-phase flow. • Equation for the effective diffusivity for gas dissolved in water through porous media, utilizing the Millington and Quirk model to define the tortuosity as a function of the degree of saturation, S w and the porosity, n [22].

•
Pall and Moshenin et al. [23] modification to the Kozeny-Carman equation [24] for the relationship between the intrinsic permeability of water and the porosity of the soil. • Mualem's model for the relative permeabilities of gas and water [9,20].

Governing Equations for Two-Phase Flow through a Linear Elastic Geomaterial
The governing equations as derived by Dagher et al. [9] and Nguyen and Le [20] for describing advective-diffusive visco-capillary controlled two-phase flow for a linear-elastic geomaterials are presented below. The governing equation for the conservation of water mass can be expressed by Equation (9), and solving the right-hand side, this simplifies to Equation (10), where ρ w density of water phase (kg m −3 ) p w is the porewater pressure (Pa) k ij intrinsic permeability tensor of the porous medium (m 2 ) k r,w relative permeability of the water phase (unitless) µ w dynamic viscosity of the water phase (Pa s or kg m −1 s −1 ) g is the acceleration due to gravity (m s −2 ) n porosity (m 3 voids · m −3 total) S w is the degree of saturation of water s is the matric suction p g − p w (Pa) K w is the bulk modulus of the water phase (Pa s or kg m −1 s −1 ) Note that in this study, the permeability is assumed to be isotropic, therefore k ij = k, however we keep the tensorial notation in the governing equation for the sake of generalization.

Conservation of Gas Mass
The governing equation for the conservation of gas mass can be expressed by Equation (11), and solving the right-hand side, this simplifies to Equation (12), where ρ g density of the gas phase (kg m −3 ) p g is the poregas pressure (Pa) k r,g relative permeability of the gas phase (unitless) µ g dynamic viscosity of the gas phase (Pa s or kg m −1 s −1 ) H is Henry's coefficient (kg species A m −3 in aqueous phase kg −1 species A m 3 in gas phase) D e is the effective diffusivity of gas dissolved in water through porous media (m 2 s −1 ) K g is the bulk modulus of the gas phase (Pa) It should be noted that Equation (12) does not require the assumption that the gas behaves as an ideal gas. In the derivation of the expression of gas density and bulk modulus one might need to invoke the ideal gas law or other laws that more accurately predict the gas density and compressibility at different pressures and thermal conditions.

Conservation of Momentum (Quasi-Static Equilibrium)
For a soil specimen in quasi-static equilibrium the total equilibrium equation for a cubical soil element can be expressed in tensor notation by Equation (18).
where, σ ij is the stress tensor (Pa) F v,i is the volumetric body force tensor (kg m −2 s −2 ) ∂σ ij ∂x j represents the change in normal and shear stresses across the soil element (kg m −2 s −2 ) Substituting Equation (3) into Equation (13), the governing equation for the conservation of momentum for a porous geomaterial can be expressed by Equation (14), For a linear poro-elastic geomaterial, applying the constitutive relations described by Equations (6) through (8), the governing equation for the conservation of momentum can be expressed by Equation (15), where G is the shear modulus (Pa), λ is the Lamé constant (Pa).

Analysis Overview
Verification of the proposed mathematical model was performed by deriving analytical solutions for the 1D steady-state and 1D transient gas flow problem depicted in Figure 1, through the simplification of the governing equations described by Equations (1) and (3). In the above problem, flow and deformation can occur only in the longitudinal direction x. A three-dimensional (3D) hydro-mechanical (HM) coupled multiphysics numerical model was built using the commercially available code COMSOL ® Multiphysics ® (COMSOL ® ) to solve the simplified governing equations and constitutive relations described below numerically. The analytical solution was then used to verify the numerical model in COMSOL ® by comparing the numerical results to the results of the analytical solutions. flow and deformation can occur only in the longitudinal direction x. A three-dimensional (3D) hydromechanical (HM) coupled multiphysics numerical model was built using the commercially available code COMSOL ® Multiphysics ® (COMSOL ® ) to solve the simplified governing equations and constitutive relations described below numerically. The analytical solution was then used to verify the numerical model in COMSOL ® by comparing the numerical results to the results of the analytical solutions.  The following assumptions were made for the simplification of the mathematical model: It should be noted that these assumptions, as well as the downstream boundary condition set to atmospheric pressure ( Figure 1) was done for simplification of the verification analysis. In a bentonite-based shaft seal design, this may be a reasonable assumption when one end of the shaft seal is exposed to surface. However, if applied to bentonite buffer surrounding the waste canister, the boundary pressure conditions would be different and would affect the fluid properties at the higher temperatures and pressures of a DGR environment.

One-Dimensional (1D) Steady-State Gas Flow Problem
Based on the above assumptions, the governing equation for the conservation of gas mass for 1D steady-state flow through a linear elastic porous medium expressed by Equation (11) simplifies to the following expression, d dx where, k g is the gas permeability (m 2 ) Assuming an immobile water phase and ignoring the effects of gravity, the governing equation for the conservation of momentum expressed by Equation (14), simplifies to the following equation for a 1D gas flow problem, dσ dx Solving Equation (16), the analytical solution for the pore-gas pressure and displacement as a function of distance can be expressed by, where P inj is the gas injection pressure (Pa) P atm is atmospheric pressure (Pa) L is the length of the soil specimen (m) The displacement field along the longitudinal x-direction is given by: where ν is the Poisson ratio (unitless) E is Young's Modulus (Pa) The corresponding strain and stress can be calculated by the following expressions, Detailed derivation of the analytical solution for the 1D steady-state gas flow problem is provided in Appendix A.

1D Transient Gas Flow Problem
Based on the assumptions, the governing equation for the conservation of gas mass for 1D transient gas flow through a linear elastic porous medium expressed by Equation (11) simplifies to the following, By assuming the gas density and degree of saturation remains constant in both time and space, the governing equation for the conservation of momentum expressed by Equation (19) simplifies to the following, Solving for dn, this simplifies further to, Since strain is nil in the y and z directions, it can be shown from Hooke's Law that Minerals 2020, 10, 30 Substituting into Equation (25) and re-arranging, Integrating both sides of the governing equation for the conservation of momentum expressed by Equation (17), and re-arranging for σ xx , As Equation (29) is analogous to the equation for transient diffusion in a plane sheet, the analytical solution as expressed by J. Crank [25] for constant surface concentrations and constant initial conditions is expressed by the following equation, (1−S w )(1−χ) (m 2 s −1 ) Substituting p g (x,t) from Equation (30) into Equation (28) and then substituting Equation (28) into Equation (26) and integrating both sides, the analytical solution for the displacement can be expressed by, The strain and stress can now be calculated analytically by the following expressions, Detailed derivation of the analytical solution for the 1D transient gas flow problem is provided in Appendix A.

Mean Stress and Deviatoric Stress
The mean stress can be calculated by the following expression, The deviatoric stress can be calculated by the following expression, where,

Model Geometry and Mesh
For the numerical model, the geometry and mesh are presented in Figure 2. The numerical model consists of approximately 11,000 elements and has 162,000 degrees of freedom.

Model Geometry and Mesh
For the numerical model, the geometry and mesh are presented in Figure 2. The numerical model consists of approximately 11,000 elements and has 162,000 degrees of freedom.

Results and Discussion
Verification of the numerical model results against the analytical solution for a 1D steady-state gas flow problem and a 1D transient gas flow problem are presented below. Table 1 provides the material properties, initial conditions, and boundary conditions that were used in both the numerical model and analytical solution. The material properties correspond to those obtained experimentally from Daniels and Harrington [26] or calculated using the constitutive relations derived by Dagher et al. [9]. The initial and boundary conditions were selected based on the experimental setup of Daniels and Harrington [26], with the injection poregas pressure corresponding to the peak injection pressure observed in that experiment.

Results and Discussion
Verification of the numerical model results against the analytical solution for a 1D steady-state gas flow problem and a 1D transient gas flow problem are presented below. Table 1 provides the material properties, initial conditions, and boundary conditions that were used in both the numerical model and analytical solution. The material properties correspond to those obtained experimentally from Daniels and Harrington [26] or calculated using the constitutive relations derived by Dagher et al. [9]. The initial and boundary conditions were selected based on the experimental setup of Daniels and Harrington [26], with the injection poregas pressure corresponding to the peak injection pressure observed in that experiment.

Steady-State Solution
For the 1D steady-state gas flow problem, the poregas pressure and displacement in the x-direction were calculated using analytical solution presented by Equations (18) and (19), respectively. The numerical and analytical results for poregas pressure and displacement over the length of the specimen is provided in Figures 3 and 4, respectively. The results show a perfect match between the simulation results of the numerical model performed in COMSOL ® and the analytical solution.

Steady-State Solution
For the 1D steady-state gas flow problem, the poregas pressure and displacement in the xdirection were calculated using analytical solution presented by Equations (18) and (19), respectively. The numerical and analytical results for poregas pressure and displacement over the length of the specimen is provided in Figures 3 and 4, respectively. The results show a perfect match between the simulation results of the numerical model performed in COMSOL ® and the analytical solution.     Similarly, the mean effective stress, p , and deviatoric effective stress, q , were calculated by first calculating the strain and resulting principal stresses described by Equations (20) to (22), and then Equations (33) and (36). The numerical and analytical results for mean effective stress, and deviatoric effective stress over the length of the specimen are provided in Figures 5 and 6, respectively. The results, also show a perfect match between the numerical results performed in COMSOL ® and the analytical solution, thereby demonstrating verification of the numerical model being used in COMSOL ® . Similarly, the mean effective stress, p′, and deviatoric effective stress, q′, were calculated by first calculating the strain and resulting principal stresses described by Equations (20) to (22), and then Equations (33) and (36). The numerical and analytical results for mean effective stress, and deviatoric effective stress over the length of the specimen are provided in Figures 5 and 6, respectively. The results, also show a perfect match between the numerical results performed in COMSOL ® and the analytical solution, thereby demonstrating verification of the numerical model being used in COMSOL ® .

Transient Solution
For the 1D transient gas flow problem, the poregas pressure and displacement were calculated using the analytical solution presented by Equations (30) and (31). The numerical and analytical

Transient Solution
For the 1D transient gas flow problem, the poregas pressure and displacement were calculated using the analytical solution presented by Equations (30) and (31). The numerical and analytical results for poregas pressure and displacement over the length of the specimen for various times are provided in Figures 7 and 8, respectively. The results show a perfect match between the simulation results of the numerical model performed in COMSOL ® and the analytical solution. For the poregas pressure, there is an initial response at t = 0 s, corresponding to the sudden change in boundary conditions. Over time, the system tends towards the steady-state solution. Concurrently, as the poregas pressure increases along the column, there is a resulting increase in displacement, until steady-state is reached. Steady-state was achieved at about 5000 s, and matches the solutions provided in Section 3.4.1.
It should be noted that oscillations occurred in results obtained from both the analytical solution and numerical model at time t = 0 s, and this was expected. For the analytical solution, these oscillations are due to truncation of the infinite series at n = 70 and m = 40. As these parameters tend towards infinity, if higher truncation limits were used, the oscillations would effectively dissipate. For the numerical analysis, these oscillations are an artefact of mesh size, time step, and error tolerance. Increasing mesh size, and reducing the time-step and tolerance, would help to reduce the numerical oscillations.
Similarly, the mean effective stress, p , and deviatoric effective stress, q , were calculated by first calculating the strain described by Equation (32) and resulting principal stresses described by Equations (26) and (22), and then by Equations (33) and (36). The numerical and analytical results for mean effective stress, and deviatoric effective stress over the length of the specimen are provided in Figures 9  and 10, respectively. The results, also show a perfect match between the numerical results performed in COMSOL ® and the analytical solution, thereby demonstrating verification of the numerical model being used in COMSOL ® . and numerical model at time t = 0 s, and this was expected. For the analytical solution, these oscillations are due to truncation of the infinite series at n = 70 and m = 40. As these parameters tend towards infinity, if higher truncation limits were used, the oscillations would effectively dissipate. For the numerical analysis, these oscillations are an artefact of mesh size, time step, and error tolerance. Increasing mesh size, and reducing the time-step and tolerance, would help to reduce the numerical oscillations.  Similarly, the mean effective stress, p′, and deviatoric effective stress, q′, were calculated by first calculating the strain described by Equation (32) and resulting principal stresses described by Equations (26) and (22), and then by Equations (33) and (36). The numerical and analytical results for mean effective stress, and deviatoric effective stress over the length of the specimen are provided in Figures 9 and 10, respectively. The results, also show a perfect match between the numerical results performed in COMSOL ® and the analytical solution, thereby demonstrating verification of the  Additionally, an analysis of the p-q stress path at five different points along the specimen over time was conducted and the results are presented in Figure 11. Again, the numerical model results  Additionally, an analysis of the p-q stress path at five different points along the specimen over time was conducted and the results are presented in Figure 11. Again, the numerical model results Additionally, an analysis of the p-q stress path at five different points along the specimen over time was conducted and the results are presented in Figure 11. Again, the numerical model results and results of the analytical solution show a close match. The analysis shows that at the center of the specimen, there is no change in stresses as expected. The analysis also shows that at the two ends, there is also no change in stresses as described by the analytical solution, however, the numerical model does show a small change in stress from 0 s to t > 0 s. Again, this can likely be attributed to the truncation of the infinite series. At distances 0.03 m and 0.09 m, the effective deviatoric stress increases linearly as a function of the effective mean stress until steady-state is reached. This is expected as a linear elastic mechanical model was adopted. It should be noted that at the injection pressure side, the mean effective stress is in tension, while at the atmospheric pressure side, the mean effective stress is in compression. and results of the analytical solution show a close match. The analysis shows that at the center of the specimen, there is no change in stresses as expected. The analysis also shows that at the two ends, there is also no change in stresses as described by the analytical solution, however, the numerical model does show a small change in stress from 0 s to t > 0 s. Again, this can likely be attributed to the truncation of the infinite series. At distances 0.03 m and 0.09 m, the effective deviatoric stress increases linearly as a function of the effective mean stress until steady-state is reached. This is expected as a linear elastic mechanical model was adopted. It should be noted that at the injection pressure side, the mean effective stress is in tension, while at the atmospheric pressure side, the mean effective stress is in compression.

Conclusions
An important component in the design and long-term safety assessment of a DGR is the longterm performance of bentonite seals as barriers against gas migration. As gas generates from the degradation of organic waste and/or corrosion of metals, at some critical gas pressure, dilation of the bentonite could occur resulting in the creation of preferential flow pathways and a source of radionuclide exposure to people and the environment. Development and use of numerical models are of key importance in the understanding of processes involved and in their use in long-term safety assessments.
In the authors' previous work, a fully-coupled hydro-mechanical (HM) linear-elastic mathematical model for advective-diffusive visco-capillary controlled two-phase flow through a geomaterial was derived [9]. This study expands upon that work through the contribution of model verification. As part of this study, analytical solutions for a 1D steady-state gas flow and a 1D transient gas flow problem were derived using the framework of poromechanics. Using the finite element method, the numerical model was used to simulate 1D flow through a confined cylindrical sample of near-saturated low-permeable soil under a constant volume boundary stress condition. Verification of the numerical model was performed by comparing the pore-gas pressure evolution and stress evolution to that of the results of the analytical solutions. For both the steady-state and

Conclusions
An important component in the design and long-term safety assessment of a DGR is the long-term performance of bentonite seals as barriers against gas migration. As gas generates from the degradation of organic waste and/or corrosion of metals, at some critical gas pressure, dilation of the bentonite could occur resulting in the creation of preferential flow pathways and a source of radionuclide exposure to people and the environment. Development and use of numerical models are of key importance in the understanding of processes involved and in their use in long-term safety assessments.
In the authors' previous work, a fully-coupled hydro-mechanical (HM) linear-elastic mathematical model for advective-diffusive visco-capillary controlled two-phase flow through a geomaterial was derived [9]. This study expands upon that work through the contribution of model verification. As part of this study, analytical solutions for a 1D steady-state gas flow and a 1D transient gas flow problem were derived using the framework of poromechanics. Using the finite element method, the numerical model was used to simulate 1D flow through a confined cylindrical sample of near-saturated low-permeable soil under a constant volume boundary stress condition. Verification of the numerical model was performed by comparing the pore-gas pressure evolution and stress evolution to that of the results of the analytical solutions. For both the steady-state and transient problems, the poregas pressure, displacement, mean effective stress, deviatoric effective stress, and p'-q' stress path results obtained from the numerical model were compared to those of the analytical solutions. For each parameter, the results of the numerical model closely matched those of the analytical solutions, providing added confidence that the proposed mathematical model is being correctly implemented in the COMSOL ® FEM software version 5.4.
Future studies will attempt to improve upon the model complexity and investigate enhanced mechanisms for two-phase flow. In the companion paper, a process simulation and enhanced two-phase flow analysis study will be conducted, whereby three enhanced mechanisms for two-phase flow will be introduced into the model. Specifically, the companion paper will look at heterogeneity, the Klinkenberg "slip flow" effect, and the effect of swelling and dessication. An analysis of the contribution of each to flow behaviour and the potential for the formation of gas fingers will be conducted.   Taipower, for their financial and technical support of the work described in this paper. The statements made in the paper are, however, solely those of the authors and do not necessarily reflect those of the funding organizations.

Conflicts of Interest:
The authors declare no conflict of interest.

A.1. Problem Statement
Derive the steady-state and transient solutions for 1D gas flow through a linear elastic porous media of length L, fixed at both ends, u(0,t) and u(L,t) = 0, with an initial uniform poregas pressure, p g (x,0) = p 0 , and with both maintained at different constant gas pressures, p g (0,t) = p inj and p g (L,t) = p atm , Minerals 2019, 9, x FOR PEER REVIEW 17 of 29 transient problems, the poregas pressure, displacement, mean effective stress, deviatoric effective stress, and p'-q' stress path results obtained from the numerical model were compared to those of the analytical solutions. For each parameter, the results of the numerical model closely matched those of the analytical solutions, providing added confidence that the proposed mathematical model is being correctly implemented in the COMSOL ® FEM software version 5.4. Future studies will attempt to improve upon the model complexity and investigate enhanced mechanisms for two-phase flow. In the companion paper, a process simulation and enhanced twophase flow analysis study will be conducted, whereby three enhanced mechanisms for two-phase flow will be introduced into the model. Specifically, the companion paper will look at heterogeneity, the Klinkenberg "slip flow" effect, and the effect of swelling and dessication. An analysis of the contribution of each to flow behaviour and the potential for the formation of gas fingers will be conducted.

Conflicts of Interest:
The authors declare no conflict of interest.

A.1. Problem Statement
Derive the steady-state and transient solutions for 1D gas flow through a linear elastic porous media of length L, fixed at both ends, u(0,t) and u(L,t) = 0, with an initial uniform poregas pressure, pg(x,0) = p0, and with both maintained at different constant gas pressures, pg(0,t) = pinj and pg(L,t) = patm, Figure A1. One-dimensional (1D) gas flow through a porous media of length, L. Figure A1. One-dimensional (1D) gas flow through a porous media of length, L.
The following assumptions were made for the simplification of the mathematical model: The governing equation for the conservation of momentum (quasi-static equilibrium) for steady-state gas-flow through a linear elastic porous media can be expressed by Equation (A2), Assuming an immobile water phase, this simplifies to Equation (A3),

A.2.2. Formulation of the Analytical Solution to the Conservation of Mass Equation
Re-arranging Equation (A1) and solving the double integral for p g , From our known boundary conditions, at x = 0, p g (0) = P inj and we can solve for C 1 Substituting (A6) back into (A5) and from our known BCs at x = L, p g (L) = P atm , we can solve for C 1 P atm = C 1 (L) + P inj → C 1 = P atm − P inj L (A7) Substituting Equations (A6) and (A7) into Equation (A5), the analytical solution for the governing equation for the conservation of momentum can be expressed by Equation (A8)

A.2.3. Formulation of the Analytical Solution to the Conservation of Momentum Equation
From Hooke's Law the normal principal strains and shear strains are provided in Equations (A9) to (A11), and Equations (A12) to (A14) respectively, For the 1D flow problem, ε yy and ε zz = 0; v and w = 0; and du dy = 0 and du dz = 0, these equations reduce to the following, Re-arranging Equations (A16) and (A17), Equating Equation (A21) and Equation (A22), Substituting Equation (A26) into Equation (A21) and re-arranging for σ yy as a function of σ xx , Substituting (A27) and (A26) into Equation (A15), we can solve for ε xx , Re-arranging in terms of σ xx , where for a 1D flow problem, it can be shown from Hooke's Law that and the Bulk modulus for 1D flow problem, K 1D , can be expressed as, Substituting (A4) and (A31) into (A3) yields, Re-arranging, d 2 u Solving the double integral, From our known BCs, at x = 0, u(0) = 0 and we can solve for C 2 From our known BCs, x = L, u(L) = 0 and we can solve for C 2 Substituting Equations (A38) and (A39) into Equation (A37), the analytical solution for the governing equation for the conservation of momentum can be expressed by Equation (A40) and simplified further to Equation (A41)

A.2.4. Final Analytical Solution
The final analytical solution for 1D steady-state gas flow through a linear elastic porous media can be expressed by Equations (A42) and (A43), The strain and stress can now be calculated analytically by the following expressions,

A.3.1 Governing Equations
The governing equation for the conservation of gas mass for transient (i.e., time-dependent) flow through a linear elastic porous medium can be expressed by Equation (A47), d dx By assuming the gas density and degree of saturation remains constant Equation (A47) simplifies to Equation (A48), and solving for dn, this simplifies further to, The governing equation for the conservation of momentum (quasi-static equilibrium) for steady-state gas-flow through a linear elastic porous media assuming an immobile water phase was previously described by Equation (A52) as, Integrating both sides and re-arranging for σ ,

A.3.2 Formulation of the Analytical Solution to the Conservation of Mass Equation
The form of Equation (A53) is analogous to that of non-steady state diffusion in a plane sheet presented in Crank [25]. The analytical solution for non-steady state diffusion in a plane sheet as derived by Crank [25] for constant surface concentrations and constant initial conditions is expressed by the following equation, C 2 cos nπ−C 1 n sin nπx L exp − Dn 2 π 2 t L 2 where C(x, t) is the concentration at a given distance, x, and time, t, along the plane sheet C 1 is the concentration at boundary x = 0 C 2 is the concentration at boundary x = L C 0 is the initial uniform concentration Considering our boundary and initial conditions provided in the problem statement.