Modelling of unsaturated gas flow by Thebes code: Validation tests

. This paper presents the simulations replicating two well-documented benchmarks on coupled gas-liquid flow in unsaturated soil. The results serve as validation and verification of the formulation of the gas flow in unsaturated geomaterials in the newly developed THMC coupled FEM code Thebes. The paper first discusses the basis of the compositional method and the role of the dry air mass balance equation in the theoretical framework. The fundamental constitutive assumptions related to gas flow, as adopted in the Thebes code, are also discussed in details. Afterwards, the paper discusses simulation of a two-phase infiltration test in unsaturated sand, as well as a one dimensional drainage test. The numerical results of these two examples show that the code is able to capture the main features associated with the gas flow in unsaturated soil. The possible future improvements, both related to the theoretical framework and the numerical implementation, are discussed at the closure of the paper.


Introduction
Bentonite as a sealing material for nuclear waste has to withstand severe environmental conditions. In such engineering applications, the assessment of bentonite performance requires sufficient understanding of the influence of the dominant coupled thermo-hydromechanical-chemical (THMC) processes on the bentonite behaviour. In addition, the analyses should consider the effect of the gas that could be potentially generated due to the long term corrosion of the spent nuclear waste canisters. The generated gas pressure, may change the whole THMC system as it may affect the pore structure and consequently the thermal conductivity, hydraulic conductivity and mechanical response. Also, depending on the gas components, additional chemical reactions could be triggered. Therefore, the capability of modelling gas flow is necessary for any code that is designed to be used for simulating nuclear waste sealing related problems.
In fact, the importance of modelling gas flow and migration is not restricted to nuclear waste disposal applications but covers several other equally serious problems, including, for example design of the clay barriers to control the gas migration accompanying the degradation of organic matters in case of landfills.
This paper presents verification of the capabilities of Thebes code related to the coupled effects of the gas flow on the thermo-hydro-mechanical-chemical behaviour of geomaterials. The paper introduces first the adopted mass balance equations and their numerical implementation, which are followed by the two validation examples.

Mass balance
In Thebes code, the theoretical derivation of mass balance equations for each component or specie is based on the compositional method [1]. In the current version, the code assumes that the porous medium is composed of water, air, solid particles and salt. These components can be found in different states. For example, water is found as liquid and as vapour in the gas phase whereas besides the dry air in the pores, dissolved dry air can be found in the liquid phase. The detailed derivation of the balance equations in Thebes code can be found in [2,3]. For the purposes of this paper, only water and dry air mass balance equations are presented in the following sections.

Water mass balance equation
The mass balance equation of the water component is given as follows [3]: where is the porosity, and are liquid and gas degree of saturation, respectively. The condition + = 1 should be always satisfied. The densities of liquid water and water vapour are denoted by and . The fluid advective flow is assumed to obey Darcy's law through the Darcian velocities and for liquid and gas phases, respectively. The non-advective flux of vapour in the gas phase ( ) is estimated in accordance with Philips and De Vries theory [4]. The mechanical coupling is accounted for through the variation of porosity and volumetric strain ! over time .
Solution of Equation (1) requires constitutive models to estimate mechanical effects, the retention properties, the hydraulic conductivity and diffusivity of liquid and gas in the studied geomaterial. Assumptions related to thermal flow and salt transport are also needed if the tackled problem is non-isothermal and chemical effects to be studied. For a deeper look on the derivation and the adopted constitutive assumptions, the reader is referred to [2,3].

Dry air mass balance equation
The derivation of dry air mass balance equation follows from the sum of dry air mass in the geomaterial voids and the mass of dissolved dry air in the liquid water. The corresponding equation is written as [3]: This equation considers similar transport mechanisms to that for the water component. In addition to the quantities identified in Equation (1), this equation requires an estimation of Henry's volumetric coefficient of solubility / to be able to calculate the mass of the dissolved dry air in the liquid water. In the formulation, it is assumed that the density of dry air in the gas phase . equals the dissolved dry air density in the liquid phase . and they both denoted simply by . in Equation (2).
The estimation of . = . follows from the ideal gas law where the total gas density is given as: and In the above formulae the total gas density is the sum of the densities of its components, namely water vapour and dry air . . The symbols 4 refers to the reference water vapour density at temperature H. The symbols 3, E . , E and F denote the universal gas constant, the molar mass of dry air, molar mass of vapour and gravity acceleration. The relative humidity 3/ is a direct function of the suction head I defined as the difference between the pore gas pressure head ℎ and the pore liquid pressure head ℎ .

Numerical implementation
The coupled thermo-hydro-mechanical-chemical balance equations are discretised spatially using the finite element method [3]. This results in a set of first order differential equations of the form: where J and M are matrices that corresponds to the governing balance equations. The vector N collects the externally applied THMC loads including that which stay constant or changing over time. The dot • means derivation with respect to time. The vector K contains the primary unknowns to solve for: where P Q, R S T , R S U , V W and X̂ are the nodal values of displacement, liquid pressure head, gas pressure head, temperature and concentration, respectively. Depending on the processes involved in the modelled problem, the required balance equations to be solved can be active or not. For the solution over time, these equations are discretised further with the finite differences method transferring them to a set of nonlinear algebraic equations. These equations are solved iteratively for the time step [ + 1, knowing the initial conditions at time step [, using Newton-Raphson method after building the following residuals \: The iterations continuous until the relative error bδK d^_ e^_ b/bK d^_ e^_ b reaches a value less than a prescribed tolerance, where δK d^_ e^_ = K d^_ e^_ − K d e^_ is the improved increment of K ]^_ value at iteration g + 1 and time step In what follows a tolerance of 1.0x10 -4 is adopted for displacement, gas and water pressure head calculations. This is somehow tight tolerance but at the same time is a good test for the code convergence. Note that the code allows for different tolerances to be assigned for different nodal primary variables, when necessary.

Infiltration into bounded sand column (coupled air-water flow) [5,6]
Touma and Vauclin [5] carried out infiltration tests on partially saturated sandy soil. The purpose was to study the effect of different boundary conditions with respect to the air flow, on the simultaneous infiltration process of water. The bottom boundary of the studied sand column was open in one set of tests and closed in the other. The water infiltration prescribed at the top boundary varied between constant head and constant flux over the time. The tests showed clearly that the development of pore air pressure had significantly affected the infiltrated water movement.
In this paper, only one experiment will be used for the validation purposes, where the sand column was closed at the bottom boundary for both water and air flow. The top boundary was subject to a constant water infiltration rate over time but kept open to air flow. The reason for choosing this experiment is that would allow to see the code capability in capturing the development of pore air pressure and also how that would affect the infiltrated water profile. In addition, previous numerical simulation of this test is already available in literature [6] which can be used as an additional reference for comparison.

Geometry and soil properties
The infiltration test was performed on sand column with a total height of 1.0 m. The bottom boundary was kept closed for flow whereas the top boundary was open for air flow and in the same time was subject to a constant water infiltration rate of q w =8.3 cm/h (2.3x10 -5 m/s). The initial suction head had hydrostatic distribution with prevailing atmospheric air (gas) pressure head h a = 0.0 kPa. No mechanical effects are considered or measured in the test and therefore the porosity stays constant throughout the analysis at its measured initial value = 0.312.
The measured water retention curve is fitted with van Genuchten formula [7] which gives the water degree of saturation as: where the suction head I is estimated in meters. The constants F o , F s and F w are fitting parameters. The symbols j.k and lmj are the degree of saturation at full and residual saturation state, respectively. For the current analysis, the fitting produces the following values: where the water permeability } is estimated in [m/s] and • = is the volumetric water content. Equation (10) yields a value of } j.k = 42.82 × 10 *€ m/s at full saturation. As for air flow, the degree of air saturation is calculated from the condition + = 1 at any time step of the test. The air permeability } was given by Touma and Vauclin [5] as: where } •l‚ = 7.778 × 10 *6 m/s is the air permeability at the fully dry state and the fitting parameters A and B have the values ƒ = 3.86 × 10 * † and ‡ = −2.4. As the problem is isothermal, a constant value of H = 0.0135 is used for Henry's coefficient. The total duration of the experiment is 6000 seconds (1.667 hours).  Figure 2 shows that the code is able to replicate the experimental measurements in terms of volumetric water content at different time steps of the test. More interesting, for the purposes of this paper, are the results presented in Figure 3 that show a comparison to the reference solution given by Celia and Binning [6] for exactly the same benchmark. As can be seen, the code is able to duplicate the results for the development of water pressure head and gas (air) pressure head evolution over time. Tiny differences in the predicted pressure heads is noticed at the late stages of the test. That could be explained by the different expressions used to estimate the density of gas in Thebes code and in the work in [6]. Celia and Binning assumed that the gas (air) density follows the following formula:

Numerical results
where 4 and ℎ 4 are reference gas (air) density and gas (air) pressure head, respectively. On the other hand, Thebes code uses more complicated formula to estimate the total gas density as discussed in Section 2.2. This has several consequences on the final form of the mass balance equation of the dry air component, considering the derivatives with respect to time for the involved components densities. Furthermore, the liquid water density is not constant in Thebes code but a function of pore water pressure, temperature and salt concentration. The differences in the results in this example is small as the problem is isothermal which cuts out the effect of temperature on vapour generation. The full discussion on this issue is out of the scope of this paper, however, the interested reader is referred to [3] for deeper details.  For the sake of completeness, Figure 4 shows a comparison between the predicted pore water pressure head profile in the case of coupled (two-phase flow) versus uncoupled analysis (single phase flow). The comparison illustrates how the developed gas pressure hindered the infiltration of water and reduced the pore space available for water, resulting in development of considerable positive pore water if compare to the single-phase analysis with the assumption of zero pore air pressure. The effect gets more obvious with the development of time and increasing gas pressure values.

Drainage of sand column (coupled air-watermechanical behaviour) [8,9]
In this example, the mechanical effect is taken into account when studying the coupled air-water flow. The simulation replicates the one dimensional drainage test on sand column as conducted by Liakopoulos [9] and uses as a reference, a benchmark problem by Jommi et al. [8] for testing the performance of different coupled THM codes. The hydraulic properties of the sand for both water and air flow are well described by the original work of Liakopoulos, however, the mechanical properties were not provided. To assess the mechanical effects, linear elasticity was assumed by Jommi et al., therefore, linear elasticity is also assumed in what follows to allow for meaningful comparison with the numerical results provided in [8].

Geometry and soil properties
The experiment is carried out on an initially fully saturated sand column with initial hydrostatic distribution of positive pore water pressure and zero air pressure. The total height of the column is 1.0 m. The measurements start when the bottom boundary of the column is made open for air and water flow resulting in drainage of water under gravitational flow. According to Jommi et al. [8], the sand has an initial porosity = 0.29, elastic modulus ‹ = 1300 gOE•, Poisson's ratio ν = 0.4 and density of solid particles j = 2000 gF/y 6 . Despite of the relatively low values of the assumed elasticity modulus and particle density for a sand, they have been adopted in this analysis to enable the comparison. The water retention curve is given as: The water permeability function is: and gas permeability function: where the effective degree of saturation m = − lmj / j.k − lmj with j.k = 1.0 and lmj = 0.2. The saturated water permeability } j.k = 4.4127 × 10 *€ m/s and the air permeability at fully dry state } •l‚ = 2.4515 × 10 *< m/s. Also following the modelling specifications provided by Jommi et al., the air permeability is constrained with a lower bound value of } w's = 10 *< × } •l‚ to enforce the continuity of gas phase once the sand approaches full water saturation state. Figure 5 shows the finite element mesh, mechanical, water flow and gas flow boundary conditions as used to simulate the experiment. The finite element mesh consists of 75 quadrilateral 4-noded finite elements with 4 Gauss integration points per element. The number of elements are kept close to that used by the contributors in Jommi et al. [8]. Nevertheless, a separate study showed that increasing the refinement of the mesh do not yield noticeable change in the results.

Finite element model
For mechanical balance solution and to consider the effect of unsaturation during elasticity, a special version of the code that operates with Bishop's effective stresses is used where ' " = ' − " + " − " with ', ' " , " and " are the total stress, Bishop's effective stress, pore gas pressure and pore water pressure, respectively. This is done to keep consistent with Jommi et al. [8] work and for the sake of comparison. It is worth mentioning, however, that Thebes code in its standard version uses the two independent stress variables concept, namely net stress ' − " and suction " − " , during mechanical balance calculations.
The initial stress field is generated using gravity loading and the resulted displacements are set to zero before the experiment simulation starts. The total duration of the experiment is 120 minutes which is divided into 1800 time steps (∆t = 4.0 s). Figure 6 depicts the pore water pressure head evolution over the experiment time as predicted by Thebes code and compared to the reference numerical solution provided by Jommi et al. [8]. The two solution match perfectly but they both deviate considerably from the measured data by Liakopoulos [9] after 30.0 min of the onset of the physical test. That relates directly to the uncertainty in the estimation of the mechanical properties of the material and the assumption of linear elasticity. The authors expect that better match to the experimental data would be possible if more accurate mechanical properties and better constitutive model were used.

Review and conclusions
The paper presented two validation benchmarks to check the capabilities of Thebes code to simulate the coupled gas-liquid flow in porous media. The numerical results prove that the code replicates well the reported experimental data and numerical solutions by other research groups. The benchmarks, however, concentrated on sandy soils and the gas pressure that develops stayed at relatively low values when compared to the values that might develop in non-isothermal cases such as the case of bentonite used for sealing nuclear waste. Furthermore, only linear elasticity is employed here. Including elasto-plasticity will pose extra problems that need to be carefully addressed and validated.
Having this elementary check successfully done, further validations is needed in the case of nonisothermal situations with more advanced mechanical soil model. That would definitely help in figuring out the limitations of the used theoretical and numerical framework and opens the door for possible improvements. This could be in form of developing constitutive models that consider the effect of pore structure evolution under coupled THMC loads, for example.