Numerical investigation on water exchange of shale samples

Interest in the hydraulic and mechanical characterization of shales has grown in recent years, because of their application in the context of energy geotechnics. In the frame of nuclear waste disposal shales are considered as host formations for the placements of nuclear waste at high depths. In the frame of hydrocarbon production they are considered as unconventional reservoirs, from which extracting natural gas. Understanding how fluids flow through shales is then a key aspect for both fields of application. This paper focuses on the analysis of the transport of water vapour through laboratory samples. After reviewing the balance and flow laws that govern the transport of fluid in unsaturated porous media, a simplified model is put forward. The model was implemented in a commercial finite element code, and it was used to reproduce the results of a literature study on wetting and drying of Opalinus Clay shale samples, imposed through the vapour equilibrium technique. Back analysis of the water content and volume strains of these specimens suggests that existing models underestimate the actual flow rate of water vapour which takes place at low suctions. The current interpretation also seems to be consistent with microstructural investigations on the interconnection between large pores of this material.


Introduction
In the last years the underground geological storage is considered the best solution for the disposal of nuclear waste produced by power plants of different countries including Switzerland, Belgium and France [1]. Repository systems are placed in a stable rock formation, to protect the environment and the people from the radioactive waste that they contain. Shales are studied as host rock formation, because of their sealing properties.
A detailed knowledge of transport processes and mechanisms in shales is thus a fundamental requirement for a safe disposal. Careful attention must be payed to the mass transport phenomena occurring in the partially saturated conditions at which these materials are found.
Experimental studies aim to know how the geomechanical properties of Opalinus Clay shale change with the water content, in light of engineering aspects such as the excavation and the ventilation of the tunnel or the temperature decay of the canister. Much attention has been paid to the determination of water retention behaviour of Opalinus Clay shale, e.g. [2][3].
This study aims at gaining an insight on the transport of water vapour through unsaturated Opalinus Clay shale. It builds on available experimental results of vapour equilibrium tests. Suction equalisation, and related transient processes, are modelled numerically accounting both for the effects of suction on water transport and on volume change [4]. The evolution of the pathways available for gas transport is investigated by the definition at the macroscopic scale of a tortuosity factor, whose dependency on liquid degree of saturation is expressed either with a theoretical expression from the literature or through a dedicated back-analysis.

Governing equations for water transport in unsaturated porous media
Porous geomaterials, as shales, can absorb moisture in vapour or in liquid form from the environment and transport it into their matrix, through different mechanisms. The fundamental equations governing the flow and liquid and gaseous water in porous media are given below.

Mass balance
The differential equations for the mass balance of the fluid were written considering water and dry air as chemical species [5]. Water molecules are present both in the liquid and in the gas phase (as water vapour) and the mass balance accounts for content changes and transport terms both in the liquid and in the gas phase: where n is porosity, S l and S g are the liquid and gas degree of saturation respectively, c l w and c g w are the mass concentration of water in the liquid and gas phase, v l and v g are the velocities of the water and of the gas phase, J hl w and J hg w are the hydrodynamics fluxes of the water species in the liquid and gas phase, Γ l w and Γ g w are water sources in the two phases. The latter terms were neglected as chemical reactions are thought not to take place.
Analogously, the mass balance for the dry air can be written as: with c l a and c g a mass concentration of air in the liquid and in the gas phase; J hl a and J hg a the hydrodynamics fluxes of dry air in the liquid and gas phase; Γ l a and Γ g a sources of dry air in the liquid and gas phase. The porous space is filled with gas and one single liquid phase, thus the saturation rule applies: It is here assumed that the concentration of air into the liquid phase c l a is very small and that it is negligible compared both to the concentration of air in the gas phase c g a and to the concentration of water in the liquid phase c l w , so that c l w = ρ l , with ρ l mass density of the liquid phase. The simulated process refers to the equilibration of samples which are inserted in closed jars, where they are exposed to humid air, and no significant gradients of air pressure are assumed. These considerations allow neglecting the mass balance equation for the air species, introducing equation (3) into equation (1) neglecting v g and J hl w one obtains:

Water content change in time
The derivative of the water mass in time (left-hand side of equation (4)) can be expanded introducing the dependency on total suction ψ of the degree of saturation of the liquid, of porosity and water vapour concentration.
The ratio between the current water vapour concentration and the saturated one c g w | sat is the relative humidity, h r . This is related to total suction through the psychrometric law where M w is the molecular mass of pure water, R is the universal gas constant and T is the absolute temperature. c g w | sat is related to saturated vapour pressure of water p v,sat by the ideal gas law: where according with [6], adopting MPa and Kelvin as units: By introducing equation (5) into equation (4) it follows: is the specific moisture capacity and D = J H is the slope of the change in porosity due to a unit change in total suction.

Flow laws
The flow equation describing the transport of liquid water is Darcy's law, which for an unsaturated material reads: where K is the saturated hydraulic conductivity of the material, k r is the relative permeability, g is gravitational acceleration, u w is water pressure and z is the coordinate in the vertical direction. Assuming constant air pressure and negligible gradients of osmotic suction, allows then writing: The flow of water vapour in a free gas phase can be written as Modified Fick's law [5]: where R 8 is the molecular diffusivity of water vapour in the free gas phase. Accounting for the tortuosity of the porous medium, it follows: where τ τ τ τ is a tensor which accounts for the effects of the tortuosity of the porous network on the transport of the water vapour and R 8 is the diffusion coefficient of vapour in air.
Introducing equation (5) in equation (12) and summing the vapour diffusive flux to the darcian flow of liquid water where the height term was neglected, it follows that: where U * is the effective diffusivity of vapour through soil with Y = 1 − cross-sectional fraction of the soil area available for vapour flow.

Partial differential equation for water flow
The PDE that governs the water flow modelled in this study in terms of total suction, in 3D dimensions is:

Analysis of the vapour diffusion process
The transport of water vapour through the Opalinus Clay shale was studied on basis of the differential equation (15). In particular, the analysis was aimed at investigating the role of the tortuosity, under the assumption of isotropic behaviour, i.e.: where I is the identity tensor. In the frame of equation (16) and accordingly with [5], the coefficient τ has a geometrical meaning, as it is the ratio between the macroscopic straight-line distance of the diffusion path and the effective length of the diffusion path through the soil. The evolution of the τ parameter along drying and wetting steps was determined through back analysis of experimental results of Vapour Equilibrium Technique tests (VET). The obtained values were compared to the predictions of a literature model. The study is based on the experimental data available in [2].

Basic geotechnical properties of Opalinus Clay shale
The Shaly Opalinus Clay shale studied was extracted from the Mont Terri Underground Laboratory in northern Switzerland at a depth of about 300m. The specific gravity of the samples was Gs = 2.75, their initial degree of liquid saturation and porosity were S l = 0.37 and n = 0.19

Experimental methodology for the tested shale
In order to characterise the hydro-mechanical behaviour of the Opalinus Clay shale in [2], the coupling between the water retention and deformative response was studied performing a wetting/drying cycle with the Vapour Equilibrium Technique. Saturated salt solutions were used to imposed total suction ψ = 139, 86, 39, 22.6 and 9.8 MPa. For each suction, two cylindrical specimens were used (diameter 25 mm, height 20 mm). The initial total suction was measured with the WP4C dewpoint potentiometer and it was found to be ψ 0 = 96 MPa. The specimens were placed inside a desiccator under isothermal (25°C) and free stress conditions. One of them was equipped with strain gauges along the axial and circumferential directions to assess the volumetric behaviour. Another one was weighted regularly to monitor the mass variation, with the aim of determining the changes in water content. At the end of the test the specimens were oven dried at 105°C to measure their dry mass. For each imposed suction, the specimen weight initially changed rapidly, while tending to become stable at larger times, which was assumed to indicate equalization of specimen suction with the imposed one. Afterwards, suction was changed again for the following drying or wetting step.

Water retention and volumetric behaviour along drying and wetting
The values of the liquid degree of saturation S l and porosity n obtained at the end of each equalization stage were used to characterise the water retention and the mechanical behaviour as required by the equation (15).
Limitedly for the hydro-mechanical analysis of this study, the water retention was interpreted with a semilogarithm relationship. Although the observed hysteresis is very small, two different fitting relationships were determined, one for the initial wetting path and another one for the following drying path (figure 1).
As observed in [2][3][7-8] the volumetric response of Opalinus Clay shale is dependent on suction changes, mostly because of the high clay content of this material [4]. For modelling purposes, also the relationship between porosity and suction was interpreted adopting a semi-logarithmic relationship, shown in figure 2.

Numerical simulations
The insight on the tortuosity coefficient τ relied on backanalyses performed solving equation (15) with the Finite Element software COMSOL (General Form PDE (g), Mathematics module).

Fig. 2. Volumetric behaviour of Shaly Opalinus
Clay under the effects of total suction changes. Experimental data from [2] and fit with a semi-logarithmic relationship.

Model geometry, boundary and initial conditions
Due to the cylindrical shapes of the specimens, a 2D axisymmetric geometry was used in the simulations. According with experiments in [2], the model specimen had a radius of 12.5 mm and a height of 20 mm.
The complete sequence of the suction equalization steps imposed along the drying and the wetting path was modelled. The evolution in time and space of total suction during each stage was obtained by solving the Partial Differential Equation Furthermore, since no flow of liquid water occurs in free air, the transport of water mass crossing the specimen border was imposed to be due only to water vapour flux. In the simulations, this condition was imposed by assigning null hydraulic conductivity on the superficial layers of elements, which were also kept very thin in order not to impose artefacts on the transport of water within the sample, where the liquid can actually flow.

Model parameters
The following values were used for the model constants: =1000 kg/m 3 , c =0.018 kg/mol; R 8 =8.4·10 -5 m 2 /s and T=298.15 K. The hydraulic conductivity determined in [8], K=10 -14 m/s was used, while the relative permeability was assumed to follow the relationship Z % = ' I + (I + Md! (I + Md! / e (consistently with [5]), with %f0 = 0.2 residual liquid degree of saturation. The fitting relationships in figure 1 and figure 2 were used to set the dependency of S l and n on suction, both for the wetting and drying paths.
A first set of numerical simulations was run adopting the Millington and Quirk relationship [9], which states that τ is a function of both degree of gas saturation and porosity: Equation (17) is widely used in the literature and it is also implemented as a default suggestion in the software COMSOL.
Another set of simulations was run with the aim of determining the parameter τ, and its evolution along the drying and wetting process, through back analysis of the experimental data. First, each drying or wetting step was back-analysed, determining a constant value of τ which allowed the best reproduction of the pertinent experimental data. All the values of τ that were determined this way were then referred to the average degree of liquid saturation of the corresponding experimental step, as shown in figure 3. The numerical results suggest that the value of the tortuosity factor might depend on the direction of the hydraulic path. However, as a first analysis it was decided to neglect hysteretic effects and the following relationship was adopted : where i and k are non-dimensional parameters. The value of the parameters adopted in this study is i = 0.53 and k = 1.5 suctions (ψ ≥39 MPa), but the agreement between the numerical prediction and the experiments is poorer at lower suctions. In fact, there is an underestimation of the rate of the evolution of the water content and the volumetric strains in the last step of wetting and the first step of drying.
A satisfactory reproduction of the results of all the experimental stages is obtained instead with equation (18), with the exception of the last drying step. The latter is poorly reproduced by both models; however, the inaccuracy mostly consists in not reproducing the residual volume strain of the experiments, which is likely to depend on the mechanical model rather than on τ. Figure 6 provides the evolution of the effective vapour diffusivity term Dv* (equation (14)) with S l , as predicted with both relationship (17)     This difference in values has to be appointed to the trend of the tortuosity factor found through back analysis, that increases slightly with the degree of saturation, in contrast with the Millington and Quirk model.

Concluding remarks
This note presents a numerical investigation on the transport of water in shales at high total suctions. It stems from a previous experimental characterization of the water retention of Opalinus Clay shale, and it mainly focuses on the tortuosity of the water vapour path, modelled through an isotropic coefficient τ.
A wetting-drying cycle, imposed with the vapour equilibrium technique, was simulated adopting first the relationship between τ, porosity and liquid degree of saturation of Millington and Quirk. A good reproduction of the evolution with time of water content and volume strains was obtained at very high suctions, both along the wetting and the drying paths. However, the numerical predictions were not very good at lower suctions, where the Millington and Quirk relationship largely overestimates the equilibration times. Results of a preliminary back analysis suggested that during the tests τ slightly increased with the liquid degree of saturation. Such a trend is opposite to the one of the Millington and Quirk formula and it might be explained by the development of small micro-cracks as the water content increases. Possibly because the volume changes that occurred during the tests were small, in the back analysis it was not necessary to account for any dependency of τ on porosity.
Of τ from back-analysis appear to be consistent with the outcomes of a microstructural study on the spatial distribution of the interconnection between macropores (pore radii r>10nm), based on FIB/SEM nanotomography images (FIB-nt) and presented in [10]. In the Opalinus Clay shale, pores with r>10nm represent about 20% of the total porosity, which makes results in [10] relevant for the gas transport when S l ≅ 80 %, i.e. for the lowest suction step simulated in the present work. Very interestingly, according to the FIB-nt the ratio of the effective length of the diffusion path divided by the macroscopic straight-line distance of diffusion (i.e. 1/τ in equation (16)) is between 2 and 3, which perfectly matches with the value 1/τ = 3 obtained with equation (18). On the contrary, at such a degree of liquid saturation and porosity, the Millington and Quirk model would predict 1/τ = 34, which is definitely out of range. Since the diffusion process depends strongly on the microstructural features of the material [11], the Millington and Quirk model might be inadequate for this kind of material where the interconnected pore-structure evolves with the hydraulic state. Even if the macropores seem to be disconnected, they might be linked by smaller pores, and this connection is a key element for the diffusion transport mechanism [10]. Although the latter aspect is not studied directly in this paper, it might contribute to justify the moderate increase of the tortuosity coefficient and of the vapour diffusivity parameter with suggested by the results of the numerical analyses.