Modeling of thermal pressurization in tight claystone using sequential THM coupling: Benchmarking and validation against in-situ heating experiments in COx claystone

We apply thermoporoelasticity and a sequentially coupling technique for modeling thermally-driven coupled Thermo-Hydro-Mechanical (THM) processes in tight claystone. A THM benchmark case with a corresponding analytic solution for thermoporoelasticity under a constant heat loading verifies the model. Thereafter, two in situ heating experiments are simulated for model validation: a smaller-scale heating experiment (TED experiment) and a larger-scale experiment (ALC experiment) in Callovo-Oxfordian (COx) claystone at the Meuse/ Haute-Marne underground research laboratory in France. The model exhibits good performance to match the observed temperature and pore pressure evolution for the smaller-scale TED experiment. For the larger-scale ALC experiment, general trends of thermal-pressurization are captured in the modeling, but pressure is under-estimated at some monitoring points during cool-down. This indicates that the THM response in the field may be affected by the variability of rock’s properties or irreversible or time-dependent mechanical processes that are not included in the current thermoporoelastic model. The main contributions of this work are as follows: (1) we verify and validate the numerical simulator, TOUGH-FLAC, to be a valuable coupled THM modeling tool; (2) prove that the laboratory determined material parameters can be used as reference values for upscaling ex- periments. However, to better identify and quantify THM processes with modeling of in situ tests, more em-phasize should be dedicated to obtaining high-quality mechanical deformation data.


Introduction
Claystone is currently considered as a potential host material in nuclear waste disposal. Claystone has favorable properties, such as high absorption capability and low permeability, to retard subsurface contaminant transport, and prevent the hazardous pollution from reaching humans or ecological systems. Since the host rock around heat releasing nuclear waste would be subjected to strongly coupled Thermo-Hydro-Mechanical (THM) processes, it is necessary to study how these coupled effects could impact the repository performance over the long term, i.e. up to 100,000 years (Rutqvist et al., 2014). A number of in situ heating experiments to study coupled THM processes in claystone have been conducted at several underground research laboratories over the past decades. These include in situ experiments on Callovo-Oxfordian (COx) claystone in France (Conil et al., 2012;Armand et al., 2017a), Boom clay in Belgium (François et al., 2009;Chen et al., 2011), and Opalinus clay in Switzerland (Gens et al., 2007;Garitte et al., 2017). Thermally-induced changes in pore fluid pressure can be significant in tight claystone due to its low permeability and the difference between the thermal expansion coefficients of the fluid and that of claystone (Muñoz et al., 2009;Ghabezloo and Sulem, 2010;Rutqvist et al., 2014). This increase in the pore pressure would lead to a reduction of the effective stress and loss of strength of materials, which could result in shear failure or hydraulic fracturing (Ghabezloo and Sulem, 2010). Thermal pressurization has also been studied in earthquake science to describe the dynamic rupture propagation of faults during earthquake nucleation (Schmitt et al., 2011;Andrews, 2002).
In this work, we apply coupled THM numerical modeling to study thermal-pressurization phenomena at the Meuse/Haute-Marne (MHM) Underground Research Laboratory (URL), in Bure, France (Armand et al., 2017a). COx claystone is considered as a potential host rock for a nuclear waste repository in France by ANDRA (the French national radioactive waste management agency), and coupled THM behavior of the COx claystone is of great importance for the design and https://doi.org/10.1016/j.tust.2020.103428 Received 10 December 2019; Received in revised form 4 March 2020; Accepted 26 April 2020 performance assessment of the repository (Armand et al., 2017a). Apart from the in situ studies mentioned above, the THM behavior of COx claystone has been extensively studied by laboratory testing. These studied include undrained thermal compression tests to understand the temperature dependency of rock properties (Ghabezloo and Sulem, 2010;Ghabezloo and Sulem, 2009;Mohajerani et al., 2012), and tests to investigate Hydro-Mechanical (HM) parameters under unsaturated conditions (Charlier et al., 2013). Other laboratory studies include those of Zhang et al. (2017) to examine thermal effects on COx claystone, including thermal expansion, thermal-pressurization, temperature influence on mechanical properties, as well as fracture permeability and sealing. Poroelastic parameters of COx claystone have been experimentally investigated by Belmokhtar et al. (2017a,b), including Biot's coefficient, drained and undrained bulk moduli, as well as thermal volumetric and creep behavior. Recently, Braun et al. (2019) proposed tests involving thermal and mechanical loading to measure drained and undrained parameters of COx claystone. Such experimental work and in situ studies at the MHM URL has led to substantial knowledge about the THM behavior of COx claystone and a set of proposed material parameters for coupled processes modeling of the MHM URL heating experiments in COx claystone (Armand et al., 2013(Armand et al., , 2017aConil et al., 2020).
Numerical modeling is necessary for making performance assessment calculations of a nuclear waste repository for time periods far beyond any experiments. The model should be verified against analytical solutions for correctness and should be validated at the field scale by simulating in situ experiments under multiple processes, in this case thermally-driven coupled THM processes in tight claystone. A number of coupled modeling approaches exist, including so-called fully coupled method, in which all coupled equations are solved simultaneously, in contrast to iterative coupled approach, in which equations are solved sequentially (Kim et al., 2011;Rutqvist, 2017). In this work, we apply a modular approach originally proposed by Settari and Mourits (1998) for linking an existing reservoir simulator (fluid flow and heat transport) with a geomehanical simulator. This method has been mostly investigated and applied for modeling coupled HM processes through a carefully derived pore-volume coupling, or porosity correction (Kim et al., 2011). Here we apply this approach based on the theory of thermoporoelasticiy within the framework of the TOUGH-FLAC simulator, with verification against analytic solutions and validation against field experiments.
This work is part of the international DECOVALEX-2019 project, in which modeling of strongly coupled THM processes of COx claystone is a modeling task (Birkholzer et al., 2019). In the following Section 2, we start from the fundamental thermoporoelasticity theory with a general form of the Helmholz free energy, and implement corresponding equations into a coupled simulator using a sequential coupling technique. To assure the correctness of the implementation, we conduct a careful verification against an analytic solution for thermoporoelasticity (Section 3). These includes corrected and extended analytical solutions to consider the full spectrum of temperature, pressure, displacement, stress and effects of Biot's coefficient. Having verified the simulator, we then proceed to model two in situ heating experiments, the TED and ALC experiments, performed in COx claystone at the MHM URL, in Bure, France (Section 4). TED is a small-scale experiment ("propriétés et effects en TEmpérature Deux" in French, "a second experiment that studied the thermal properties of the COx" in English), whereas ALC is a large-scale experiment ("ALvéole Chauffante" in French, "heating cells" in English). Parameter studies are performed to understand the impact of the different parameters on strongly coupled THM responses that are observed in the field tests. We conclude and provide recommendations for further studies to improve the state of knowledge and our abilities of making confident model predictions of coupled THM processes in tight claystone.

Thermoporoelasticity
The fundamental theory used here is based on thermoporoelasticity, which extends thermoelasticity to porous continua by considering an underlying thermoelastic skeleton. The general form of the Helmholtz free energy s for a linear saturated thermoporoelastic material is defined as (van Duijn et al., 2019;Coussy, 2004 where s0 is the free energy at the reference state with free strain; is the strain tensor; is the tensor of skeleton tangent elastic stiffness modulus; is the density of the porous medium; C p is the specific heat capacity; T 0 and T are the absolute temperature at the reference state and current state; N is the Biot's tangent modulus linking the pressure variation and the porosity variation; p is the current pore pressure; p 0 is the pore pressure at the reference state; b is the Biot's tangent tensor; is the tensor of skeleton tangent thermal dilation coefficients; and 3 is the volumetric thermal dilation coefficient related to the porosity.
Within the framework of thermodynamics, constitutive equations of porous media can be derived from the Helmholtz free energy. The stress in solid skeleton is obtained as The porosity of porous medium is conjugated to the pore pressure in thermodynamics: The porosity can also be linked to the fluid part: Combining Eqs. (3) and (4), the fluid continuity equation is where m f is the fluid mass content; f is the fluid density; is the fluid tangent bulk modulus; and = + m f , where 3 is the volumetric thermal dilation coefficient related to the porosity and 3 f is the fluid tangent coefficient of volumetric thermal dilation (Coussy, 2004;Xu and Prévost, 2016).
The fluid flux is calculated by Darcy's law: where F f is the fluid flux vector, k is the intrinsic permeability, µ f is viscosity of the fluid, and g is the gravitational acceleration.
The flow of heat in the soil is assumed to be governed by Fourier's law: Here F h is the heat flux vector, and t is the thermal conductivity.
H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 2 The heat storage term in the porous media is where C s is the specific heat of the solid skeleton and u f is the specific internal energy in the fluid. The heat flux studied here is assumed to include conduction and convection: Here t is the thermal conductivity, and h f is the specific enthalpy in the fluid.
In the solid skeleton, the momentum balance equation is solved according to: where · is the divergence operator, v s is the velocity of the solid skeleton, and t is the time.
If the material is assumed isotropic, the second order tensor b k , , and reduce to scalars, b k , , and . Moreover, the bulk modulus K can be used to calculate the mean stress.

simulator and Coupling Approach
We implement thermoporoelasticity and model coupled THM processes in the framework of TOUGH-FLAC numerical simulator, which links the multiphase fluid flow and heat transport simulator, TOUGH2 (Pruess et al., 2012), with the finite-difference geomechanical code, FLAC3D (Itasca, 2009). The simulator was firstly developed by Rutqvist et al. (2002) and has been used for coupled THM processes modeling of a wide range application requiring multiphase fluid flow and geomechanics, such as nuclear waste disposal, geologic carbon sequestration, geothermal energy and hydrocarbon exploration (Rutqvist, 2011;Rutqvist, 2017). The current version of TOUGH-FLAC includes a sequential coupling scheme corresponding to a fixed stress-split method (Kim et al., 2011;Blanco-Martín et al., 2017), which is here applied to investigate thermally-driven coupled THM in tight claystone. In this method, fluid flow equations are solved first under fixed stress in TOUGH2; then pressure and temperature are passed to FLAC and prescribed during mechanical simulations. A porosity correction c is derived from the constitutive equations of porous media for this scheme implementation (Kim et al., 2012): where K s is the bulk modulus of solid grains, = s is the linear thermal dilation coefficient related to the solid grains. v is the mean total stress, and v is the volumetric strain. The porosity change in Eq. (11) is implemented in TOUGH2 to solve the flow transport at the current time step, while the porosity correction  Volumetric coefficient of thermal expansion of water [1/K] = × 4 10 w 4 * Note: the analytical solution in Booker and Savvidou (1985) assumes the water is incompressible.
H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 in Eq. (12) is computed in FLAC3D from previous time step and then passed to TOUGH2. Eq. (11) is different from that in Kim et al. (2012) with our derivation leading 3 s , instead of b 3 s . The correctness of Eq. (11) is verified through an excellent agreement with the analytical solution as presented in the following Section 3.

Verification of model implementation
In this section, we conduct benchmark simulations against analytical solutions to verify the correctness of the code implementation of the sequential fixed-stress THM coupling algorithm. The basic benchmark test considers thermally-driven THM coupled processes of an infinite homogeneous saturated porous medium around a point heat source. Booker and Savvidou (1985) and Smith and Booker (1993) provide the analytical solution for this problem and all equations are summarized in A, including some corrections of the original equations. Moreover, in the analytical solution presented in A, we retain the option to vary the value of Biot's coefficient, which is an important parameter in this study of COx claystone. The analytical solution is based on the hypothesis that the pore water and the solid grains are incompressible, so the term p d M 1 in Eq. (5) is zero. As a result of that,  Booker and Savvidou (1985).
H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 the variation of fluid mass is solely due to mechanical deformations and thermal expansion of the fluid and the porous media.

Verification against analytical solution
A 3D TOUGH-FLAC model is generated as a cube, 30 m × 30 m × 30 m, as shown in Fig. 1. Considering three symmetry planes, only 1/8 of the domain is simulated. The initial temperature is set to 0°C, while the pore pressure and stresses are both set to 0 Pa. Regarding thermal and hydraulic conditions and symmetry conditions, the three symmetry planes are impermeable and adiabatic. At the outer boundaries, the temperature and pore pressure are set to 0°C and 0 Pa, respectively. At the heat source, a constant heat power of = Q 700 W (700/8 = 87.5 W for the 1/8 symmetric model) is instantaneously applied at = t 0. For mechanical conditions, all boundaries are free to move except the symmetry planes where zero displacement conditions are applied normal to the boundaries.
A homogeneous and isotropic material is considered in this benchmark. The model parameters including two components, the solid phase and pore water, are listed in Table 1. These parameters were defined based on site investigations of the COx claystone at MHM laboratory, though they are simplified as isotropic. Comparison between modeling results and analytical solutions for temperature, pore pressure, displacement and total stress evolution up to 100,000 h (about 11 years) is provided at monitoring points as listed in Table 2.
The simulation results are in excellent agreement with the analytical solution for temperature, fluid pressure, displacement and total stress as shown in Fig. 2. As Fig. 2a displays, the calculated temperature at four H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 5 positions match the analytical solutions accurately. Fig. 2b shows that the thermal-pressurization at all points along with the temperature rise is captured and matches the solutions accurately. After reaching the peak pressure, the fluid pressure slowly dissipates to zero. Fig. 2c presents the displacement evolution at point P4, indicating that the claystone expands in all three directions by the combined effect of temperature and pressure increases based on the simulation of thermoporoelasticity. After about 200 h, the displacements begin to decline, as a result of decreasing pressure. The displacements do, however, not completely rebound to zero as the temperature still remains elevated to the end. An excellent agreement with the analytical solution is achieved for normal and shear stresses, although there are slight differences for normal stresses towards the end (Fig. 2(d) and (e)). These deviations are caused from the fact that the infinite domain is assumed in the analytical solution, while the current model's domain is only 30 m × 30 m × 30 m, and the free stress boundary is imposed at the outer boundaries. These deviations can be reduced if the domain size increases.
The overall excellent agreement with the analytical solution verifies the implementation of the sequential coupling scheme for strongly  H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 coupled thermally-driven THM processes and properties corresponding to tight COx claystone. This includes the porosity correction made through Eqs. (11) and (12), and therefore verifies the corrected term made for Eq. (9). Here we have to consider the sequential coupling scheme in which TOUGH2 calculates porosity change based on the current pore pressure and temperature (Eq. (11)), while the porosity correction is computed from previous stress state (Eq. (12)), which is one step behind. Consequently, highest accuracy can be assured by choosing sufficiently small step size.

Biot's coefficient effects
Biot developed the governing equations to couple three dimensional fluid flow with mechanical deformation for linear elastic porous media (Biot, 1941). Biot's coefficient is a key parameter in poroelasticity theory to couple the flow transport and mechanical analysis, and strongly depends on the microstructure of the porous medium (Tan and Konietzky, 2014;Lion et al., 2005;Salimzadeh et al., 2018). We can easily extend the work by Booker and Savvidou (1985) to account for Biot's coefficient, and verify it with the work of Smith and Booker (1993), which presents a general form of analytical solutions for coupled THM problems. The newly developed solutions are summarized in A as well. Here we keep all other parameters the same as listed in Table 1, but only change Biot's coefficient. Fig. 3 displays the results with = b 0.6, 0.8, and 1.0 with excellent agreement between numerical simulation results and the analytical solution. A reduction of Biot's coefficient corresponds to that the rock grains becomes softer, while the coupling between mechanical deformation and fluid flow becomes weaker. With a reduction in Biot's coefficient, the magnitude of deformation and stress in the porous medium decreases, while the porosity increases less and thereby induces higher pore pressure. This benchmark case verifies that the current model can capture the effect of Biot's coefficient in coupled THM processes, which will enable us to confidently simulate and study the impact of Biot's coefficient on the observed responses at the MHM URL in situ experiments.

Validation against in-situ experiments
With the coupled simulator verified in Section 3 for predicting thermally-driven coupled THM responses in tight claystone, the next step is to validate the model by simulating the two in situ heating experiments in COx claystone at the MHM URL. The two in situ experiments are (Armand et al., 2017b): • The TED experiment: a borehole heating experiment focused on the THM behavior of the undisturbed claystone; • The ALC experiment: a heating experiment focused on the THM behavior of the claystone surrounding a nuclear waste emplacement micro-tunnel.
Both in situ experiments involve heating of the COx claystone to investigate thermal-pressurization and geomechanical responses sujected to coupled THM processes. The coupled modeling of these experiments are conducted in these steps: • A model prediction of the TED borehole heating experiment using H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 material properties of the COx, as recommended best parameters based on previous laboratory tests on COx claystone.
• Comparing the predicted results with that of measurements at TED in situ test and making necessary adjustments (calibration) of material parameters to better match modeling and measurements.
• Using the TED validated model to perform a blind-prediction (field data not known to the modelers) of the coupled THM responses at the larger-scale ALC experiment.
• Comparing the predicted results with that of measurements at the ALC and investigating any discrepancies that might exist between numerical results and experimental data.
Using this stepwise modeling approach, we address some fundamental questions related to repository development in claystone, such as how representative laboratory determined properties are for predicting in situ coupled processes, and how useful smaller-scale heating tests are for calibrating material properties at a larger repository emplacement micro-tunnel scale. The model predictions are evaluated by comparison to measured temperature and pressure responses, whereas the mechanical deformation measurements were not of sufficient quality to make useful comparisons to simulated results. The modeling of these two experiments are presented in the following Section 4.1 for the TED experiment and 4.2 for the ALC experiment.

TED experiment
The TED experiment was conducted over 3 years between 2010 and 2013 and involved three heaters in three parallel boreholes separated of about 2.7 m (Fig. 4). The three heaters, 4 m long, were installed at the end of 160 mm diameter and 16 m long boreholes, drilled from a tunnel (the GED tunnel) and parallel to the maximum horizontal stress. This arrangement represents a similar configuration to high-level nuclear waste cells in the French waste disposal concept, with emplacement in parallel micro-tunnels, but at a smaller scale. The TED experiment was heavily instrumented with 108 temperature sensors in the rock mass, 69 temperature sensors in the 3 heater boreholes, 18 piezometers, 2 extensometers and inclinometers, and 10 temperature sensors recording the temperature at the level of the main drift. The temperature measurements recorded during the TED experiment show that the COx claystone has an anisotropic thermal conductivity; at the same distance from the heater, the temperature increase is higher along the bedding plane direction than in the perpendicular direction. Observations of pore pressure also show that its evolution depends on the location with respect to the bedding; following a thermal power increase, the pore pressure increases faster in the direction parallel to bedding than in the perpendicular direction.
We model the THM response of COx claystone in the TED experiment with material parameters based on previous laboratory tests and field observations at the MHM laboratory (Armand et al., 2017b;Conil et al., 2020). We then conduct an initial prediction of the TED field test and compare against field monitoring data. The TED field test is discretized into a cube with a side length of 50 m centered in height at = z 0, i.e. exactly as the GED tunnel center at a depth of 490 m (Fig. 5). The three heaters are embedded at the center of the domain and are surrounded by refined grids. It is assumed that the whole domain remains saturated during the experiment. The boundary and initial conditions are summarized in Table 3. The drilling of the heater boreholes and the extensometer boreholes are explicitly modeled with fixed atmospheric pressure leading to water flow from the host rock into these boreholes. For the other boreholes, they were drilled to place the sensors, then were backfilled to ensure sealing and low compressibility, and to reproduce the low permeability of the COx. During the experiment, only the heater boreholes and the measurement boreholes  H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 (TED1230 and TED1231) showed to be draining. Thus, their influence on the results is minimized and treated as COx in the numerical model. Besides that, we simulate the excavation of the GED tunnel, the drilling of monitoring boreholes, and account for an Excavation Damaged Zone (EDZ) around the GED tunnel with increased permeability. The simulation starts with instant excavation of the GED tunnel by applying a cylindrical boundary representing the surface of tunnel, followed by the drilling of other boreholes by deactivating the elements of boreholes in the mechanical simulator. Then, 506 days after the excavation, the heating phase starts and is running for about 1251 days. The time zero, corresponding to the excavation of the GED tunnel, is at April 21, 2008. Heater 1 (TED1201) was turned on first, using three steps to reach the planned heat power, 600 W. After operating heater 1 for 400 days, heaters 2 (TED1202) and 3 (TED1203) were turned-on and increased to the power of 600 W in three steps. The raw data of heat power shown in Fig. 6 were input into the model. A cylindrical 1-m thick EDZ is assumed around the GED tunnel with enhanced permeability ( = × k 1 10 18 m 2 in Fig. 5(a)). During the heating stage, the measured temperature at the GED tunnel plotted as OHZ1290 and TED1270 in Fig. 7 is smoothed to a pseudo-sinusoidal function and applied as the temperature boundary conditions on the GED tunnel wall.
Water properties in the TOUGH-FLAC simulation are calculated from the steam table equations (International Formulation Committee, 1967), which is standard in TOUGH2. As mentioned, we first conducted a model prediction using recommended properties from previous Fig. 11. The ALC experiment at Bure with various monitor boreholes and micro-tunnel ALC1604 that are used by modeling teams in DEOVALEX-2019, Task E, for predictive and interpretative modeling (Armand et al., 2017b). H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 laboratory and sites investigations of COx claystone. These included anisotropic thermal conductivity, = 1.96 W/m/K, = 1.26 W/m/K, and permeability, = × k 6 10 20 m 2 , = × k 3 10 20 m 2 . However, with these parameters, we could not match the in situ experimental observations exactly with our model. Thus, we recalibrated the parameters to achieve a better match between simulated and measured temperature and pressure data. The in situ calibrated thermal conductivity values are = 2.05 W/m/K, = 1.15 W/m/K, and calibrated permeability values are = × k 3 10 20 m 2 , = × k 0.7 10 20 m 2 , as listed in Table 4. The in situ calibrated values for thermal conductivity and permeability are about the same magnitude as the laboratory determined, but the degree of anisotropy is somewhat more prominent. All other parameters utilized in the current model are the same as the laboratory experiment measurements and summarized in Table 4 as well.
In the TED experiment, 6 sensors are placed at different boreholes near heaters to measure the temperature evolution during the heating stage. Fig. 9 displays the simulated temperature results and their comparison with experimental data. A good agreement between model simulations and observed temperature is achieved at sensors 1210_05, 1219_05, 1250_01, and 1251_01, which are located 0.61 m, 0.57 m, 0.32 m, and 0.42 m respectively from the heater 1-TED1201 (Fig. 8), while at sensors 1253_01 and 1258_01, located farther away (1.14 m and 1.36 from the heater 1 as shown in Fig. 8), the model overestimates the temperature by 3-5 degree. However, the data measured by these two sensors are reported as unreliable and therefore the overall result is that the temperature evolution can be accurately modeled using TOUGH-FLAC with the parameters listed in Table 4. Fig. 10(a) represents the simulation results and experimental observations of pore pressure at the positions of sensors 1253_01 and 1258_01, which are located on the side of and below the heater 1, respectively, as displayed in Fig. 8. Since the horizontal permeability is 4 times higher than the vertical one, and sensor 1253_01 is closer to the heater and extensometer TED1230, the pore pressure at 1253_01 decreases more than at 1258_01. From the figure, our simulated pore pressure at boreholes 1253 and 1258 is in good agreement with experimental data, although some disagreements of pore pressure at sensor 1258_01 are noted during the beginning of the heating phase. The early pressure reduction is due to the excavation process which is well captured in the modeling for 1253, but not so for 1258. Fig. 10(b) displays the simulated and observed pore pressure at 5 sensor positions in borehole 1240. These 5 sensors are placed on the side of the heating zone and parallel to the heater boreholes as plotted in Fig. 8. Compared with experimental data, the general trend of pore pressure change has been captured in the model simulation, though the peak pressure towards the end of experiment tends to be underestimated.
To conclude, the model simulated with calibrated parameters captures accurately the temperature evolution as well as the general trends of pore pressure as observed in the field. In the simulation, we also applied sophisticated conditions, such as excavation of tunnel and boreholes, and EDZ around the tunnel and drainage boundary at boreholes. The temperature and pressure changes also induce significant deformation and stress changes. However, no reliable measurements of mechanical deformations and stress are available from the TED experiment for comparison to the simulation results (Conil et al., 2020).

ALC experiment
The ALC experiment started its heating phase in 2013 and is an ongoing in situ heating test performed in the MHM URL (Armand et al., 2017b). The experiment is a full scale representation of a single highlevel waste cell in COx claystone according to ANDRA's 2009 reference concept of high-level radioactive nuclear waste disposal (Fig. 11). The ALC1604 micro-tunnel was drilled 25 m into the host rock from the GAN drift, and is oriented parallel to the major horizontal stress, H . The heated part in the ALC experiment is located in the ALC1604 microtunnel between 10 and 25 m along its length (Fig. 11) and is made up of five heating elements. Each heating element is 3 m long and has a diameter of 508 mm. The ALC experiment is instrumented with temperature sensors, relative humidity sensors, piezometers, strain gauges, and displacement sensors. The temperature measurements recorded during the ALC experiment show a similar phenomenon as in the TED experiment, indicating an anisotropic thermal conductivity. Observations of pore pressure show that its evolution depends on the location with respect to the bedding, and strong HM coupling induces opposite pore pressure change near ALC1604 after its excavation. In the vertical direction, the volumetric strain is positive (volumetric expansion), induing the pore pressure decrease, while in the horizontal direction, the volumetric strain is negative (the volume decreases), driving pore pressure increase.

Simulation results of the ALC experiment
The geometry of the model domain is a cube with a side length of 50 m centered in height at = z 0, with a half of the GAN tunnel excavated along the y-direction, and a half of the GRD tunnel excavated H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 along the x-direction (Fig. 12(c)). The five heating elements are located in the ALC1604 between 10 and 25 m along x-direction. These heaters are discretized with refined elements. The ALC1604 micro-tunnel is supported with a steel casing which is explicitly modeled, including a gap between the casing and the host rock (Fig. 13). An additional draining borehole, ALC4004, is explicitly simulated with excavation and drainage as shown in Fig. 12(c). The ALC experiment is conducted close to the location of the TED experiment, so we start an initial blindprediction on the THM response of COx claystone with the calibrated material parameters from the TED experiment (Table 4).
In the simulation, the model starts with instant excavation of GAN and GRD tunnels, followed by drilling of other boreholes. Then at 458 days after the excavation, the heating phase starts and is running for about 1500 days. The time zero corresponds to the time for the H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 excavation of tunnels, which occurred on November 1, 2011. A heating test at a very low power (33 W/m, 495 W in total) was conducted between January 31 and February 15, 2013. The main heating phase started on April 18, 2013, at a constant nominal power of 220 W/m (3300 W in total) for the 15 m occupied by the heater elements (i.e., 660 W per element). The initial pore water pressure is considered uniform at 4.7 MPa in the entire domain when the excavation starts. The initial stress is anisotropic with same values as applied in the TED experiment. That is, the vertical stress is set to 12.7 MPa being the intermediate principal stress, whereas the major horizontal stress H is 16.1 MPa (y-direction) and the horizontal minor stress h is 12.4 MPa (x-direction). The gap between casing and rock is assigned air-like properties after the heaters are placed inside the micro-tunnel. Most importantly, no mechanical resistance is applied, and the thermal conductivity of air is assumed. Around 400 days after the excavation, the field measurements indicated that the gap was sealed. Thus, in the model, the gap zone is replaced with claystone properties when this closure occurs.
In the field, 6 sensors are placed at different boreholes near heaters as displayed in Fig. 14 to measure the temperature evolution during the heating stage. Sensors 1617_01, and 1617_02 are located above the ALC1604 micro-tunnel. Sensors 1616_03 and 1616_05 are placed on the side of the ALC1604, but 1616_03 is near the heating elements, while 1616_05 is close to the GAN tunnel. Sensors 4005_02 and 4005_04 are installed at borehole 4005, which is on the other side of the ALC1604 micro-tunnel and above the ALC4004 extensometer borehole. Fig. 15 displays the temperature results compared with experimental data recorded by the sensors. As the figure shows, the model prediction matches the experimental data well, though some small discrepancies occur at some sensors located farther away from the heat source.
We compare numerical results and observations of pore pressure in Fig. 16. The pore pressure decreases after the initial excavation of the GAN and GRD tunnels, which is due to the drainage conditions applied on tunnel surfaces. Then, around 0 day, corresponding to October 23, 2012, when the ALC1604 micro-tunnel is drilled, HM responses due to excavation occur. Volumetric expansion due to excavation of ALC1604 occurs at sensors 1617_01 and 1617_02, which are above ALC1604. As a result, pore pressure at these two positions drops. On the side of the tunnel, volumetric compression occurs at sensors 1616_02 and 1616_05, which raises the pore pressure. After the heating starts, the pore pressure increases due to thermal-pressurization in all sensors except at sensor 1616_05, which is close to the GAN tunnel, and affected by the atmospheric pressure boundary on the surface of the tunnel. At sensor 4005_04, the pore pressure also grows with temperature, but since it is far from the heaters, the temperature changes are small and the peak of pore pressure is relatively low. The effect of gap closure is negligible on the temperature evolution curve based on the analysis. Fig. 17 presents a plane view of the pore pressure evolution below the plane = z 0. At the beginning of the computation, i.e. before time = t 0, the pore pressure is almost uniformly distributed, except near GAN and GRD tunnels, where the temperature fluctuation in the tunnels has a slight effect on the pore pressure. After the excavation of the ALC1604 micro-tunnel and ALC4004 extensometer borehole, the pore pressure in these two parts changes to the atmospheric pressure for drainage. In addition, the excavation of ALC1604 induces volumetric compression and local pressure build-up on each side of this microtunnel. At 258 days, which is 81 days after the main heating phase started, the water pressure starts to increase due to the thermally-induced pressurization. Tracking the results in Fig. 17, from 258 to 1402 days, it can be observed that the pore pressure near the heated ALC1604 micro-tunnel grows and propagates farther into the host rock. The magnitude of pore pressure rises sharply during the first 200-500 days of the heating, then starts to dissipate back towards the hydrostatic state.
Overall, with the calibrated parameters from the smaller-scale TED in situ experiment, the results showed excellent model prediction of temperature, and the modeling could capture the pore pressure change due to volumetric expansion or contraction at different locations in the field induced by the excavation of micro-tunnel. However, the thermally-induced pressure increases during heating phase were underestimated, especially for the predictions at sensors 1617_01 and 1617_02 that are located relative far above the ALC1604 micro-tunnel.

Parametric studies on Biot's coefficient and permeability
In this section, we present parametric studies on Biot's coefficient and permeability of COx claystone to investigate the corresponding effects of these variables and try to improve the model prediction against the experimental observations, especially the pore pressure at sensors 1617_01 and 1617_02. Three Biot's coefficients ( = b 0.6, 0.7, and 1.0) are assigned to the COx claystone in the parametric study. All other parameters are kept the same as listed in Table 4. Based on the results presented from Fig. 18(a) to (f), in general, with a decrease of Biot's coefficient, the magnitude of peak pore pressure is reduced. However, the change of Biot's coefficient does not affect the trend of the pore pressure evolution. Thus, only changing the values of the Biot's  coefficient does not improve the prediction of pore pressure. Fig. 19(a)-(f) illustrate the pore pressure evolution with the parametric study on permeability. In general, a reduction of permeability results in a more effective thermal-pressurization, inducing a higher peak pressure occurring at a slightly later time. Through this study, we can expect that lower permeability will help to capture the pore pressure evolution trend at sensors 1617_01 and 1617_02. However, it will worsen the prediction of pore pressure at sensor 1616_02. Besides that, too small permeability is not realistic for COx claystone. Thus, it seems impossible to achieve a good agreement at all monitoring locations considering homogenous and elastic porous media in the coupled THM model. This may indicate effects of heterogeneous permeability or some irreversible or time-dependent processes. Based on the measurements, the permeability's coefficient of variation (CV, the ratio of the standard deviation to the mean value) is about 0.98, higher than other properties, such as 0.19 for the thermal conductivity or 0.07 for the specific heat capacity (Plúa et al., 2020). Thus, heterogeneity of pearmeability has more impacts over the space, but it is not included in the current linear thermoporoeastic model.

Concluding remarks
We applied thermoporoelasticity theory and a numerical model to investigate the thermally-induced coupled THM processes in claystone. The fundamental theory extends thermoelasticity to porous continua, and a general form of the Helmholtz free energy for a linear saturated thermoporoelastic material is used to derive the conjugated state variables, such as the stress, and porosity. We implemented the equations into the TOUGH-FLAC simulator using a sequential coupling scheme, known as the fixed stress-split method to couple the multiphase fluid flow reservoir simulator with geomechanical analysis. Some conclusions can be drawn based on the simulation results: • A THM benchmark with an analytic solution of thermally-induced coupled THM processes under a constant heat load verified the correctness of the model implementation; • A borehole in situ heating experiments, the TED experiment, was simulated, for validation and in situ calibration, showing that a good agreement could be achieved in temperature and fluid pressure evolutions with parameters derived from small-scale laboratory experiments. In-situ calibrations indicated somewhat higher degree of anisotropy of thermal conductivity and permeability.
• A larger-scale in situ heating experiment of a simulated nuclear waste emplacement micro-tunnel, the ALC experiment, was simulated as a blind prediction using the same material parameters as in the previous TED experiment. The results showed excellent model prediction of temperature, while the thermally-induced pressure increases were somewhat underestimated, especially after 1000 days.
A parametric study was performed to investigate the cause and potential for improving the prediction of the pressure evolution during cool-down, including changes in Biot's coefficient and permeability. However, it is impossible to achieve good agreement at all monitoring locations considering homogenous and elastic porous media in the coupled THM model. This may indicate some irreversible or time-dependent processes that are not included in the current linear thermoporoeastic model. In order to resolve some of these questions, there is a need for additional higher-quality mechanical data for detailed comparison of simulated and observed mechanical responses. New fiberoptic monitoring might enable such higher quality mechanical

Acknowledgments
DECOVALEX is an international research project comprising participants from industry, government and academia, focusing on development of understanding, models and codes in complex coupled problems in sub-surface geological and engineering applications; DECOVALEX-2019 is the current phase of the project. The authors appreciate and thank the DECOVALEX-2019 Funding Organizations Andra, BGR/UFZ, CNSC, US DOE, ENSI, JAEA, IRSN, KAERI, NWMO, RWM, SÚRAO, SSM and Taipower for their financial and technical support of the work described in this paper. Funding for LBNL's H. Xu, et al. Tunnelling and Underground Space Technology 103 (2020) 103428 modeling work was provided by the Spent Fuel and Waste Science and Technology, Office of Nuclear Energy, of the U.S. Department of Energy under Contract Number DE-AC02-05CH11231 with Lawrence Berkeley National Laboratory. The statements made in the paper are, however, solely those of the authors and do not necessarily reflect those of the Funding Organizations. (A.8) Note that in Eqs. (A.7) and (A.8), we corrected two typos in these equations presented in Booker and Savvidou (1985). The term on the right side of the equations should be g f (3 ) instead of f g (3 ). In these equations, q is the heat input density, Q is the point heat power, and is the thermal diffusivity of the porous medium given by c v is the coefficient of consolidation (Coussy, 2004). In general, it can be expressed as where K H is the hydraulic conductivity, L and G L are the Lame constants of the soil skeleton, and w is the specific weight of water.
u is equivalent volumetric thermal expansion coefficients of porous medium, can be computed as in which f and s are the volumetric thermal expansion coefficients of pore water and solid grains, respectively. R is the distance to the heat source: The functions f g f , , and g are defined such as: with X Y , , and Z such as In the work by Booker and Savvidou (1985), they assume the Biot's coefficient equals to 1, and the fluid is incompressible. Thus, In this way, Eqs. (A.10) and (A.18) can be simplified as presented in the reference (Booker and Savvidou, 1985).