Characterization of Boom Clay anisotropic THM behaviour based on two heating tests at different scales in the HADES URL
-
Published:August 24, 2023
- Standard View
- Open the PDF for in another window
-
CiteCitation
Guangjing Chen, Xiangling Li, Arnaud Dizier, Jan Verstricht, Xavier Sillen, Severine Levasseur, 2023. "Characterization of Boom Clay anisotropic THM behaviour based on two heating tests at different scales in the HADES URL", Geological Disposal of Radioactive Waste in Deep Clay Formations: 40 Years of RD&D in the Belgian URL HADES, X. L. Li, M. Van Geet, C. Bruggeman, M. De Craen
Download citation file:
- Share
Abstract
To examine the impact of the heat generated by high-level radioactive waste on Boom Clay, two heater tests have been launched in the HADES underground research facility: the small-scale ATLAS Heater Test and the large-scale PRACLAY Heater Test. The major objective of these tests is to confirm and refine the thermo-hydro-mechanical (THM) constitutive models and associated parameter values obtained from a laboratory characterization programme. This paper presents the observations from the ATLAS and PRACLAY heater tests and the combined numerical modelling of these tests. To characterize the excavation damaged zone in the clay around these tests, a mechanical model with a strain-dependent elastic modulus is introduced for the Boom Clay. The consistency between the observations from laboratory tests and in-situ tests and the outcomes from the numerical models strengthen the confidence in our understanding of the THM behaviour of Boom Clay. They also enabled us to validate the mechanical model and produce a set of anisotropic THM property values for both intact and damaged Boom Clay.
The Belgian radioactive waste management organization ONDRAF/NIRAS studies Boom Clay as a potential host rock for the disposal of high-level radioactive waste. Constructing a deep geological repository will lead to an excavation damaged/disturbed zone (EDZ/EdZ; Tsang et al. 2005) in the clay around the repository. In addition, coupled thermo-hydro-mechanical (THM) perturbations will be induced by the heat emitted by the high-level radioactive waste. These disturbances may affect the clay's ability to retard radionuclide transport. One of the challenges for the disposal of heat-emitting waste is to estimate the impact of the EDZ/EdZ and THM perturbations on the short- and long-term performance of the deep geological repository. This requires a good understanding of the processes that occur in the host formation and their evolution over time, for which research, development and demonstration have been going on in Belgium for more than 40 years.
To study the thermal impact of heat-emitting waste on the Boom Clay THM properties and processes, a large number of laboratory characterization programmes (e.g. Baldi et al. 1988; Sultan 1997; Delage et al. 2000; Le 2008; Li et al. 2011; Lima 2011) and a series of in-situ heater tests of different scales were conducted in the HADES Underground Research Laboratory (URL), which is excavated in the Boom Clay at a depth of 223 m in Mol, Belgium. The two most important heater tests are discussed in this paper: the small-scale ATLAS Heater Test and the large-scale PRACLAY Heater Test. The ATLAS Heater Test was carried out in four phases from 1992 to 2012 (De Bruyn and Labat 2002; Chen et al. 2011, 2023). The PRACLAY Heater Test has been running since November 2014 (Li et al. 2010; Van Marcke et al. 2013; Dizier et al. 2021; Chen et al. 2021a, b). Figure 1 shows the location of these two tests in the HADES URL.
The ATLAS and PRACLAY heater tests in the HADES Underground Research Laboratory (URL).
The ATLAS and PRACLAY heater tests in the HADES Underground Research Laboratory (URL).
The small-scale ATLAS Heater Test allows characterization of the far-field THM behaviour of a real repository. The heater diameter is small and the heating load is moderate and progressively changed, so the heater essentially functions as a line heat source and the THM disturbance of the clay is limited. In addition, the monitoring points are sufficiently far away from the heat source that the clay is less disturbed by the installation of the heater and thermal loading at the monitoring points. Therefore, the THM responses at the monitoring points mainly depend on and are more representative of the far-field THM behaviour of the Boom Clay, and the ATLAS Heater Test does not allow an insight to be gained into the Boom Clay THM behaviour in the near field of disposal galleries. The large-scale PRACLAY Heater Test aims to examine, at a scale and in conditions more representative of a repository, the combined impact of hydro-mechanical disturbances caused by the gallery construction and a thermal load on the Boom Clay. In this way, the PRACLAY Heater Test not only allows the gap left by the ATLAS Heater Test to be filled and provides insight into the effect of the temperature on the EDZ, but it also offers an opportunity to confirm and refine at a large scale our understanding of the Boom Clay THM behaviour obtained from the previous laboratory characterization programmes and small-scale in-situ tests.
This paper focuses on the numerical interpretation of the ATLAS and PRACLAY heater tests to confirm and refine the anisotropic THM parameter values and the mechanical model of the Boom Clay on the basis of the existing knowledge previously obtained from the laboratory characterization programmes and prior numerical modelling of in-situ heater tests. The present numerical modelling interprets the measured pore pressure and temperature from the most recent phase of the ATLAS Heater Test (i.e. ATLAS IV) and from the PRACLAY Heater Test which has been going on for 5 years. The present interpretation of both tests was realized in a stepwise and combined way: it was first performed for the small-scale ATLAS IV Heater Test, and then up-scaled to the large-scale PRACLAY Heater Test to investigate the up-scaling effect on the Boom Clay behaviour.
This paper first introduces the history, design, setup, phases, instrumentation, thermal loading and main observations of the small-scale ATLAS Heater Test and the large-scale PRACLAY Heater Test, respectively. Then it presents the combined numerical interpretation of the two tests. It outlines the theoretical framework of the modelling, provides a brief overview of the THM property values of the Boom Clay and presents the numerical modelling carried out for both tests. This also includes an illustration of how these simulations supported the interpretation of the test results and how they enabled derivation of the THM parameter values. The paper ends with conclusions and perspectives.
Two in-situ heater tests
The small-scale ATLAS Heater Test
Between 1993 and 2012, the small-scale ATLAS Heater Test was carried out in four phases (i.e. ATLAS I–IV). The goal of the test was to assess the hydro-mechanical effects of the thermal transient generated by the high-level radioactive waste on the Boom Clay. The ATLAS Heater Test provided valuable insights to characterize the THM behaviour of Boom Clay, in particular its far-field behaviour, and it significantly supported the preparation, prediction and interpretation of the large-scale PRACLAY Heater Test.
Setup and phases
The original test setup for the ATLAS Heater Test was established in 1992 by SCK CEN within the framework of the European project Interclay II (1990–1994; Knowles et al. 1996). The history of the ATLAS Heater Test is briefly summarized in Table 1, and the test setup is illustrated in Figure 2.
The ATLAS Heater Test setup (red segment in TD89E represents the ATLAS heater): (a) top view; and (b) front view.
The ATLAS Heater Test setup (red segment in TD89E represents the ATLAS heater): (a) top view; and (b) front view.
History of the ATLAS Heater Test
Year | Test setup and phases |
---|---|
1992 | Drilling and installation of the heating borehole (TD89E) and instrumented boreholes (TD85E and TD93E) |
1993–96 | First heating phase (ATLAS I) |
1996–97 | Second heating phase (ATLAS II) |
1997 | Power shutdown |
2006 | Drilling and installation of instrumented boreholes TD97E and TD98E |
2007–08 | Third heating phase (ATLAS III) |
2008 | Power shutdown |
2010 | Drilling and installation of instrumented borehole TD90Iu |
2011–12 | Fourth heating phase (ATLAS IV) |
2012 | Power shutdown |
Year | Test setup and phases |
---|---|
1992 | Drilling and installation of the heating borehole (TD89E) and instrumented boreholes (TD85E and TD93E) |
1993–96 | First heating phase (ATLAS I) |
1996–97 | Second heating phase (ATLAS II) |
1997 | Power shutdown |
2006 | Drilling and installation of instrumented boreholes TD97E and TD98E |
2007–08 | Third heating phase (ATLAS III) |
2008 | Power shutdown |
2010 | Drilling and installation of instrumented borehole TD90Iu |
2011–12 | Fourth heating phase (ATLAS IV) |
2012 | Power shutdown |
The 19 m long ATLAS heater borehole TD89E with a diameter of 230 mm was drilled in 1992 from the Test Drift (Fig. 1). The borehole has a steel casing with an external diameter of 190 mm and an internal diameter of 160 mm. The 8 m long heater was installed at a distance between 11 and 19 m. Four U-profiled sections are attached to the inner wall of the borehole casing, which are evenly distributed over the perimeter. Grooved aluminium strips containing the heater cables are fitted into these sections.
The test setup for the first two phases, ATLAS I (1993–96) and ATLAS II (1996–97), consists of two horizontal observation boreholes (TD85E and TD93E). In 2006, the setup was expanded by drilling two additional observations boreholes: one horizontal (TD98E) and one inclined downward (TD97E). The heater was subsequently switched on from April 2007 to April 2008. This heating phase is called ATLAS III. To further examine the anisotropic THM response, an additional inclined observation borehole TD90Iu was drilled above the heater borehole TD89E in 2010. The fourth heating phase of the ATLAS Heater Test, ATLAS IV, lasted from 2011 to 2012.
Instrumentation
Table 2 gives a brief summary of the sensors in the five observation boreholes. These boreholes allow monitoring of the evolution of the temperature and pore water pressure in the clay up to a radius of c. 4.5 m from the heater.
Summary of the six boreholes and main instrumented sensors for ATLAS Heater Test
Borehole | Casing diameters (mm) | Material | Sensor type | Distance to Test Drift lining intrados (m) | Sensor code |
---|---|---|---|---|---|
TD89E | 190/160 | Stainless steel | |||
TD85E | 89/78 | Stainless steel | Thermistor | 14.82 | Th-TD85E |
Piezometer | 14.62 | PP-TD85E | |||
TD93E | 89/78 | Stainless steel | Thermistor | 14.86 | Th-TD93E |
Piezometer | 14.66 | PP-TD93E | |||
TD98E | 101.6/94.4 | Stainless steel | Thermocouple | 20/17/16/14/13/10/9 | TC-TD98E-1/3/4/6/7/9/10 |
Thermocouple | 19/15/11 | TC-TD98E-2/5/8 | |||
Piezometer | 19/15/11 | PP-TD98E-1/2/3 | |||
TD97E | 80/40 | PVC | Thermocouple | 21–10 (evenly distributed) | TC-TD97E-1-12 |
TD90Iu | 88.9/81.7 | Stainless steel | Thermocouple | 23/19/15/11/7 | TC-TD90Iu-1/2/3/4/5 |
Piezometer | 23/19/15/11/7 | PP-TD90Iu-1/2/3/4/5 |
Borehole | Casing diameters (mm) | Material | Sensor type | Distance to Test Drift lining intrados (m) | Sensor code |
---|---|---|---|---|---|
TD89E | 190/160 | Stainless steel | |||
TD85E | 89/78 | Stainless steel | Thermistor | 14.82 | Th-TD85E |
Piezometer | 14.62 | PP-TD85E | |||
TD93E | 89/78 | Stainless steel | Thermistor | 14.86 | Th-TD93E |
Piezometer | 14.66 | PP-TD93E | |||
TD98E | 101.6/94.4 | Stainless steel | Thermocouple | 20/17/16/14/13/10/9 | TC-TD98E-1/3/4/6/7/9/10 |
Thermocouple | 19/15/11 | TC-TD98E-2/5/8 | |||
Piezometer | 19/15/11 | PP-TD98E-1/2/3 | |||
TD97E | 80/40 | PVC | Thermocouple | 21–10 (evenly distributed) | TC-TD97E-1-12 |
TD90Iu | 88.9/81.7 | Stainless steel | Thermocouple | 23/19/15/11/7 | TC-TD90Iu-1/2/3/4/5 |
Piezometer | 23/19/15/11/7 | PP-TD90Iu-1/2/3/4/5 |
The two observation boreholes, TD85E and TD93E, are drilled at either side of the heater borehole TD89E. At the end of both boreholes, sensors were installed at the same depth as the heater centre to measure the temperature and pore water pressure in the horizontal plane.
Borehole TD97E is 21 m long and has a downward inclination of 10° and a horizontal deviation of 10° to the left. It has 12 temperature sensors which are evenly distributed along depths of 10–21 m. Borehole TD98E is 20 m long and is located at a horizontal distance of about 2.66 m from the central borehole. It has 10 temperature sensors distributed along depths of 9–20 m, and three pore water pressure sensors at depths of 11, 15 and 19 m, respectively. Borehole TD90Iu is 23 m long and has an upward inclination of 8.6°. It has five pore water pressure sensors integrated with five thermocouples.
Thermal loading
The heating power during ATLAS I and ATLAS II is shown in Figure 3a. A constant heating power of 900 W was applied in ATLAS I from July 1993 to June 1996. The power was then increased to 1800 W for ATLAS II. It was maintained at that level from June 1996 to May 1997 (De Bruyn and Labat 2002). This was followed by the shutdown of the heater in May 1997 and the subsequent natural cooling.
The heating power during the four phases of the ATLAS Heater Test: (a) ATLAS I and II; and (b) ATLAS III and IV.
The heating power during the four phases of the ATLAS Heater Test: (a) ATLAS I and II; and (b) ATLAS III and IV.
The heating power during ATLAS III and IV is shown in Figure 3b. ATLAS III ran from April 2007 to April 2008 and ATLAS IV from October 2011 to November 2012. Both tests used nearly identical heating and cooling cycles in order to facilitate the interpretation and comparison of the results. The heating power was increased incrementally from 400 to 900 to 1400 W, after which the test was stopped. This means that the Boom Clay was subjected to a higher thermal load and for a longer time during ATLAS II than during the subsequent ATLAS III and ATLAS IV tests. The main difference between ATLAS III and ATLAS IV is the duration of the heating step of 1400 W. This lasted 27 days longer in ATLAS IV.
Main observations from ATLAS IV
The findings from ATLAS I-III can be found in François et al. (2009) and Chen et al. (2011). This section focuses on the observations of ATLAS IV.
Figure 4 shows the temperature evolution over time during ATLAS IV at five sensors located in the mid plane perpendicular to the heater and crossing the heater centre. These sensors have a radial distance from the heater centre ranging from 1.34 to 3.34 m, with a maximum temperature variation of 7–25°C, respectively. The anisotropic thermal behaviour of the Boom Clay can be directly observed using two sensors: TC-TD98E-5 and TC-TD90Iu-3. The former is in the horizontal direction, with a radius of 2.66 m to the heater centre, while the latter is at a similar distance (2.64 m) in the vertical direction. The measured temperature in the horizontal direction is higher than that in the vertical direction, with a maximum difference of up to 3.5°C (i.e. 35% of the maximum temperature variation at TC-TD90Iu-3).
Temperature variation at five sensors in the mid plane of the ATLAS heater.
Temperature variation at five sensors in the mid plane of the ATLAS heater.
Figure 5a and b shows the pore water pressure variations measured in ATLAS IV, respectively, at five sensors in the same horizontal plane as the heater and five above the heater. Similar to ATLAS III (Chen et al. 2011), a temporary decrease in the pore water pressure at the start of each heating step was observed in ATLAS IV at pore water pressure sensors in the horizontal plane. The opposite response, a temporary increase, was seen when the heater was switched off. Different observations were made from sensors above the heater from Figure 5b:
For each step that the heating power increased, the pore water pressure increased immediately and it increased at a higher rate than it did at the end of the previous heating step.
After the power was switched off, the pore water pressure at three sensors above the heater decreased immediately.
Pore water pressure evolution in ATLAS IV Heater Test: (a) at five sensors in the same horizontal plane as the heater; and (b) at five sensors above the heater.
Pore water pressure evolution in ATLAS IV Heater Test: (a) at five sensors in the same horizontal plane as the heater; and (b) at five sensors above the heater.
These differences in the pore water pressure response also emerged in numerical simulations performed for ATLAS III. They result from the instantaneous mechanical reaction owing to mechanical anisotropy and hydro-mechanical coupling (Chen et al. 2011). Unfortunately, this behaviour that was observed in the numerical simulations could not be verified in ATLAS III owing to the lack of sensors above or below the heater. With the pore water pressure sensors installed in the new upward borehole drilled for ATLAS IV, field measurements could confirm this phenomenon, which was predicted by numerical simulations.
ATLAS III and ATLAS IV use almost identical heating and cooling cycles, which facilitates comparison of their measurements. Figure 6a compares pore water pressure variation at three representative sensors (i.e. PP-TD93E, PP-TD98E-1 and PP-TD98E-2), and Figure 6b compares temperature variation between ATLAS III and ATLAS IV at three representative sensors (i.e. Th-TD93E, TC-TD98E-5 and TC-TD97E-6). ATLAS IV experienced a power failure of 12 h after heating for about 210 days, and its heating step of 1400 W is 27 days longer than that of ATLAS III, without which identical responses of both temperature and pore water pressure can be expected. This interesting observation indicates that THM disturbances from ATLAS III and ATLAS IV had no irreversible impact on the Boom Clay behaviour, and the Boom Clay behaviour can reasonably be assumed to be elastic during ATLAS III and ATLAS IV.
Comparison of pore water pressure and temperature between ATLAS III and ATLAS IV measured at three sensors: (a) pore water pressure; and (b) temperature.
Comparison of pore water pressure and temperature between ATLAS III and ATLAS IV measured at three sensors: (a) pore water pressure; and (b) temperature.
The large-scale PRACLAY Heater Test
The main goal of the PRACLAY Heater Test is to investigate the combined impacts on the Boom Clay of hydro-mechanical disturbances caused by gallery construction and a large-scale thermal load that mimics the heat emitted by high-level radioactive waste.
Design
Since it is not possible to reproduce the timescale (hundreds or thousands of years) and the spatial scale (hundreds of metres) of the thermal and hydraulic conditions of a real repository (Sillen and Marivoet 2007, Dizier 2020), some basic design criteria were defined by scoping calculations (Li et al. 2010):
Heating a 30 m long section of a gallery is estimated to be sufficiently long for the middle of that section to be representative of the actual repository.
The PRACLAY Heater Test was designed to be carried out such in a way that the THM conditions are more penalizing than expected in a real repository. To reach this goal, the thermal and hydraulic boundary conditions are controlled as follows:
The temperature at the interface between the Boom Clay and gallery lining along the heated section is kept at a constant temperature (target temperature). Both the target temperature and the temperature increase rate to reach this target temperature have to be higher than those in an actual repository.
A quasi-undrained hydraulic boundary for the Boom Clay along the heated gallery has to be applied to generate more penalizing hydro-mechanical conditions in the Boom Clay.
Setup and phases
Figure 7 shows a general layout of the setup in the PG. The setup mainly consists of the gallery lining, a seal including an insulation door, the heaters and backfill sand. A brief introduction of the test setup is given below. More details can be found in Van Marcke et al. (2013) and Dizier et al. (2016).
The gallery is supported by three kinds of linings:
a segmented concrete lining, with an internal diameter of 1.9 m, and an external diameter of 2.5 m. The lining consists of 81 concrete lining rings; the rings are 0.5 m long and consist of nine segments each (eight normal segments, S1–S8, and one key segment, S9);
the steel tunnelling shield that remained in place at the end of the construction work to support the 2.5 m long gallery section between the last lining ring and the end plug; and
the concrete end plug supporting the end of the PG.
The PG is heated using a primary heater consisting of electrical cables placed at a distance of 10 cm from the gallery intrados (the inner side of the gallery lining). A secondary heater is placed in the centre of the gallery and acts as a backup for the primary heater. That primary heater is split into three zones that can be controlled independently. To efficiently transfer the heat from the heating elements towards the clay massif, the heated section of the gallery was backfilled with water-saturated sand before the heater was switched on. To minimize the heat loss through the seal and to have a better control of the thermal boundary condition during the heating phase, a thermal insulation door was installed in front of the seal.
The main test phases are the following (see also Fig. 8):
construction of the PG in October and November 2007;
the waiting phase during which pore water drainage into the PG occurred (about 4 years);
saturation and pressurization of the 33.5 m long backfilled section of the PG to 1.0 MPa from December 2011 to November 2014 (about 3 years);
start-up of the heating phase from November 2014 to August 2015 (about 9.5 months):
250 W m−1 from 3 November 2014 to 7 January 2015 (65 days);
350 W m−1 from 7 January 2015 to 3 March 2015 (55 days)
450 W m−1 from 3 March 2015 to 17 August 2015 (167 days)
Stationary heating phase for 10 years from August 2015 until August 2025 during which a constant temperature of 80°C is kept at the lining extrados (the outer side of the gallery lining)
Instrumentation
A large instrumentation network with about 1100 sensors was set up around the PG to monitor the THM response in the main test components, especially in the Boom Clay. Several boreholes from the Connecting Gallery and the PG were instrumented with piezometers (pore water pressure monitoring) complemented with thermocouples (temperature monitoring; Fig. 9). Several concrete lining rings were instrumented as well to monitor the temperature inside the lining (e.g. temperature at extrados of rings 37, 50, 55 and 81, shown in Fig. 7).
Layout of the monitoring boreholes around the PRACLAY Gallery. The sensor numbering in the boreholes starts from the end.
Layout of the monitoring boreholes around the PRACLAY Gallery. The sensor numbering in the boreholes starts from the end.
Thermal loading
The heating power is applied to the three zones of the primary heater located c. 10 cm from the inside of the gallery lining (see Fig. 7). During the stationary heating phase, the heating power is controlled by the following two temperature indicators:
temperature indicator 1, Tind1 – the average of the temperature at the extrados of rings 37, 50 and 55, all located in zone 2 (zone 2 is the longest of the three zones, accounting for 88% of the total length);
temperature indicator 2, Tind2 – the average of the temperature at the extrados of ring 81 located in zone 3.
Heater power in the three zones of the primary heater in the PRACLAY Gallery.
Heater power in the three zones of the primary heater in the PRACLAY Gallery.
Two temperature indicators and pore water pressure in PRACLAY Gallery (see the definitions of Tind1 and Tind2 in this section, P_PG is the pore pressure in the PRACLAY gallery, PG).
Two temperature indicators and pore water pressure in PRACLAY Gallery (see the definitions of Tind1 and Tind2 in this section, P_PG is the pore pressure in the PRACLAY gallery, PG).
From then on, the stationary heating phase started. The power in zone 3 was first kept at 450 W m−1 until 6 June 2016, when Tind2 reached 80°C as well. The heating power in all zones was then decreased stepwise to maintain both temperature indicators at 80°C.
Observations
As observed in Figure 11, the target temperature of 80°C was reached at the interface between the gallery lining and the Boom Clay within about 9.5 months. Owing to the rapid increase in temperature and the high permeability of the sand in the backfilled gallery, the pore water pressure in the gallery increased rapidly and homogeneously from 1.0 to 2.9 MPa during this start-up heating phase. To maintain the target temperature during the stationary heating phase, the heating power was reduced stepwise. This induced a small drop in the pore water pressure followed by quite a stable pore water pressure evolution.
Because of the low permeability and small deformation of the Boom Clay, the heat dissipates primarily through conduction. This causes a temperature gradient in the Boom Clay. Figure 12a shows the temperature variation along four horizontal boreholes after 4 years of stationary heating. The temperature variation reached 50.9°C in the middle of CG35E (0.75 m from PG) and 5.7°C at CG49E (14.75 m from PG).
Variation of temperature and pore water pressure along four horizontal boreholes after 4 years of stationary heating: (a) temperature variation; and (b) pore water pressure variation.
Variation of temperature and pore water pressure along four horizontal boreholes after 4 years of stationary heating: (a) temperature variation; and (b) pore water pressure variation.
Because the thermal expansion coefficient of pore water is considerably higher than that of the Boom Clay's solid skeleton, the temperature gradient induces excess pore water pressures and a hydraulic gradient in the clay. Figure 12b shows the thermally induced pore water pressure variation along four horizontal boreholes after 4 years of stationary heating. The gradient in pore water pressure variation over a distance of 14 m is about 1.2 MPa (between CG35E and CG49E). At 14.75 m from the PG, the pore pressure was disturbed by about 0.57 MPa.
Similar to the anisotropic thermal behaviour of the Boom Clay observed in ATLAS III and ATLAS IV, the thermal anisotropy is clearly observed in the PRACLAY Heater Test. Figure 13a shows the measured temperature profiles in the mid plane of the PRACLAY Heater Test. The temperature in the horizontal direction is measured in four CG boreholes (CG35E, CG38E, CG42E and CG49E) and the temperature in the vertical direction in borehole PG50D (see Fig. 9 for both directions). The measured temperature profile in the horizontal direction is always higher than that in the vertical direction. This confirms that the Boom Clay thermal conductivity parallel to the bedding plane is higher than that perpendicular to the bedding plane.
Thermal anisotropy in mid plane (some measurements are not available owing to sensor failure): (a) radial profiles of temperature in mid plane; and (b) time evolution of temperature in the mid plane.
Thermal anisotropy in mid plane (some measurements are not available owing to sensor failure): (a) radial profiles of temperature in mid plane; and (b) time evolution of temperature in the mid plane.
In the mid plane of the PRACLAY Heater Test, two pairs of sensors (i.e. pair 1, CG35E-6 and PG50D-9; pair 2, CG42E-2 and PG50D-5) are selected to study the thermal anisotropy. The two sensors in each pair are more or less at the same radial distance from the heater, but one sensor lies in a horizontal plane while the other lies in a vertical plane with the heater. Figure 13b presents the evolution of temperature over time at these sensors. For pair 1, the measured temperature at CG35E-6 is higher than that at PG50D-9 although the latter is 0.25 m closer to the lining extrados. Since the boundary temperature at the interface Boom Clay/lining is isotropic and both sensors of pair 1 are close to the interface, the thermal anisotropy at this pair is not so significant. The two temperature sensors of pair 2 have almost the same radius, and a maximum temperature difference up to 5.0°C is observed. This again demonstrates the thermal anisotropy of the Boom Clay.
Figure 14a presents both horizontal and vertical profiles of the pore water pressure measured in the mid plane of the PRACLAY Heater Test at three different times: just before heating, at the end of the start-up heating and after 4 years of stationary heating. An anisotropic response can be clearly observed at three different times.
Anisotropic response of pore water pressure in mid plane: (a) radial profiles of pore water pressure in mid plane; and (b) time evolution of pore water pressure in mid plane.
Anisotropic response of pore water pressure in mid plane: (a) radial profiles of pore water pressure in mid plane; and (b) time evolution of pore water pressure in mid plane.
Figure 14b presents the pore water pressure evolution over time at three pairs of sensors in the mid plane of the PRACLAY Heater Test. These include the aforementioned two pairs (i.e. CG35E-6 and PG50D-9; CG42E-2 and PG50D-5) and a third pair (CG55E-4 and CG30Iu-3) at a distance of 20 m from the heater, but again in different directions to the heater.
Before switching on the heater, the backfilled PG was pressurized from atmospheric pressure to 1.0 MPa. During this pressurization period, the pore water pressure in the Boom Clay was measured by sensors below and above the backfilled gallery. They showed that the pore water pressures increased at different rates depending on their distances from the gallery. At sensors in the same horizontal plane as the backfilled gallery, pore water pressures close to the gallery (e.g. 2.24 m) behaved in the same way as those at the sensors in a vertical plane. The pore water pressures at sensors far away from gallery (≥8.97 m), however, kept decreasing.
After the heater was switched on, the pore water pressures at sensors below and above the heater continued to increase at higher rates. At sensors in the same horizontal plane far away from gallery (≥8.97 m), the pore water pressures continued to decrease at a higher rate. This decrease took longer for sensors further away from the gallery. It is remarkable that although CG30Iu-3 and CG55E-4 (pair 3) are located about 20 m from the gallery, their pore water pressures responded almost immediately to the switch-on of the heater. These observations are consistent with those made from the smaller-scale ATLAS Heater Test, which can be interpreted by the anisotropic mechanical behaviour of Boom Clay and hydro-mechanical coupling.
Total pressure sensors are installed at the ends of some boreholes around the PG. They measure the stress in four directions: above, below, to the left and to the right of the casing. Figure 15a shows the total pressures measured at the end of CG42E, which is 9.0 m from the gallery axis and in the same horizontal plane as the heater. Figure 15b gives the total pressures measured at the end of CG30E, which is 3.5 m below the gallery axis. The following can be observed at the start of heating:
The vertical total pressures measured in CG42E decreased immediately after heating, but began to rise after a few months. The horizontal total pressures measured here increased immediately after heating.
The vertical total pressures measured in CG30E increased immediately after heating, but the horizontal total pressures decreased immediately after heating and then began to increase after a few weeks.
Anisotropic responses of total pressures: (a) CG42E (in the horizontal plane); and (b) CG30E (below the PRACLAY Gallery).
Anisotropic responses of total pressures: (a) CG42E (in the horizontal plane); and (b) CG30E (below the PRACLAY Gallery).
The total stress sensors located closer to the heater were thermally disturbed earlier, resulting in an earlier thermally induced stress increase. This explains why the duration of the decrease at CG30E is shorter than that at CG42E. The anisotropic response of the total stress observed at the two borehole ends in the PRACLAY Heater Test confirms the anisotropic mechanical behaviour of the Boom Clay.
Combined numerical interpretation of two in-situ heater tests
This section presents the numerical modelling carried out for both heater tests. A brief account of the balance equations and the Boom Clay THM properties is given for completeness. A more detailed description can be found elsewhere (Olivella et al. 1996; Villar et al. 2020; Chen et al. 2021a; Li et al. 2022).
Balance equations
The solution of the coupled THM problem requires the simultaneous solution of the balance equations. The balance equations are established for the porous medium as a whole. The materials involved in the numerical models are considered as two-phase saturated porous media.
The state variables in equations (1)–(3) are the temperature T, the liquid pressure Pl and the solid displacement vector u.
Boom Clay thermo-hydraulic parameter values
The Boom Clay is a cross-anisotropic material: its main THM parameter values parallel to the bedding plane (horizontal plane in this study) are higher than those perpendicular to the bedding plane (vertical plane in this study). The value of specific heat C is reported to be in a range of 1400–1550 J (kg K)−1 (Buyens and Put 1984; Djeran et al. 1994; Van Cauteren 1994) with low uncertainty (Garitte et al. 2014). The thermal conductivity values of the Boom Clay have been estimated from laboratory tests (Djeran et al. 1994; Lima 2011; Dao 2015) and from numerical thermal interpretation of the in-situ heater tests in HADES URL (Chen et al. 2011, 2023; Garitte et al. 2014). The thermal conductivity parallel to the bedding plane (λh) ranged approximatively from 1.34 to 1.90 W (m K)−1 and that perpendicular to the bedding plane (λv) from 0.82 to 1.31 W (m K)−1. More information about the thermal parameter values can be found in Chen et al. (2023) and Li et al. (2022). The values {λh, λv} = {1.9, 1.2} W (m K)−1 will be used in this study because they were derived from a large amount of temperature measurements from these two tests. In addition, its equivalent value: , is close to those derived from the other heater tests in HADES URL (i.e. CERBERUS and CACTUS; Van Cauteren 1994) and well-controlled laboratory tests. Over the years, the intrinsic permeability K of the Boom Clay around the HADES URL has been measured from a large number of laboratory tests on core samples and in-situ tests. The critical review of these laboratory and in-situ hydraulic conductivity measurements carried out by Yu et al. (2013) shows that nearly all the K values measured for the Boom Clay are in the range of approximately 1.5 × 10−19 to 8.0 × 10−19 m2 with an anisotropic ratio Kh/Kv value of about 2.5. The intrinsic permeability values used or derived from the prior and present modelling fall in this range.
Boom Clay mechanical constitutive laws and parameter values
The Boom Clay mechanical model used in the numerical interpretation of both tests is an elasto-plasticity model based on Drucker–Prager yield criteria, and the elastic deformation is calculated by Young's elastic modulus and Poisson's ratio.
Yield criteria, plastic potential function and hardening law
Elasticity
The elastic deformations are usually calculated based on the Hook's law of linear elasticity. In the case of cross-anisotropic conditions, the linear elasticity is defined by five independent parameters:
Young's elastic modulus in the bedding plane, Eh;
Young's elastic modulus in the plane perpendicular to the bedding plane, Ev;
shear elastic modulus in the plane perpendicular to the bedding plane, Gv;
Poisson's ratio in the bedding plane, νhh;
Poisson's ratio for the effect of horizontal strain on vertical strain νvh.
Typical stiffness variation with the strain and approximate strain limits for reliable measurement (Mair 1993 and Clayton 2011).
Typical stiffness variation with the strain and approximate strain limits for reliable measurement (Mair 1993 and Clayton 2011).
The Boom Clay is characterized by a highly non-linear stress–strain response with a very clear trend of stiffness variation with strain level (Villar et al. 2020; Li et al. 2022). Using the bender element method, Lima (2011) and Dao (2015) investigated the anisotropic shear modulus of the natural far-field Boom Clay samples. Their measured values vary between 900 and 1450 MPa depending on the direction relative to the bedding plane. These values are close to the magnitude of the shear modulus based on the in-situ seismic measurements performed in the HADES URL (Demanet et al. 1994; Schuster et al. 2001). These values could represent the Boom Clay stiffness in the far field with very small strain.
Triaxial laboratory tests on the Boom Clay resulted in tangent modulus at large strain with a distributed range of 75–351 MPa (Ouvry 1983; Horseman et al. 1987; Coll 2005; Le 2008; Charlier et al. 2010; Villar et al. 2020). During the excavation of the galleries in the HADES URL, large deformations were measured in the near field. Based on in-situ measurements made during the excavation of the Test Drift, the drained Boom Clay Young's modulus was estimated to be approximately equal to 300 MPa (Mair et al. 1992). This value was confirmed by the back analysis of the convergence measured during the excavation of the Connecting Gallery (Bernier et al. 2007) and by laboratory tests (Charlier et al. 2010). These values are most likely representative of the elastic behaviour of the Boom Clay in the near field with large strain.
The numerical modelling code
The computer code CODE_BRIGHT is used for all the numerical modelling introduced in this paper. The code was developed by the Technical University of Catalonia (Olivella et al. 1996). The code takes into account all the major THM processes, balance equations, constitutive laws and equilibrium constraints in the system of both heater tests. In this code, the system of partial differential equations is discretized in space (finite elements) and time (finite difference). The discretization in time is linear and the implicit scheme uses two intermediate points, tk+ε and tk+θ between the initial tk and final tk+1 times. Finally, since the problems are non-linear, the Newton–Raphson method is applied to solve the problem in an iterative way.
Numerical modelling of the ATLAS Heater Test
The small-scale ATLAS Heater Test provides an ideal test setup for characterizing the Boom Clay THM properties as it has a simple geometry, well-defined boundary conditions and simple test phases. A summary of the main features of the prior and the present modelling of the ATLAS I-IV Heater Test can be found in Table 3.
A summary of the main features of the numerical modelling performed for ATLAS and PRACLAY Heater tests
Prior modelling | Present modelling | ||||||||
---|---|---|---|---|---|---|---|---|---|
François et al. (2009) | Chen et al. (2011) | Chen et al. (2021a) | Chen et al. (2021b) | Chen et al. (2023) | This study | ||||
Heater test | ATLAS I&II | ATLAS III | PRACLAY | PRACLAY | ATLAS IV | PRACLAY | ATLAS IV | PRACLAY | |
THM model | 2D-axis THM | 3D THM | 2D-PS THM | 2D-axis THM | 3D T | 3D THM | 2D-PS THM | ||
Initial stress | Isotropic | Anisotropic | Anisotropic | Isotropic | – | Anisotropic | Anisotropic | ||
THM property | Isotropic | Anisotropic | Anisotropic | Isotropic | – | Anisotropic | Anisotropic | ||
Excavation modelling | No | No | Yes | Yes | – | No | Yes | ||
Thermal B.C. | Temperature | Heater power | Heater power | Heater power | Heater power | Temperature | Heater power | Temperature | |
Main THM parameter values | T | λ = 1.7 W (mK)−1 | {λh, λv} = {1.65, 1.31} W (m K)−1 | λ = 1.53 W (m K)−1 | {λh, λv = } = {1.9, 1.2} W (m K)−1 | {λh, λv} = {1.9, 1.2} W (m K)−1 | |||
H | K = 2.0 × 10−19 m2 | {Kh, Kv} = {4.0, 2.0} × 10−19 m2 | {Kh, Kv} = {6.0, 3.0} × 10−19 m2 | K = 4.5 × 10−19 m2 | – | {Kh, Kv} = {3.0, 1.5} × 10−19 m2 | |||
αl* | 3.5 × 10−4 C−1 | 3.0 × 10−4 C−1 | αl (T) | – | αl (T) | ||||
M | ACMEG-T elasto-thermoplasticity | Drucker–Prager yield criteria | – | Drucker–Prager yield criteria | |||||
{Eh, Ev}EDZ&far-field = {1400, 700} MPa | {Eh, Ev}EDZ = {400, 200} MPa, {Eh, Ev}far-field = {1400, 700} MPa | EEDZ = 300 MPa, Efar-field = 1050 MPa | – | {Eh, Ev}EDZ&far-field = {Eh(εi), Ev(εi)} | |||||
Agreement with measurements | Fairly good | Good | Overestimate | Reasonably good | Good | Good | |||
Main limitations | Without considering anisotropic stress and THM properties; αl is assumed constant | Both αl and Young's modulus are assumed constant | Using power B.C. in this model; Two distinct E values are used for EDZ and far-field | Without considering anisotropic stress and THM properties; Two distinct E values are used for EDZ and far field | Thermal numerical interpretation can be improved by considering various uncertainties | PRACLAY modelling can be improved by a 3D THM modelling to study the impact of connecting gallery |
Prior modelling | Present modelling | ||||||||
---|---|---|---|---|---|---|---|---|---|
François et al. (2009) | Chen et al. (2011) | Chen et al. (2021a) | Chen et al. (2021b) | Chen et al. (2023) | This study | ||||
Heater test | ATLAS I&II | ATLAS III | PRACLAY | PRACLAY | ATLAS IV | PRACLAY | ATLAS IV | PRACLAY | |
THM model | 2D-axis THM | 3D THM | 2D-PS THM | 2D-axis THM | 3D T | 3D THM | 2D-PS THM | ||
Initial stress | Isotropic | Anisotropic | Anisotropic | Isotropic | – | Anisotropic | Anisotropic | ||
THM property | Isotropic | Anisotropic | Anisotropic | Isotropic | – | Anisotropic | Anisotropic | ||
Excavation modelling | No | No | Yes | Yes | – | No | Yes | ||
Thermal B.C. | Temperature | Heater power | Heater power | Heater power | Heater power | Temperature | Heater power | Temperature | |
Main THM parameter values | T | λ = 1.7 W (mK)−1 | {λh, λv} = {1.65, 1.31} W (m K)−1 | λ = 1.53 W (m K)−1 | {λh, λv = } = {1.9, 1.2} W (m K)−1 | {λh, λv} = {1.9, 1.2} W (m K)−1 | |||
H | K = 2.0 × 10−19 m2 | {Kh, Kv} = {4.0, 2.0} × 10−19 m2 | {Kh, Kv} = {6.0, 3.0} × 10−19 m2 | K = 4.5 × 10−19 m2 | – | {Kh, Kv} = {3.0, 1.5} × 10−19 m2 | |||
αl* | 3.5 × 10−4 C−1 | 3.0 × 10−4 C−1 | αl (T) | – | αl (T) | ||||
M | ACMEG-T elasto-thermoplasticity | Drucker–Prager yield criteria | – | Drucker–Prager yield criteria | |||||
{Eh, Ev}EDZ&far-field = {1400, 700} MPa | {Eh, Ev}EDZ = {400, 200} MPa, {Eh, Ev}far-field = {1400, 700} MPa | EEDZ = 300 MPa, Efar-field = 1050 MPa | – | {Eh, Ev}EDZ&far-field = {Eh(εi), Ev(εi)} | |||||
Agreement with measurements | Fairly good | Good | Overestimate | Reasonably good | Good | Good | |||
Main limitations | Without considering anisotropic stress and THM properties; αl is assumed constant | Both αl and Young's modulus are assumed constant | Using power B.C. in this model; Two distinct E values are used for EDZ and far-field | Without considering anisotropic stress and THM properties; Two distinct E values are used for EDZ and far field | Thermal numerical interpretation can be improved by considering various uncertainties | PRACLAY modelling can be improved by a 3D THM modelling to study the impact of connecting gallery |
*αl is the volumetric thermal expansion coefficient of pore water.
Prior numerical modelling of ATLAS I-IV
For ATLAS I and II, François et al. (2009) used the ACMEG-T elasto-thermoplastic constitutive model for the Boom Clay, using both a 2D plane strain and a 2D axisymmetric THM model. The authors concluded that a 2D axisymmetric model is needed to account for axial flow and axial deformation. The 2D axisymmetric model, however, does not take account of THM anisotropy. Furthermore, the authors highlighted the significant role of thermo-plasticity on the global THM response of the clay formation, as shown by the laboratory tests (Baldi et al. 1988; Sultan 1997). The model did not reproduce the limited measurements from ATLAS I and II with reasonable accuracy.
Chen et al. (2011) modelled ATLAS III using a 3D coupled THM model and a Drucker–Prager elasto-plasticity model taking into account the cross-anisotropic initial stresses and THM properties of the Boom Clay. The model successfully reproduced measurements at 24 temperature sensors and five pore water pressure sensors. This allowed the derivation of a set of cross-anisotropic THM parameter values (see Table 3) and confirmed the THM anisotropy of the Boom Clay. However, the model did not take into account the temperature dependency of the thermal expansion coefficient of the pore water (αl) and assumed a constant value for this parameter for ATLAS I, II and III.
Recently, a purely thermal 3D modelling of ATLAS IV was performed by Chen et al. (2023). Thanks to an additional inclined observation borehole and very accurate estimations of the thermocouple positions, more precise thermal conductivity values could be obtained: {λh, λv} = {1.9, 1.2} W (m K)−1. This study highlights that the derived thermal conductivity values depend on the position of the sensors with respect to the heater. The aforementioned limitations and the updated thermal conductivity will be taken into account in the 3D coupled THM modelling for ATLAS IV which is discussed hereafter.
3D THM model of ATLAS IV
The 3D THM model developed for ATLAS IV is very similar to the one developed for ATLAS III (Chen et al. 2011), except that the former considers the thermal expansion coefficient of the pore water as temperature dependent, as suggested for Boom Clay pore water by Baldi et al. (1988), Delage et al. (2000), Ghabezloo and Sulem (2009) and Chen et al. (2021a, b). For the sake of completeness and for convenience of the reader, the main features of the model are briefly presented below.
The simulated domain measures 100 m in the radial direction and 119 m in the axial direction of the central heater (Fig. 17a). The 19 m long steel tube in the central heater borehole is included in the geometry, and the heater is attached to its innermost part, which is 8 m long. Compared with ATLAS III, ATLAS IV uses a finer mesh with 17412 nodes and 15 264 hexahedral elements (Fig. 17b).
(a) Geometry, materials and initial conditions; and (b) mesh of 3D thermo-hydro-mechanical (THM) model for ATLAS IV.
(a) Geometry, materials and initial conditions; and (b) mesh of 3D thermo-hydro-mechanical (THM) model for ATLAS IV.
The conditions at the depth of the HADES URL – an in-situ temperature of 16.0–16.5°C, a pore water pressure of 2.2–2.3 MPa and a vertical stress of 4.5 MPa at the far-field Boom Clay – are given (Bernier et al. 2007; Charlier et al. 2010; Dizier et al. 2016; Chen et al. 2021a, b). Since the model is mainly intended to interpret the variation of temperature and pore pressure within a limited distance above and below URL, a uniform initial temperature of 16.5°C, a pore water pressure of 2.25 MPa and a vertical total stress of 4.5 MPa with a coefficient of earth pressure at rest equal to 0.7 are assumed for the whole model domain (Fig. 17a).
All the boundary surfaces are assumed to be adiabatic and impermeable. ATLAS IV started by re-activating the heater on 18 October 2011. The power was first increased to 400 W and then stepwise to 900 and 1400 W. On 29 November 2012, it was instantaneously shut down to observe the effects of the cooling transient (see Fig. 3b).
The mechanical boundary conditions are defined as follows (see Fig. 17):
the displacement in the Z-direction is fixed on the surface AOCE;
the displacement in the X-direction is fixed on the surface AOBD;
the displacements in both the X- and Z-directions are fixed on the surface BCED;
the displacement in the Y-direction is fixed on the surface OCB; and
a normal boundary stress of 3.825 MPa is applied on the surface ADE.
The model developed for ATLAS IV uses the temperature-dependent volumetric thermal expansion coefficient αl for the pore liquid and the refined Boom Clay thermal conductivity values {λh, λv} = {1.90, 1.2} W (m K)−1. The other main THM parameter values are discussed below.
Boom Clay THM behaviour based on parametric sensitivity analysis of ATLAS IV
A parametric sensitivity analysis can support the interpretation of the in-situ measurements and help to identify the key parameters. During the ATLAS I and ATLAS II tests, the Boom Clay around the heater had been disturbed by a higher and longer thermal loading than during the ATLAS III and ATLAS IV tests. Therefore, no thermo-plasticity and only reversible thermal and mechanical strain are expected in ATLAS III and ATLAS IV (Cui et al. 2000). This hypothesis is confirmed by Figure 6: the pore water pressure variations are the same during ATLAS III and IV and the Boom Clay behaviour is only thermo-elastic during both tests. Therefore, one can assume that the pore water pressures during ATLAS IV are primarily controlled by the following THM parameters:
temperature-dependent volumetric thermal expansion coefficient αl(T) of the pore water, as in Chen et al. (2021a, b);
cross-anisotropic thermal conductivity as estimated by Chen et al. (2023) – {λh, λv} = {1.90, 1.2} W (m K)−1;
cross-anisotropic permeability values Kh and Kv;
five independent cross-anisotropic elasticity parameter values – Eh, Ev, Gv, νhh and νvh.
For ATLAS IV, a large number of numerical modelling cases were calculated (Chen et al. 2021a). The results from four representative cases (see Table 4) are chosen and compared with measurements at three filters: TD93E, TD98E-P2 and TD90Iu-P3 (see Fig. 18). These three filters are all located in the mid plane, which is perpendicular to the heater and at the centre of the heater. The former two filters are in the same horizontal plane as the heater while the third one is placed above the heater.
Sensitivity analysis using 3D THM model for ATLAS IV: (a) Case 1 (base case), (b) Case 2 (Kh/Kv = 3.0/1.5 × 10−19 m2), (c) Case 3 (Eh/Ev/Gv = 1600/800/320 MPa) and (d) Case 4 (Eh0/Ev0/Gv0 = 1600/800/440 MPa, modified Boom Clay elastoplastic model).
Sensitivity analysis using 3D THM model for ATLAS IV: (a) Case 1 (base case), (b) Case 2 (Kh/Kv = 3.0/1.5 × 10−19 m2), (c) Case 3 (Eh/Ev/Gv = 1600/800/320 MPa) and (d) Case 4 (Eh0/Ev0/Gv0 = 1600/800/440 MPa, modified Boom Clay elastoplastic model).
Four representative cases for the sensitivity analysis of ATLAS IV
Case no. | Kh (10−19 m2) | Kv (10−19 m2) | Eh (MPa) | Ev (MPa) | Gv (MPa) | νhh (–) | νvh (–) |
---|---|---|---|---|---|---|---|
1 (Base) | 4.0 | 2.0 | 1400 | 700 | 280 | 0.25 | 0.125 |
2 | 3.0 | 1.5 | 1400 | 700 | 280 | 0.25 | 0.125 |
3 | 4.0 | 2.0 | 1600 | 800 | 320 | 0.25 | 0.125 |
4 | 3.0 | 1.5 | 1600 (Eh0) 320 (Ehf) | 800 (Ev0) 160 (Evf) | 440 (Gv0) 88 (Gvf) | 0.25 | 0.125 |
Case no. | Kh (10−19 m2) | Kv (10−19 m2) | Eh (MPa) | Ev (MPa) | Gv (MPa) | νhh (–) | νvh (–) |
---|---|---|---|---|---|---|---|
1 (Base) | 4.0 | 2.0 | 1400 | 700 | 280 | 0.25 | 0.125 |
2 | 3.0 | 1.5 | 1400 | 700 | 280 | 0.25 | 0.125 |
3 | 4.0 | 2.0 | 1600 | 800 | 320 | 0.25 | 0.125 |
4 | 3.0 | 1.5 | 1600 (Eh0) 320 (Ehf) | 800 (Ev0) 160 (Evf) | 440 (Gv0) 88 (Gvf) | 0.25 | 0.125 |
The hydro-mechanical parameter values for the base case are taken from Table 3. They were calibrated using 24 temperature sensors and five pore water pressure sensors from ATLAS III, using a constant value of αl = 3.0 × 10−4°C−1. When combining this set of hydro-mechanical (HM) parameters with a more realistic temperature-dependent αl(T) in Case 1, the pore water pressure variations are obviously underestimated by the model (Fig. 18a). This indicates the significant impact of αl(T) on the pore water pressure and highlights the importance of re-evaluating the Boom Clay THM parameter values derived from ATLAS III.
Case 2 differs from Case 1 only in terms of the intrinsic permeability. Figure 18b shows some discrepancy between the observations and the modelled pore water pressure variation of Case 2, particularly during the cooling phase. When comparing the modelling results of Case 1 and Case 2, it becomes clear that the pore water pressure variation is highly dependent on the permeability.
The modelled pore water pressure variations during both the heating and cooling phases increase as Eh, Ev and Gv increase from 1400, 700 and 280 MPa (Case 1) to 1600, 800 and 320 MPa (Case 3), respectively. This is shown in Figure 18a for Case 1 and in Figure 18c for Case 3. A clear discrepancy between the observations and the outcomes of the Case 3 modelling can be seen in Figure 18c. Based on the modelling results from Cases 1–3, one can observe that the model using constant elastic moduli values does not allow proper capture of the pore water pressure variation during both the heating and cooling phases.
Case 4 adopts the Boom Clay mechanical model with a non-linear elasticity modulus proposed in equation (8) to interpret the measured pore water pressures. By applying trial and error on a large number of values of three cross-anisotropic moduli (Eh, Ev and Gv), it was found that when Eh0, Ev0 and Gv0 are equal to 1600, 800 and 440 MPa, and Ehf, Evf and Gvf are equal to 320, 160 and 88 MPa, the model outcomes match well with the measurements at the three filters (Fig. 18d). Furthermore, the pore water pressure variations at all 10 filters were calculated using the model and the associated parameters of Case 4. Figure 19 shows these calculations in comparison with measurements. Figure 19a presents the pore water pressure variations at the five filters in three horizontal boreholes (TD85E, TD93E and TD98E), and Figure 19b shows the pore water pressure variations at the five filters in the inclined upward borehole (TD90Iu). A good agreement of the pore water pressures, in terms of both trend and magnitude, is obtained for all the filters. The set of cross-anisotropic THM parameter values used in Case 4 for the Boom Clay is summarized in Table 5. These fall in the ranges outlined in the sections ‘Boom Clay thermo-hydraulic parameter values’ and ‘Boom Clay mechanical constitutive laws and parameter values’.
Comparison of modelled (full lines) and measured (dashed lines) pore water pressure variation for Case 4: (a) pore water pressure variation at TD85E, TD93E and TD98E; and (b) pore water pressure variation at TD90Iu.
Comparison of modelled (full lines) and measured (dashed lines) pore water pressure variation for Case 4: (a) pore water pressure variation at TD85E, TD93E and TD98E; and (b) pore water pressure variation at TD90Iu.
Anisotropic thermo-hydro-mechanical parameter values for a non-linear elastoplastic model of the Boom Clay
Thermal conductivity (W m−1 K) | Intrinsic permeability (10−19 m2) | Initial and final elastic modulus (MPa) | |||||||
---|---|---|---|---|---|---|---|---|---|
λh | λv | Kh | Kv | Eh0 | Ev0 | Gv0 | Ehf | Evf | Gvf |
1.90 | 1.20 | 3.0 | 1.5 | 1600 | 800 | 440 | 320 | 160 | 88 |
Thermal conductivity (W m−1 K) | Intrinsic permeability (10−19 m2) | Initial and final elastic modulus (MPa) | |||||||
---|---|---|---|---|---|---|---|---|---|
λh | λv | Kh | Kv | Eh0 | Ev0 | Gv0 | Ehf | Evf | Gvf |
1.90 | 1.20 | 3.0 | 1.5 | 1600 | 800 | 440 | 320 | 160 | 88 |
Numerical modelling of the PRACLAY Heater Test
The geometry, materials, conditions and test phases of the PRACLAY Heater Test are more complex than those of the ATLAS Heater Test. Representing this complexity with sufficient accuracy by a full 3D-coupled THM modelling is computationally expensive. Therefore, for the time being, the THM modelling of PRACLAY Heater Test is limited to two-dimensional approaches. A summary of the main features of the prior and the present modelling of the PRACLAY Heater test can be found in Table 3.
Prior numerical modelling
Before switching on the PRACLAY heater, blind prediction modelling was performed by both a 2D plane strain (2D-PS) THM model and a 2D axisymmetric (2D-Axis) THM model (Chen et al. 2021a). The geometry of the 2D-PS model is a middle cross-section of the heated gallery and Boom Clay, perpendicular to the PRACLAY axis, which is representative of the mid plane of the PRACLAY Heater Test (Dizier et al. 2016). This model has the advantage of accounting for the anisotropic initial stress and anisotropic THM behaviour of the Boom Clay. However, it cannot consider heat flux, liquid flow and displacement in the direction parallel to the PG axis. When the heater power of zone 2 as shown in Figure 10 is used as the thermal boundary condition, the model overestimated the temperature and pore water pressure, in particular in the long term.
Owing to the axisymmetric nature of the test setup geometry around the PRACLAY axis, a 2D-Axis THM model allows all main materials to be considered, as well as the axial heat flux, fluid flow and deformation. Unfortunately, the anisotropic initial stress and anisotropic THM behaviour of the Boom Clay cannot be examined. In the PRACLAY Heater Test, the temperature at the interface between the lining and Boom Clay is designed and controlled to be isotropic. The pore water pressure at the interface is also isotropic owing to the highly permeable backfill sand and lining. Because of these isotropic thermo-hydraulic (TH) boundary conditions, the anisotropic effects on the TH responses of the Boom Clay in the near field of the PG and during the first years of heating are limited. As a first order of approximation, the anisotropy of the initial stress and the THM behaviour of the Boom Clay can be neglected when modelling the first years of heating (Chen et al. 2021b).
The prior 2D-Axis modelling of the PRACLAY Heater Test showed that the general Boom Clay THM responses in the PRACLAY Heater Test are well captured by the numerical blind predictions, confirming our already good knowledge of the Boom Clay THM behaviour (Chen et al. 2021a, b). Yet in the aforementioned modelling, two distinct elastic modulus values were assigned to the near- and far-field Boom Clay in the Drucker–Prager elasto-plasticity model, which is not entirely representative of the physical reality. This limitation will be addressed by adopting the non-linear elasticity given in equation (8), which could represent the continuous alteration of the elastic modulus from the near to the far field.
Recently, a purely thermal 3D modelling has also been performed for the PRACLAY Heater Test (Chen et al. 2023) using the same pair of thermal conductivity values {λh, λv} = {1.9, 1.2} W (m K)−1 derived from the thermal modelling of ATLAS IV. The temperatures measured within the EDZ as well as in the far field from the PRACLAY Heater Test were well reproduced using this pair of thermal conductivity values. This indicates that, in terms of heat transport properties, the impact of the excavation and THM perturbation on the thermal conductivity of the Boom Clay around PRACLAY is not significant.
By addressing the above limitations and using the anisotropic THM property values of Boom Clay derived from the ATLAS Heater Test, a 2D-PS THM model described below will focus on the interpretation of the observations in the mid plane of the PRACLAY Heater Test.
2D-PS THM model
The geometry of the 2D-PS model is a middle cross-section of the heated PG and Boom Clay perpendicular to the PG axis. Without taking into account the gravitational effect and owing to the symmetrical nature of the problem, only a quarter of the entire cross-section is modelled. The simulated region measures 100 m in both directions. Three materials – the Boom Clay, concrete lining and backfill sand – are included in the model. The main test phases considered in this model were introduced in the previous section ‘Setup and phases’.
Figure 20 shows the geometry and materials of the model. It assumes a 200 m thick Boom Clay formation. This is larger than the actual thickness of 100 m at the Mol site. This means that the aquifers above and below the Boom Clay are not taken into account. The geometry is discretized by four-node quadrilateral structured elements and the mesh has 3201 nodes and 3075 elements with reduced/selective integration and modified B-bar (B is the strain–displacement matrix used in the finite element method, Olivella et al. 1996).
Geometry, materials, initial and boundary conditions in the 2D plane strain coupled THM model.
Geometry, materials, initial and boundary conditions in the 2D plane strain coupled THM model.
All materials are assumed to have an initial temperature of 16.5°C. The Boom Clay has an initial pore water pressure of 2.3 MPa, while the lining and backfill sand have an initial pore water pressure of 0.1 MPa. The initial vertical total stress in the Boom Clay perpendicular to the bedding plane is 4.5 MPa with a coefficient of earth pressure at rest equal to 0.7. The initial stresses in all the other materials are isotropic, with a value of 0.1 MPa.
The temperature along the boundaries AB and BC (Fig. 20) is fixed at 16.5°C because they are far enough from the heater and will not be disturbed during the 10 years when the heater test is running. The symmetrical boundaries OC and OA can be assumed adiabatic. The pore water pressure at the boundaries AB and BC is assumed constant at 2.3 MPa and the symmetrical boundaries OC and OA can be assumed to be impermeable.
During the excavation, it is assumed that the pore water pressure at the excavation surface decreases from 2.3 MPa to atmospheric pressure (i.e. 0.1 MPa). Before the start of the heater test, the pore water pressure in the PG was first artificially increased (by injecting water) from 0.1 to 0.5 MPa over a period of three months and from 0.5 to 1.0 MPa in the following 2 years and 5 months.
The top boundary AB is subjected to a vertical total stress of 4.5 MPa and the lateral boundary BC to a horizontal total stress of 3.825 MPa. The symmetric boundary OA is horizontally fixed (no horizontal displacement) and the bottom symmetric boundary OC is vertically fixed (no vertical displacement). A radial convergence of 6 cm is assumed to take place linearly during the excavation period. Gravity is not taken into account.
By using the measured temperatures at the interface between the Boom Clay and the lining as the thermal boundary condition and the measured pore water pressure inside the PG as the hydraulic boundary condition, the modelled temperature and pore water pressure in the Boom Clay essentially depend on the THM behaviour of Boom Clay. This reduces the potential impact of the model's inherent deficiency to consider heat flux, liquid flow and displacement parallel to the PG axis. This remains consistent with the goal of the PRACLAY Heater Test to characterize the Boom Clay THM behaviour.
The Drucker–Prager elasto-plasticity model is applied for the Boom Clay using the non-linear elasticity modulus proposed in equation (8). The Boom Clay THM parameter values derived from ATLAS IV (see Table 5) are used.
Main modelling results compared with measurements in the mid plane
The modelling results are compared with the measurements in the mid plane of the PRACLAY Heater Test. The thermal anisotropy that could be observed in the measurements (see Fig. 13) can be compared with the modelling results, shown in Figures 21a and 21b respectively. The good agreement at all sensors and at different times validates the set of thermal conductivity values {λh, λv} = {1.90, 1.20} W (m K)−1 derived from the 3D thermal interpretation of ATLAS IV (Chen et al. 2023).
Comparison between the modelled and measured temperature evolution in Boom Clay (some sensors stopped working due to failure, radial distances from sensors to gallery axis are indicated in the figures): (a) radial profiles of temperature in mid plane; and (b) time evolution of temperature in mid plane.
Comparison between the modelled and measured temperature evolution in Boom Clay (some sensors stopped working due to failure, radial distances from sensors to gallery axis are indicated in the figures): (a) radial profiles of temperature in mid plane; and (b) time evolution of temperature in mid plane.
Figures 22a and 22b compare the measured pore water pressure profiles presented in Figure 14a with the modelling results in the horizontal and vertical directions at three different times. Two cases of modelling results are presented in Figures 22a and 22b. The case with the legend ‘Model with E(ε)’ corresponds to the 2D-PS THM model adopting non-linear elasticity given in equation (8). The case with the legend ‘Model with 2 const. E’ corresponds to the case using two distinct sets of elasticity values for the near and far-field Boom Clay. The separation between the near and far field is set at a radius of 10.8 m (Dizier et al. 2016). Both the trends of horizontal and vertical profiles at three different times are captured reasonably well, although the modelled results do not perfectly match the measurements. For the modelling case with legend ‘Model with 2 const. E’, a non-smooth transition of modelled pore water pressure profiles around the radius of 10.8 m can be observed at the end of the start-up heating phase. This is due to the use of two distinct sets of elasticity values for the near and far field. By using the non-linear elasticity in equation (8) in the modelling case ‘Model with E(ε)’, smooth profiles can be observed at all times.
Comparison between modelled (dashed lines) and measured (solid lines) pore water pressure profiles in Boom Clay: (a) horizontal direction; and (b) vertical direction.
Comparison between modelled (dashed lines) and measured (solid lines) pore water pressure profiles in Boom Clay: (a) horizontal direction; and (b) vertical direction.
Figure 23a and b compares the measurements presented in Figure 14b with the corresponding modelling results in the horizontal and vertical directions. Although there is no perfect match with the measurements, the modelling results capture the observed evolution and magnitude of pore water pressures before and during heating at all locations reasonably well. No significant difference between the two modelling cases can be found from Figure 23a and b.
Comparison between modelled (full lines) and measured (dashed lines) pore water pressure evolution in Boom Clay (radial distances from sensors to gallery axis are indicated in the figures): (a) sensors in horizontal direction; and (b) sensors in vertical direction.
Comparison between modelled (full lines) and measured (dashed lines) pore water pressure evolution in Boom Clay (radial distances from sensors to gallery axis are indicated in the figures): (a) sensors in horizontal direction; and (b) sensors in vertical direction.
Using this 2D-PS THM model, the measured temperature in the mid plane of the PRACLAY Heater Test is reproduced excellently. The evolution and the magnitude of measured pore water pressure in the mid plane of the PRACLAY Heater are captured reasonably well. This validates the proposed mechanical model and the set of anisotropic THM parameter values for the Boom Clay derived from ATLAS IV.
Discussion
In this section, the relevance of the modified elastoplastic model and the anisotropic HM parameter values (see Table 5) derived from the small-scale ATLAS IV Heater Test and later confirmed by the large-scale PRACLAY Heater Test are further tested by modelling some triaxial tests performed on the Boom Clay samples. Table 6 presents the basic test conditions of three representative drained triaxial tests on the Boom Clay samples taken at the depth of the HADES URL (Coll 2005; Le 2008). The tests were performed at the room temperature and the samples were consolidated under a stress of 2.3–2.5 MPa.
Three laboratory tests on the Boom Clay samples taken at the depth of the HADES Underground Research Laboratory
Sample no. | Sample size D×h (mm) | Initial void ratio e0 | Confining pressure σ′3 (MPa) | Shear rate (µm min−1) | References |
---|---|---|---|---|---|
B07 | 39.9 × 40.6 | 0.582 | 2.3 | 0.25 | Coll (2005) |
B20 | 39.9 × 39.6 | 0.559 | 2.3 | 0.25 | Coll (2005) |
E09 | 38.4 × 74.5 | 0.61 | 2.5 | 1 | Le (2008) |
Sample no. | Sample size D×h (mm) | Initial void ratio e0 | Confining pressure σ′3 (MPa) | Shear rate (µm min−1) | References |
---|---|---|---|---|---|
B07 | 39.9 × 40.6 | 0.582 | 2.3 | 0.25 | Coll (2005) |
B20 | 39.9 × 39.6 | 0.559 | 2.3 | 0.25 | Coll (2005) |
E09 | 38.4 × 74.5 | 0.61 | 2.5 | 1 | Le (2008) |
These tests were modelled using two cases of elasticity parameter values. Case 1 used the modified mechanical model and the set of elasticity parameter values presented in Table 5, while Case 2 used the constant elasticity parameter values Ehf, Evf and Gvf given in Table 5. Figure 24 compares the modelled deviatoric stress with the laboratory test results.
Comparison of the deviatoric stress between modelling and the laboratory test results.
Comparison of the deviatoric stress between modelling and the laboratory test results.
The modelled curves from Case 1 present two distinct stages. During the initial stage with small axial strain, the sample behaves quite stiffly and remains in the elastic state with the decrease in the elastic modulus. During the later stage the sample shows highly non-linear stress–strain behaviour and experiences both modulus decrease and plastic deformation.
In general, the numerical modelling reproduces the trend and peak value of the measured deviatoric stress from laboratory tests reasonably well. Modelling Case 1 presents a stiffer initial response to shearing than the modelling Case 2 and the laboratory tests (Fig. 24). This discrepancy between modelling Case 1 and the laboratory test results could be partly explained by the disturbance of the test samples by core drilling and sample trimming. This disturbance history is not taken into account in modelling Case 1. Modelling Case 2 reproduces the initial response in the laboratory test results better. This is because it uses the constant and final elasticity modulus values which correspond to values when the Boom Clay is subject to sufficient strain.
The reasonably good reproduction and interpretation of the laboratory test results strengthen our confidence in the proposed elastoplastic model with non-linear and cross-anisotropic elasticity. It also confirms the parameter values obtained from the ATLAS IV and PRACLAY heater tests.
Conclusions and perspectives
In the HADES URL, the small-scale ATLAS Heater Test ran in four phases from 1993 to 2012, and the large-scale PRACLAY Heater Test has been running since November 2014. These tests are conducted to examine how the Boom Clay responds to the heat generated by high-level radioactive waste. The specific goals of these tests are to characterize the THM behaviour of the Boom Clay and confirm or refine the THM values obtained from laboratory tests. Numerical simulations have supported interpreting the test results and enabled to derive THM parameter values.
The small-scale ATLAS IV Heater Test is first analysed by a 3D coupled THM model. An extensive sensitivity analysis was performed demonstrating the model capability to accurately reproduce both the magnitude and the trend of the pore water pressure evolution at all filters during the entire heating and cooling cycle. These excellent results validate the proposed mechanical model for the Boom Clay and allowed the derivation of a set of anisotropic THM parameter values for the Boom Clay.
The cross-anisotropic THM characterization results of the Boom Clay from the small-scale ATLAS Heater Test are then used to support the numerical interpretation of the measurements in the mid plane of the large-scale PRACLAY Heater Test. This was done using a 2D-PS coupled THM model. This model sets the measured temperature at lining extrados in the mid plane as the thermal boundary condition and the measured pore water pressure in the PG at the lining intrados as the hydraulic boundary condition. The model provides results that closely match the measurements within both the EDZ and in the far field.
The mechanical model and the anisotropic HM parameter values derived from the small-scale ATLAS IV Heater Test and later confirmed by the large-scale PRACLAY Heater Test are further reaffirmed by reproducing both the general trend and the peak value of the laboratory test on small Boom Clay samples. It shows that measurements from both in-situ tests and laboratory tests can be well reproduced and interpreted using the unified models and associated parameter values for the Boom Clay. This strengthens our confidence in the Boom Clay THM characterization results from this study.
The anisotropic yielding behaviour of the Boom Clay demonstrated by experimental results is not introduced in the present modelling of the two heater tests. This could be further studied in the future by referring to earlier research on the anisotropic yielding behaviour of the Boom Clay (Sultan et al. 2010; François et al. 2014). More data will be generated in the years to come. This will allow further refining the THM characterization of the Boom Clay, particularly of its far field. A full 3D THM model which can account for the impacts of both THM anisotropy and Connecting Gallery is necessary to interpret the test response in the long term and far field.
Acknowledgements
The authors would like to express their appreciation to all colleagues who have contributed over the past 40 years to the ATLAS and PRACLAY Heater tests. Special thanks go to the present and former colleagues of the ESV EURIDICE technical team that installed these setups in the HADES URL.
Competing interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper, except from the agencies to which the authors are affiliated.
Author contributions
GC: investigation (lead), writing – original draft (lead); XL: investigation (equal), writing – review & editing (equal); AD: investigation (equal), writing – review & editing (equal); JV: investigation (equal), writing – review & editing (equal); XS: investigation (equal), writing – review & editing (equal); SL: investigation (equal), writing – review & editing (equal).
Funding
This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.
Data availability
The data that support the findings of this study are available from EIG EURIDICE, SCK CEN and NIRAS/ONDRAF but restrictions apply to the availability of these data, which were used under licence for the current study, and so are not publicly available. Data are, however, available from the authors upon reasonable request and with permission of EIG EURIDICE, SCK CEN and NIRAS/ONDRAF.