1. Introduction
Deep geological disposal is widely accepted to be the safest and most efficient form of disposal for high-level radioactive waste [1]. In the deep geological disposal, highlevel radioactive waste such as spent nuclear fuels is typically disposed into the bedrock at a depth of about 500 m, e.g., in Swedish SKB concept [2]. The bedrock isolates the waste from the human environment as a natural barrier. The deposition holes are filled with bentonite buffer around emplaced canisters, and disposal tunnels are backfilled. Potential rock types suitable for the host rocks of a geological repository include crystalline rock, clay, and rock salt [3]. It is necessary to ensure the stability of the disposal system for tens of thousands of years considering the half-lives of the wastes. In South Korea, the KAERI Reference Disposal System for Spent Nuclear Fuel (KRS) was developed in 2007 [4] based on the Swedish KBS-3 concept. Subsequently, KRS for High Burnup Spent Nuclear Fuel (KRSHB) and the Improved Korean Reference Disposal System (KRS+) were suggested by [5] and [6], respectively.
Coupled thermo-hydraulic-mechanical (THM) processes are the key processes that ensure the long-term performance of the deep disposal of radioactive waste. The heat released from the high-level radioactive waste is the source of heat transfer by conduction and convection, fluid flow and phases, and thermal stress in the geological repository. The increase in temperature can change the saturation and generate thermal stress in the buffer. The changes in stress and water saturation can influence the material properties, such as permeability and thermal conductivity as well as thermal, hydraulic, and mechanical processes at the nearfield rock mass [7].
Field experiments and numerical analyses have been conducted to study the coupled THM processes at the near-field rock mass of deep disposal systems over the last several decades [8]. For example, the Drift Scale Test at Yucca Mountain, USA was carried out with a heating system of canister heaters and wing heaters from 1997–2006. A numerical study was performed to predict the field measurements of thermal, hydraulic, mechanical, and chemical responses [9] as well as the thermally induced changes in rock permeability and deformation [10, 11]. The Full-scale Engineered Barriers Experiment (FEBEX) experiment at the Grisel Test Site in Switzerland started in 1994 with the aim of demonstrating the Spanish engineered barrier system in crystalline rock. Coupled THM numerical analyses of the FEBEX experiment were carried out with many simulators together with comparison studies [12-14]. In Korea, numerical studies on the coupled THM processes have been conducted using the TOUGH2-MP/FLAC3D simulator [15]. Recent work includes the assessment of the thermal behavior according to disposal spacings, THM responses in multi-layer repositories and KRS+, and the Barcelona Basic Model (BBM) [16-19].
The physical properties of underground rock masses influence the thermal, hydraulic, and mechanical behavior and their interactions with one another. It is essential to understand the effect of rock mass properties on coupled THM responses for reliable assessment of the long-term performance of the deep geological repository. In this study, the effects of rock mass properties on the thermal, hydraulic, and mechanical responses were evaluated by numerical simulations using TOUGH-FLAC. A reference numerical model was established by benchmarking the heater test at the Kamaishi mine in Japan which has a wealth of monitoring data as well as extensive validation exercises employing various numerical methods [20-21]. The current study intends to complement previous numerical studies by considering additional sensitivity studies by using a different numerical technique using TOUGH-FLAC [21]. A benchmark sensitivity study was carried out with a range of properties, including the Young’s modulus, permeability, thermal conductivity, and thermal expansion coefficient of the rock mass, that encompass those of the candidate rock types for the host rock of the repository, including crystalline rock, clay, and rock salt. The responses used for analyses were temperature, water content, displacement, and stress at the buffer and near-field rock mass.
2. Numerical Model
2.1 Heater Test at the Kamaishi Mine
The heater test at the Kamaishi mine in Japan was selected as a reference model for this numerical study [20]. The heater test was conducted to test the THM responses of a hypothetical repository at a depth of about 250 m in a fractured granodiorite rock mass. The drift was excavated in 1994, and in the following year, a test pit with a 1.7 m diameter and 5 m depth was drilled. An electric heater was emplaced and surrounded by buffer material, which filled the test pit. A concrete lid was placed on the top of the test pit. A water pool constructed on the drift floor maintained a constant water head of 40 cm and a temperature of 12.3℃. The heating phase started in 1996 and lasted for about 260 days with a heater temperature of 100℃, followed by a 180-day natural cooling period. Monitoring sensors installed in the buffer and the near-field rock mass measured the THM responses, including temperature, water saturation, pore pressure, stress, and displacement.
During an international collaborative project titled DEvelopment of COupled models and their VALidation against EXperiments (DECOVALEX), the field measurements were numerically predicted by four research teams using different finite element codes for various THM responses, including temperature, moisture content, displacement, and pressure in the buffer and near-field rock mass [21]. The numerical results achieved reasonable agreement with the measured ones, although the prediction for total stress had discrepancies from the in situ measurements. The THM experiment in the Kamaishi mine was also used for calibrations of computer codes for one of the Benchmark Tests of the DECOVALEX III project [22]. In this study, a numerical analysis was performed with reference to the THM responses in the field experiment and numerical models presented in previous studies on the DECOVALEX project [21].
2.2 Numerical Simulator
The simulation was conducted using a THM simulator called TOUGH-FLAC [23]. The simulator coupled TOUGH2 [24], a fluid and heat flow simulator of multiphase, multicomponent systems in fractured porous media, and FLAC3D [25], a geomechanical simulator. The two simulators were sequentially coupled via exchanges of thermal, hydraulic, and mechanical responses at calculation steps. Hydraulic and thermal responses such as pore pressure, saturation, and temperature were calculated in TOUGH2 and delivered to FLAC3D. With the updated hydraulic and thermal variables, mechanical calculation was conducted in FLAC3D with consideration of effective stress, thermal stress, and swelling pressure. Through the application of the updated mechanical variables to TOUGH2, hydraulic and mechanical properties were recalculated. The TOUGHFLAC simulator has been applied to various geomechanical applications including nuclear waste disposals, geothermal reservoirs, CO2 sequestration, natural gas production, and energy storage systems [26-32].
The rock mass and buffer were assumed as continuum media in this study. Two-phase flow models described by Eqs. (1)–(6) were applied in TOUGH2 after calibration with the field observations and the numerical results given in [21]. The relative liquid (krl) and gas (krg) permeability can be expressed with effective saturation (S*) as Eqs. (1) and (2) [33]:
Effective saturation is defined as
where Sl and Slr are liquid saturation and residual liquid saturation, respectively. The intrinsic permeability of liquid (kliq) and gas (kgas) phases are assumed to be related as Eq. (4):
where b is the Klinkenberg parameter [34]. Capillary pressure (Pcap) is considered to follow the van Genuchten model [35]
where Pvan and λvan are the van Genuchten parameters, and Pmax is the maximum capillary pressure. The molecular diffusion is also taken into account with the diffusion coefficient for gases. The diffusion coefficient () of component κ in phase β is related to temperature and pressure:
where (P0, T0) is the coefficient in the standard condition of P0 = 1 atm and T0 = 0℃. In this study, the diffusion coefficients in the gas phase of water and air were considered and the coefficient in the standard pressure and temperature condition was 2.13 × 10−5 m2∙s−1 [36, 37]. Thermal conductivity (λ) was assumed to have a linear relation with liquid saturation according to
where λdry and λwet is dry and wet thermal conductivity, respectively. The thermal conductivities of the buffer material and rock mass were taken from the previous numerical work on this site [21].
The mechanical process was simulated in FLAC3D using an elastic model. The main mechanical components affected by coupling with thermal and hydraulic processes are thermal stress (thermal strain) and effective stress in the rock mass and buffer and swelling pressure in the buffer. The thermal stress (ΔσT) is related to temperature changes (ΔT) according to
where αT is linear thermal expansion coefficient, and K is bulk modulus. The effective stress (σ´) is obtained from Eq. (9):
where σ, α, and Pp are total stress, Biot coefficient, and pore pressure, respectively. A linear elastic swelling model was applied in the buffer to simulate swelling strain linearly that is dependent on saturation changes. Swelling stress was calculated as a function of the changes in water saturation using the following equation:
where Δσsw is induced swelling stress, σmax,sw is the maximum swelling stress, ΔSl is the change in water saturation, and Si is initial liquid saturation [31]. The parameters for the swelling stress function were determined based on the previous numerical models of the heater test [21]. Changes in thermal stress, effective stress, and swelling stress were calculated using changes in temperature, pore pressure, and saturation transferred from TOUGH2 at each time step.
2.3 Model Setup
A numerical model was constructed by referring to the Kamaishi field experiment and the previous numerical models in [21]. A quarter-symmetric model that included a heater, buffer, a lid, and rock mass was generated. The model was 10 m × 10 m × 20 m in size (Fig. 2). Roller boundary conditions were applied to the bottom and vertical symmetry planes to prevent normal displacement. Constant normal stresses were applied to the lateral outer boundaries, and a free boundary was applied to the top of the model. Hydrostatic pressure was assigned in the initial condition with a pore pressure of 4 kPa on the surface representing the water pool. The initial temperature was set to 12°C throughout the model. Constant pore pressure and temperature were applied on the boundaries, excepting the symmetry planes. Initially, the rock mass was almost fully saturated, and the buffer was partially saturated with a saturation of 0.62. The lid and heater were non-porous to prevent fluid flow, whereas they had elastic parameters set for mechanical processes.
In the numerical simulation, the temperature of the heater was increased to 100°C at the start of the heater test and maintained for 258 days to match the actual heater test. For the following 180 days of the cooling phase, the heater was naturally cooled down. Calibration of the reference model was carried out by matching the simulated results with the measured and predicted data at monitoring points given in [21]. The parameters for the reference case are presented in Table 1. For simplicity in the parametric study, parameters were constant over a simulation, independent from thermal, hydraulic, and mechanical responses.
A series of sensitivity studies were conducted to determine the effect of rock mass properties on the THM responses in the near-field rock mass. The properties were thermal conductivity, which is a thermal property; Young’s modulus, which is a mechanical property; thermal expansion coefficient, which is related to the coupled thermomechanical behavior; and permeability, which is a hydraulic property. Granite, rock salt, and clay were selected as the bedrock of the potential deep repository. The rock mass properties for the sensitivity study were determined as shown in Table 2, covering the property values of the candidate rock types. Each simulation was performed by changing only the value of the corresponding properties while keeping those of the other properties identical to those of the reference case. Three additional simulations that combined thermal conductivity, Young’s modulus, and thermal expansion coefficient according to the rock types were also conducted. The effects of the rock mass properties were analyzed by evaluating the temperature, water content, displacement, and stress at the measuring points in the near-field rock mass and buffer.
3. Results
3.1 Reference Case
The simulation results of the reference case were compared with those of the field experiment and previous modeling results presented in [21]. The previous modeling study was carried out by four research groups: Atomic Energy Control Board in Canada (AECB), Clay Technology in Sweden (CLAY), Lawrence Berkeley National Laboratory in the USA (LBNL) and a team in Japan noted as KIPH (Kyoto and Iwate Universities, Japan Nuclear Cycle Development Institute and Hazama Corporation). The responses included in the comparison were temperature and water content at monitoring points in the buffer and near-field rock mass and fluid pressure and vertical displacements in the rock mass.
The THM responses in the reference case generally agreed well with the field measurements and previous numerical results of the Kamaishi heater test in [21], as shown in Figs. 3–6. The temperature rose with the start of the heating and stabilized in about half a month at the selected point of the buffer (BT4) and in about a month at RH16 of the rock mass. The temperature then returned to the initial value in the cooling phase (Fig. 3). The interface between the buffer and rock mass (BW1 and BW3) was saturated due to the water inflow from the rock mass, and it stayed fully saturated throughout the experiment period (Fig. 4(a)). At the point of the buffer close to the interface with the heater (BW5), the water content decreased due to evaporation, and the vapor diffused outward due to the thermal gradient during the heating phase. The bufferheater interface was slowly resaturated during the cooling phase (Fig. 4(b)). The simulated pore pressure reached hydrostatic pressure at monitoring points in the near-field rock mass, as shown in Fig. 5, although it varied from the field measurement at RH13 in Fig. 5(a). Similar observations were also found in previous modeling exercises. Ground heaving of about 0.4 mm was simulated on the surface of KBM4 at the near-field rock mass, mainly due to by the thermal expansion. The simulated value was larger than the measured one, as in the previous modeling study, because the reference model applied the same elastic property values as those in the previous modeling work (Fig. 6).
3.2 Sensitivity Study
3.2.1 Thermal Conductivity
The changes in thermal conductivity of the rock mass led to changes in temperature, water content, stress evolution, and vertical displacement at the buffer and nearfield rock mass. The temperature stabilized at a higher value in the heating phase as the thermal conductivity decreased. The difference in the quasi-steady temperature between the cases of the largest and smallest thermal conductivity was 12℃ at BT4 in the buffer. It decreased as the distance from the heater increased, becoming about 5℃ at RH16 in the near-field rock mass. As the cooling started, temperature returned to the initial temperature more slowly with smaller thermal conductivity. The area affected by the temperature increase in the near-field rock mass was largest in the case of the lowest thermal conductivity (Fig. 7).
The water content evolution at the monitoring point of the buffer-rock interface (BW1) was similar between all three cases of thermal conductivity of the rock mass. Due to water inflowing from the saturated rock mass, the monitoring point, BW1, became fully saturated about one day after the heating started in all the cases. The duration for the saturation at BW1 was so short that temperature variations could not make a noticeable difference in the saturation time of this monitoring point (Fig. 8(a)). In contrast, the values of water contents varied with the different rock thermal conductivities at BW5. The trend of water content variation was found to be the same in all three cases.
During the heating process, the water content decreased due to the drying process and the diffusion from the heater to a further area according to the thermal gradient. As the temperature returned to the original values in the cooling period, the thermal diffusion became weaker and the water content slowly recovered. Although the overall behavior of the water content seemed similar, the temperature increased more in the buffer and the nearfield rock mass in the case of lower thermal conductivity of the rock mass. This increased the energy level of the system, which may have enabled a faster inflow of water from the saturated rock mass to the buffer system. When the increase of temperature near the heater is higher due to the lower thermal conductivity, the diffusion coefficient for gas would be larger as well, as explained by Eq. (6). This accelerated the outward vapor movement and facilitated greater water saturation near the heater. Therefore, the decrease in the rock thermal conductivity played a role in delaying the reduction of the water content near the heater in the buffer system (Fig. 8(b)).
The upward vertical displacement occurred at the nearfield rock mass mainly due to the thermal expansion caused by the temperature increase. The uplift was the largest at the end of the heating phase. On the monitoring borehole KBM4 located in the near-field rock mass, the vertical compressive stress was maximized at the depth of 3.0–3.4 m in all cases, which was around the depth of the center of the heater. The increase of the vertical stress was more significant when the rock mass had lower thermal conductivity. This was related to the higher temperature at the buffer and near-field rock mass when the rock thermal conductivity was higher. The largest ground heaving was 0.53 mm with the lowest thermal conductivity, and the smallest was 0.34 mm with the highest thermal conductivity. The displacement nearly returned to the original value in the cooling process. At the end of the cooling, slight ground heaving ranging from 0.005–0.038 mm still existed at KBM4 because of the remaining temperature increase by a few tenths of a degree and the swelling pressure due to the saturation of the buffer (Fig. 9).
3.2.2 Young’s Modulus
Changes in Young’s modulus of the rock mass mainly affected the mechanical behaviors, such as stress and displacement of the rock mass (Fig. 10). The maximum vertical compressive stress at KBM4 reduced from 6.1 MPa in the reference case to 3.6 MPa in Case 2-1 and 0.5 MPa in Case 2-2, appearing at the depth of 3.0–3.4 m around the mid-depth of the heater. The vertical displacement was dominantly affected by thermal expansion and was only slightly affected by Young’s modulus of rock mass. The ground heaving at the end of the heating period was predicted to be about 0.43 mm in the reference case and Case 2-1, and 0.45 mm in Case 2-2, showing a marginal difference. During the cooling period, the ground uplift on the top of KBM4 reversed and became close to the initial state, 0.02 mm in the reference case and Case 2-1 and 0.03 mm in Case 2-2. The thermal and hydraulic behavior showed no difference between the three cases because this modeling did not consider the mechanically induced changes in hydraulic and thermal properties.
3.2.3 Thermal Expansion Coefficient
The effect of the changes in the rock thermal expansion coefficient was obvious in the mechanical behavior of the near-field rock mass (Fig. 11). A large thermal expansion coefficient resulted in greater expansion of the near-field rock mass when the temperature changes were the same. Therefore, the vertical displacement, as well as the thermal stress, increased along KBM4 in the near-field rock mass with an increasing rock thermal expansion coefficient. The maximum vertical stress in the three cases ranged from 6.1–29.1 MPa, demonstrating a 4.8-fold difference. The ground heaving at the top of the monitoring line was 2.02 mm, 0.52 mm, and 0.43 mm in the reference case and Cases 3-1 and 3-2, respectively, at the end of the heating phase. Though the mechanical responses to the thermal expansion coefficients were apparent, the coefficients’ influence on the hydraulic and thermal responses was negligible. This was because this modeling did not consider the changes in hydraulic properties, such as permeability and porosity, caused by the poroelastic effects, and no mechanically induced changes in thermal properties occurred during the simulations.
3.2.4 Permeability
The permeability of the rock mass was related to the rate of water infiltration from the saturated rock mass to the buffer system. At the monitoring point of the buffer in the vicinity of the rock mass (BW1), the water content started to increase in both cases of rock permeability. Saturation occurred sooner in the case of higher rock permeability, and it reached full saturation about one day after the start of the experiment. The water content in the model with lower rock permeability, in contrast, developed more slowly, and the simulation stopped before the monitoring point attained full saturation (Fig. 12(a)). The water content at BW5 near the heater declined faster in the heating phase and recovered faster in the cooling phase for the higher rock permeability case than the lower rock permeability case. In the case of higher rock permeability, the buffer became more pressurized due to the faster water infiltration from the rock, which led to a smaller diffusion coefficient for gases and intrinsic gas permeability, as presented in Eqs. (6) and (4), respectively. This may have caused the slower outward movement of evaporated water from the heater-buffer interface, and consequently, a smaller saturation in the early period of the heating. However, in the long run, the amount of water inflow from the rock mass to the buffer was more significant in the lower rock permeability case, and the water content became higher than in the larger permeability case after 259 days (Fig. 12(b)).
The temperature showed a small discrepancy up to about 2 between the two cases at BT4 in the early heating phase (Fig. 12(c)). The slightly larger temperature in the low permeability case was attributed to the difference in buffer thermal conductivity due to the greater liquid saturation. The effect of rock permeability on temperature was nonetheless limited in terms of magnitude and period. The sensitivity of the rock permeability to mechanical responses was predicted to be insignificant along KBM4 in the rock mass. This was presumed to be because the thermoelastic process was more dominant in this model than hydro-mechanical processes, including swelling pressure and poroelastic stress.
3.3 Rock Types
Three additional simulations were performed by combining thermal conductivity, Young’s modulus, and thermal expansion coefficient of the rock mass according to the type of rock mass: crystalline rock in the Kamaishi mine, rock salt, and clay, as presented in Table 2. In all the cases, the rock permeability was assumed to be 1 × 10−16 m2, representing a fractured rock mass. Temperature variations at the monitoring points in the buffer and rock mass were almost identical to those in Fig. 7, which were from the case in which only thermal conductivity was varied from the reference case. This suggested great sensitivity of rock thermal conductivity to temperature and minor sensitivity of Young’s modulus and thermal expansion coefficient. The temperature increase was largest in clay, followed by crystalline rock, with rock salt showing the smallest temperature increase. The water content evolutions at BW1 and BW5 in the buffer were similar to the results shown in Fig. 8. This indicates that the mechanical properties, Young’s modulus, and thermal expansion coefficient played negligible roles in fluid flow. There was no distinction between the rock types in terms of the saturation rate at the rock-buffer interface. The degree of decrease in water content at BW5 was the largest in rock salt, followed by granite and clay.
The vertical stress and displacement along KBM4 in the near-field rock mass were influenced by all three rock mass parameters: thermal conductivity, Young’s modulus, and thermal expansion coefficient (Fig. 13). The vertical compressive stress was predicted to have the smallest value in the clay model because clay has a 10 times smaller Young’s modulus than rock salt and crystalline rock, despite having the smallest thermal conductivity. The greatest vertical stress was observed in the rock salt model due to the strong effect of the thermal expansion coefficient, even though its large thermal conductivity reduced its vertical stress to some extent. The smallest vertical displacement appeared in the reference model, as it had the smallest thermal expansion coefficient and the largest Young’s modulus. In addition, the crystalline rock’s larger thermal conductivity than clay’s induced the smaller vertical displacement of the reference model. The vertical displacement was the greatest in the rock salt model due to it having the largest thermal expansion coefficient.
4. Discussion
The sensitivity studies in the current study assumed that permeability and porosity were constant due to being independent of the THM responses. However, hydraulic properties are known to be affected by mechanical states of rock masses. For instance, an increase in effective stress leads to elastic deformation of rock pores and enhances rock porosity. Rock permeability can be also affected by changes in porosity and by rock failure induced by disturbance to the stress field. Changes in hydraulic properties influence the hydraulic response, and through that, thermal and mechanical responses. In this regard, changes in mechanical properties such as Young’s modulus and thermal coefficient generate changes in mechanical responses, and in turn, changes in hydraulic properties of the rock mass. Therefore, it is supposed that the application of appropriate hydraulic properties that are dependent upon mechanical behavior would enable the prediction of more intricate but accurate thermal, hydraulic, and mechanical responses.
Thermal properties are affected by not only porosity and saturation but also temperature. The thermal conductivity of rock is known to decrease with an increase in temperature, which also causes the specific heat and thermal expansion coefficient to increase [38]. During the simulation of the reference case, rock temperature increased to a maximum of 40℃ from the initial temperature of 12℃. If the heating period had lasted longer, mimicking the actual disposal of nuclear fuel, the rock temperature would have reached a much higher value. In that case, the long-term changes in thermal properties due to the high temperature would contribute to variation in temperature evolution, and accordingly, in hydro-mechanical responses as well.
Thermal stress had a dominant effect compared to swelling pressure and poroelastic stress on the mechanical system. In particular, in the buffer, the role of thermal stress was likely to be overestimated in this study due to the use of a relatively large thermal expansion coefficient, 1 × 10−4 [1/K], and a small value of maximum swelling pressure, 1 MPa. If different values for thermal expansion coefficient and maximum swelling pressure were applied to the buffer, it is expected that the buffer would have a different stress distribution, affecting the integrity at the buffer-rock interface. An auxiliary simulation with a thermal expansion coefficient of 1 × 10−6 [1/K] at the buffer revealed the build-up of thermal stress became less significant, and the swelling pressure and pore pressure were more dominant in the buffer and at the buffer-rock interface. However, the effect of the thermal expansion coefficient of the buffer diminished significantly in the rock mass. Further study is thus required for an in-depth evaluation of mechanical responses with respect to the buffer properties.
5. Conclusion
A series of benchmark sensitivity studies of the nearfield rock mass properties were performed to analyze the THM responses in a heater test. The thermal conductivity of the rock mass influenced the evolution of temperature, water content, stress, and displacement in the near-field rock mass and buffer. Young’s modulus and the thermal expansion coefficient of the rock had a large effect on the mechanical behavior but a negligible effect on the thermal and hydraulic behaviors. The rock permeability affected the rate of water inflow from the rock mass to the buffer, resulting in a disparity in water content and temperature between the different permeability cases.
Based on sensitivity studies using the properties of crystalline rock, rock salt, and clay, which are candidate rock types for the host rock for deep geological disposal, the temperature increases were largest in clay, followed by crystalline rock and rock salt, mainly because of the thermal conductivities of the rocks. The reduction of the water content near the heater in the buffer was thus delayed more in clay, followed by crystalline rock and rock salt, due to the difference in the thermal conductivities of the rocks. The vertical stress and displacement in the rock mass were predicted to be the greatest in the rock salt model due to it showing the largest thermal expansion coefficient. The vertical compressive stress was the smallest in the clay model because it had the smallest Young’s modulus.
This numerical study demonstrated that the material properties of rock affected the thermal, hydraulic, and mechanical behaviors in the near-field rock mass as well as the buffer. Depending on the bedrock types, the deep geological repository would have different coupled THM responses and, accordingly, different long-term performances. However, this study has room for improvement in terms of long-term rock mass properties related to THM behaviors. A close study of the rock properties in various THM environments is thus required in future research. This numerical study is important fundamental research that can be used for the construction of underground research laboratories and the design of geological repositories for high-level waste.