Abstract
Disposal of radioactive waste originating from reprocessing of spent research reactor fuel typically includes stainless steel canisters with waste immobilised in a glass matrix. In a deep borehole disposal concept, waste packages could be stacked in a disposal zone at a depth of one to potentially several kilometres. This waste will generate heat for several hundreds of years. The influence of combining a natural geothermal gradient with heat from decaying nuclear waste on radionuclide transport from deep disposal boreholes is studied by implementing a coupled heat-solute mass transport modelling framework, subjected to depth-dependent temperature, pressure, and viscosity profiles. Several scenarios of waste-driven heat loads were investigated to test to what degree, if any, the additional heat affects radionuclide migration by generating convection-driven transport. Results show that the heat output and the calculated radioactivity at a hypothetical near-surface observation point are directly correlated; however, the overall impact of convection-driven transport is small due to the short duration (a few hundred years) of the heat load. Results further showed that the calculated radiation dose at the observation point was very sensitive to the magnitude of the effective diffusion parameter of the host rock. Coupled heat-solute mass transport models are necessary tools to identify influential processes regarding deep borehole disposal of heat-generating long-lived radioactive waste.
1. Introduction
Underground nuclear waste disposal may occur in mined tunnels and caverns at a depth of a few hundred meters [1, 2] or in vertical boreholes at depths up to several kilometres [3, 4]. The host rock in which these geologic disposal facilities (GDF) are being developed may be either a lower strength sedimentary rock such as clay [5] or a higher strength rock such as granite [6]. GDFs have also been developed in rock salt [7]. Deep underground GDFs are the only option considered acceptable for the final disposal of high-level radioactive waste (HLW) [8]. GDFs for HLW may include reprocessed spent nuclear fuel from commercial powerplants or spent fuel elements that are directly disposed without reprocessing. Long-lived intermediate-level waste (ILW) requiring geological disposal includes reprocessed spent nuclear fuel from research reactors and other conditioned wastes from medical isotope production [9, 10]. Reprocessed spent fuel waste packages consist of stainless steel canisters inside which the nuclear waste is immobilised in a borosilicate glass matrix [11]. It is estimated that the global radioactivity confined within borosilicate glass is on the order of 1020 Bq, with an approximate weight of 15,000 metric tonnes [12, 13]. The half-life of long-lived radionuclides in high-level waste is on the order of 105-109 y (e.g., y for 135Cs, for 79Se, and for 238U) [12]. Owing to these long half-lives, HLW requires permanent disposal such that radionuclides remain contained and isolation from humans and the environment for hundreds of thousands of years [14].
Demonstrating the safety of radioactive waste disposal involves postclosure safety assessments using numerical models and developing a synthesis of scientific, technical, administrative, and managerial arguments and evidence known as the safety case [15]. The safety case is built around short-term and long-term experimental observations and numerical modelling of key processes happening within the GDF and the surrounding geosphere and biosphere [12]. The timescale for postclosure safety assessments is typically from ten thousand years to one hundred million years, with one million years a commonly accepted timescale for quantitative assessments [14, 16].
The geological environment for such timeframes can be considered constant, an assumption which facilitates the modelling as no transient features need accounting for. Exceptions are seismic faulting that might interact with the GDF [17]. Other processes that are transient and that could potentially affect radionuclide release and transport through the geosphere, are of anthropogenic origin. One such process is heat generation from the radioactive decay of the waste [18]. Transient heat production occurs for a few hundred years for vitrified ILW [9, 10] or HLW to several thousand years for spent nuclear fuel [19]. Knowledge of the temperature evolution within the engineered barriers and the surrounding host rock is critical to assess if heat-generating waste may impact the safety functions of both the engineered and natural barriers. Indeed, high temperatures will accelerate the corrosion rate of waste forms [20, 21] and metal-based containers [22, 23] and may accelerate the degradation of backfill materials like bentonite [24–27] or cements [28]. High temperatures may also affect the properties of clay host rock [29] and may generate new or propagate existing fractures even in crystalline rock [30]. In addition to having a negative impact on engineered barrier performance, heat production may also generate convective flow and potentially accelerate radionuclide transport through the geosphere [31, 32].
The postclosure assessment strategies typically include several scenarios of projected future behaviour of the disposal facility, including release scenarios for the radionuclide source term [12, 33, 34]. In addition to scenarios describing expected repository evolution, very conservative scenarios or “what if” scenarios are typically tested. An example of the latter could be: “what if the glass matrix is instantaneously and completely dissolved and the stainless-steel disposal container corroded at the time of disposal?”. While purely hypothetical, such “what if” scenarios allow to test—using detailed modelling—the robustness of the disposal system [35, 36].
Following the release of the radionuclides from the degraded waste packages, the transport mechanism into the surrounding geosphere would be a function of the rock petrophysical properties (permeability, connected porosity), the hydraulic gradient, and geochemical conditions that determine radionuclide solubility and sorption [5, 37]. In low permeability rock, mass transport is by molecular diffusion, rather than by advection [38]. To determine whether diffusion is the dominant transport mechanism, the Peclet number (Pe) has commonly been used [39].
For conventional mined GDFs, the majority of postclosure safety assessments have considered isothermal transport of dissolved radionuclides, using simulation codes such as FRAC [40] or PORFLOW [34, 41]. Other studies modelled coupled thermal, hydraulic, and mechanical processes within the near field and host rock of the disposal environment [42–44]. Coupled processes have been studied using TOUGH [45], TOUGHREACT [46], PFLOTRAN [47], and MOOSE [48].
Simulation-based studies involving coupled mass-heat transport processes in geologic environments have been reported across several applications. Harp et al. [49] conducted a heat transport study in a HLW repository in salt with both constant and temperature-dependent thermal conductivity of intact salt. Based on the calculated heat profiles surrounding the waste canisters, recommendations were made about acceptable combinations of canister heat loads and spacings. Using a 100 m deep simulation domain, Jacquey et al. [50] studied coupled flow and heat in a synthetic faulted reservoir for the purpose of code benchmarking (the domain was too shallow to account for depth-dependency of fluid viscosity and density). [51] studied the effect of a constant natural heat flux at the bottom of their 6 km deep domain in combination with heat from decaying waste packages on fluid fluxes and radionuclide concentrations. The estimated travel distance by fluids in the seal zone owing to buoyancy was about 120 m. In a 2-dimensional study, Freeze et al. [52] modelled single-phase, isothermal groundwater flow and nonreactive tracer transport to explore the effects of density stratification and vertical flow and mixing at depth. Considering simulations up to 1,000,000 years with different imposed hydraulic gradients, the initial salinity gradient in the deep brines was shown to be persistent. Ackerer et al. [53] presented a new iterative coupling scheme for solving the transport and flow equation, showing that by solving the transport equation first, the CPU time could be significantly reduced.
To the best of our knowledge, no modelling studies exist for deep borehole disposal which include an evaluation of the effects of different heat sources on radionuclide migration in the rock environment surrounding the borehole, while also testing if these effects are important at different depths of waste disposal [54]. This study therefore includes an explicit linkage between the natural hydrostatic pressure and temperature profiles, density variations, and heat and solute mass transport. We further aim to characterise the level of complexity required for a full-scale postclosure safety assessment to accurately represent all physical processes. This requires identifying the most influential processes and provide some new insights on what processes are first-order controls on radionuclide transport in the subsurface. The aim is to reduce the CPU costs of the solution of these strongly nonlinear coupled equations commonly used in repository safety assessments, without loss of accuracy.
We assess computational affordability and feasibility of a coupled heat and solute mass transport modelling framework, assuming an instantaneous release of radionuclides from the conditioned waste in a hypothetical deep disposal borehole (up to 3 km depth). The model accounts for depth-dependent temperature, pressure, and viscosity profiles for realistic field conditions. Modelling scenarios consider a natural geothermal gradient together with heat-generation from vitrified long-lived ILW, HLW, and spent fuel to test if heat production by the nuclear waste affects radionuclide migration by generating convection-driven mass transport ([31, 55, 56], Sillen, [32]). Diffusive mass transport under these transient thermal conditions is the second process that is evaluated for its importance and sensitivity.
2. Governing Equations
Radionuclide release from a deep disposal borehole into the surrounding host rock will be simulated by considering all relevant migration processes at the Darcy scale while ignoring any mechanical alterations in the borehole itself or in the surrounding host rock (e.g., increased mechanical stresses owing to the heat load [57]). The fundamental solute mass transport equation in porous media is [42]: where is time [], is the concentration [ML-3], is retardation factor accounting for adsorption/desorption, is the dispersivity [L], is the effective advective velocity [LT-1], is the sink/source term [ML-3 T-1] (including contaminated glass dissolution [13], decay of radionuclides, and precipitation/dissolution reactions), is the pore diffusion [L2T-1], is the free water molecular diffusion [L2T-1] and is the tortuosity coefficient. In Equation (1), the advective Darcy velocity vector is: where is the effective porosity, is the groundwater flux [LT-1], is the hydraulic conductivity [LT-1], is the head gradient, is the intrinsic permeability [L2], is the water density [ML-3], is the kinematic viscosity [L2T-1], and is the pressure gradient [ML-2 T-2] (including the body force, here with the buoyancy effect). The heat transfer (thermal energy conservation) equation is: where is temperature [Θ], and are the bulk (fluid and solid) and water density [ML-3], respectively, is the (bulk) thermal conductivity [MLT−3Θ−1], and are bulk and water specific heat [L2Θ−1 T−2], and is again the Darcy velocity. Equations (1)–(3) are coupled partial differential equations that link advection and convection.
For shallow subsurface environments, it is reasonable to assume that the background temperature and density are constant with depth. However, if the depth is on the order of hundreds of meters to several kilometres, the density of water varies with pressure and temperature (which itself changes according to the geothermal gradient), and hence, additional equations that couple these parameters are introduced. Along with the continuity equation, the constitutive relationships between the pressure, density of water, and temperature are required [58, 59].
3. Conceptual Model Setup and Scenarios
We consider coupled heat and water-saturated fluid flow, including buoyancy effects, with temperature-driven variations of the fluid properties (density and viscosity) over the full depth of the simulation domain, and radionuclide decay and retardation. These processes occur in the surroundings of a deep borehole (up to 3 km deep) assuming a fixed natural geothermal gradient and a superimposed transient heat source from the disposal of heat-generating long-lived intermediate-level radioactive waste [9]. Several other processes and variables are not considered, including mechanical alternations of the rocks near the heat source (overheating [60]), geochemical reactions, heterogeneity in the rock formations, and salinity. These are the subject of future studies once the temperature effects on fluid flow and mass transport have been assessed. The modelling platform is based on TOUGHREACT v. 1.2 [46], supplemented by in-house codes to accommodate the boundary conditions.
For the simulation of the above coupled processes involving a borehole with surrounding rock, three possible options for the conceptual site model geometry were initially considered: a full 3-dimensional (3D) Cartesian domain, a 2D Cartesian domain, and a 3D axisymmetric (also 2D vertical) domain. The 3D Cartesian domain allows the modelling of three-dimensional groundwater flow accounting for nonuniform boundary conditions when regional groundwater flows are important (e.g., in shallower fractured rocks), while the geometry of the cylindrical solute and heat sources (i.e., radioactive waste cylinders emplaced in the bottom of the borehole) can be appropriately represented with a fine-spaced grid. Preliminary simulations showed a single-processing core simulation platform was not computationally capable of handling the 3D Cartesian domain with the required spatial resolution to capture depth-wise changes in the density and viscosity, accurately represent the geometry of the disposal borehole, and allow for sufficient numerical accuracy. Further discussion is provided later in this section.
A 2D Cartesian geometry (- plane) enables the consideration of two-dimensional groundwater flow accounting for nonuniform boundary conditions and requires relatively lower computational resources compared to the 3D Cartesian geometry. However, approximating the cylindrical waste canister by means of a point solute and heat source leads to inaccuracies, as the fundamental assumption is the symmetricity (or no-flow as in TOUGH group of codes) over the third dimension (-direction orthogonal to the - plane). This results in an overestimated release of radionuclides and heat [61].
A 3D axisymmetric geometry allows accurate representation of a cylindrical source, as the transport phenomena is largely radial. A disadvantage is that groundwater flow driven by external boundary conditions (e.g., strong regional hydraulic gradients) cannot be modelled due to the symmetricity around the -axis. Applications of axisymmetric models have mostly been restricted to groundwater flow towards an extraction well or away from an injection well, although coupled flow and transport examples such as dense saline groundwater in a partially penetrating well have been reported [62].
To apply a 3D axisymmetric geometry to our problem, a justification is needed to assess under what conditions the regional groundwater flow becomes immaterial to heat and mass transport and whether such conditions exist for our test case. Therefore, a basic dimensional analysis was undertaken using the length-scales for molecular diffusion (), heat conduction (), and solute advection () [63–65]: where is (effective) molecular diffusion [L2T-1], is the thermal diffusivity [L2T-1], is pore-water velocity (LT-1), and is time [T]. For deep crystalline rock, we assume (assumed to be constant for the sake of simplicity in the analysis, despite effective diffusion may vary with depth, see further), and . We further consider Equation (2) with for crystalline rock [51], , and the pressure gradient is assumed to represent a regional hydraulic gradient of 0.0025 m/m [9]. For the first 106 years, the advection length scale () is smaller than the heat () and diffusion () length scales, meaning that solute advection is not the dominant transport mechanism under the above conditions (see Figure 1). This provides a simple justification for using an axisymmetric geometry without groundwater-driven solute transport, as the assumed regional hydraulic gradient in combination with the very low rock permeability produces a negligible velocity (). As the axisymmetric model is a suitable choice for the simplified system considered here, it was preferred over the 2D and 3D cartesian models.
Based on an axisymmetric geometry, a three-layer conceptual model is adopted that includes the following stratigraphic layers: a 20 m thick regolith, underlain by an 80 m thick weathered basement rock and a 3100 m thick crystalline rock (Figure 2, left). For the crystalline rock, the upper 500 m ( to ) is assumed to have an effective diffusion coefficient of 10-10 m2/s (Shapiro, [66–68]). The rest of this layer has an effective diffusion coefficient of 10-11 m2/s. Sensitivity of radionuclide transport with respect to the diffusion coefficient of the upper 500 m will be tested using a five times larger value, i.e., .
All layers are assumed to be homogeneous and isotropic in the radial and vertical direction, with physical properties reported in Table 1. Plutonic and volcanic rocks are isotropic regarding thermal conductivity, while anisotropy is often observed for sedimentary and metamorphic rocks [69, 70]. This simplification was primarily done to minimise computational cost. Two different permeability values (10-16 and 10-17 m2) are considered for the crystalline rock, allowing testing of model sensitivity relative to permeability.
The axisymmetric domain has a radius of 3200 m and a depth of 3200 m. The side and top boundaries are of Dirichlet type (for pressure, temperature, and solute mass). A background geothermal gradient of 20°C/km is assumed, which is representative of low-heat producing areas in Central Australia [71, 72]. The simulation results (temperature and radionuclide concentration) are collected along the radial direction at the centre of the source (i.e., at the depths of 400, 900, and 2900 m). Also, at the depth of , a hypothetical domestic groundwater well is assumed (Figure 2, right). The location of the well in the radial direction is made variable between simulation cases and is put at the location with the highest concentrations; in this way, the most conservative results in terms of impact are obtained. No pumping of the well is applied; this again is conservative as pumping typically results in additional dilution with radionuclide-free groundwater.
To test the sensitivity of radionuclide transport and their potential radiological impact on humans via a groundwater well relative to the disposal depth, three borehole disposal depths are evaluated: 0.5 km, 1 km, and 3 km. The radionuclide and heat source are located at the bottom 200 meters of the borehole; the waste disposal zone is assumed to accommodate a total of 100 waste canisters, each about 1.3 m long. For the purpose of simulation, the borehole is assumed to have a bottom-hole diameter of 0.6 m, consistent with current drilling technologies [73]. The remainder of the borehole from the top of the waste disposal zone until the ground surface is assumed to be backfilled with a single material to simplify the analysis. In the current model, a cementitious backfill is used, although other materials such as bentonite and crushed rock have also been considered [51].
One of the aims of this study is to test the sensitivity of radionuclide transport with respect to the strength of the heat produced by the decaying radioactive waste. To this end, the heat source is represented by three levels of heat outputs at the time of disposal. These are 50 Watt/canister, 500 Watt/canister, and 1500 Watt/canister [9]. The lowest heat output (50 Watt/canister) is typical of vitrified waste from reprocessed research reactor fuel, 500 Watt/canister is typical of waste from reprocessing of fuel from commercial power plants (HLW), while the highest output (1500 Watt/canister) represents spent fuel from commercial nuclear power plants. These three cases use a simplified heat production function with time (Figure 3): where is the initial heat output (Watt/canister) at the time of disposal. The time evolution of thermal load was approximated by assuming the heat source decays according to the half-life of 90Sr (29 years), which together with 137Cs (30 years half-life) are important fission products of nuclear reactions and major contributors to decay heat [74].
Exact determination of extended heat output profiles is possible through, e.g., analytical approximations for the long-term decay behaviour of spent fuel and HLW [75]. For the shallowest disposal hole (0.5 km depth), an extended heat output (EHO) profile is considered to assess the effect of prolonged heat generation, which represents a more realistic heat production including the decay of long-lived radionuclides (e.g., by actinides) present in spent fuel from commercial nuclear power plants (Johnson et al. [19]) (Figure 3).
The heat output for these four cases has also been displayed per meter of waste in the disposal zone (Figure 3(b)). This way of presenting the heat output will aid with the interpretation of the results, as the length of the canister for the fourth case (EHO) is larger (i.e., 5 m) than the length of the other three canisters (assumed 2 m, which includes the total length of the canister (1.3 m) plus a short buffer zone in between two neighbouring canisters).
While long-lived radioactive waste from nuclear fuel contains a large number of radionuclides [19], only a selected set will be discussed here for the purpose of illustrating how coupled heat-transport processes may affect long-term migration. Selected radionuclides are 137Cs, 239Pu, 79Se, 99Tc, and 126Sn, with their initial activity concentration and transport properties summarized in Table 2. While the initial radioactivity is given in Becquerel (Bq), TOUGHREACT requires molar concentration as the refence; hence, the following conversion is used: where is number of moles, is the Avogadro constant, and is the half-life (s). Sorption of radionuclides to the solid phase of the different materials is described by means of the linear solid/liquid partition coefficient or distribution coefficient (Table 2).
No engineered barriers were assumed other than the cement-base seals across the 200 m long disposal zone. This is a reasonable assumption because the waste containment safety function will be primarily provided by the tight host rock, rather than by engineered barriers [9]. In other words, rather than a time-dependent corrosion of the stainless-steel disposal container and a time-dependent release of radionuclides from the vitrified waste form, an instantaneous release is assumed whereby the entire initial activity is completely dissolved and available for transport.
Initial model runs were undertaken to optimize the model grid for the 3200 m deep and 3200 m radial dimension of the simulation domain. In the absence of any regional hydraulic gradient for this axisymmetric model, any advective transport will be of very low magnitude and driven only by density variations due to changes in pressure and temperature. Initial simulations showed that a too coarse mesh would result in an inaccurate calculation of water density at the boundary cells, and therefore, an erroneous advective velocity would be produced in both vertical and radial directions (results not shown). By using an even distribution of 3.3 m long grid cells in the -direction across the entire depth, acceptable results were obtained with advective velocity excluding heat effects (in -direction) not exceeding 10-14 m/s. Also, mesh clustering in the radial direction was applied, with spatial increment defined as , where (varying from 1 to 91) is the cell index in direction. Three distance views of the structured mesh are depicted in Figure 4. Total number of grid cells is 88,320. A further mesh sensitivity analysis for selected cases showed that the mesh quality is acceptable (with respect to the convective flow pattern and a passive scalar transport; results not shown). The structured mesh could be replaced by an unstructured mesh to better allocate the grid cells over the volume, although no attempt was made here.
Changes in the density of water is crucial to properly calculate the flow field. This is most pronounced in cases where changes in the concentration of dissolved solids is significant (i.e., brine). Note that brine was excluded from the current calculations to simplify the investigations of heat effect on flow and transport. Furthermore, the single-processing version of the TOUGHREACT code applied here does not include the EOS7 module, which has the equations of state for brine and water transport. The TOUGH3 parallel-processing group of codes have limited capabilities in this regard (the EOS7r module) [76]. Nevertheless, options such as PFLOTRAN [77] or the parallel processing version of the TOUGHRACT-v3.3-OMP, TOUGHREACT-V4.0-OMP, and TOUGHREACT-Brine seem to be better candidates if a proper 3D study with heterogeneities and groundwater flow is required [78].
Fluid viscosity of sedimentary basins is mainly dependent on temperature (decreases with increasing temperature), less dependent on salinity, and even less dependent on pressure (increases with increasing salinity and pressure). According to Adams and Bachu [79], the effect of pressure ranges typically observed in sedimentary basins causes small increases in brine viscosity (<5%). In the current model, effects of both temperature and pressure on fluid viscosity are accounted for based on the International Water Property Standards 2009 Formulation for the Viscosity of Water [80]. These dependencies were also implemented to generate appropriate boundary conditions for modelling.
Assigning representative pressure values at the side boundaries that reflect the conceptual hydrogeological model is of critical importance. To this end, an in-house code, linked to the external MATLAB code XStem [58], was developed to assign the side boundary conditions. The procedure is as follows: (1)Assume initial hydrostatic pressure distribution with constant water density(2)Update density as a function of pressure and background temperature(3)Update pressure with the updated density(4)Step 2 is repeated until the results converge
Allowing the viscosity (and density) to depend on temperature across the 3000 m deep depth profile resulted in a curvi-linear decrease from at the top of the model to at 3000 m depth. The temperature increase over this distance is 60°C. The effect of temperature on the viscosity of water is included in Figure 5. Figure 5 further shows the converged values for the density and pressure over the full depth of the model domain, the difference between the hydrostatic (constant density) pressure and the converged pressure values, and finally, the depth-dependency of the hydraulic conductivity (according to Equation (2)). The interdependencies between hydrodynamic parameters are provided in Table 3.
4. Results and Discussion
4.1. Temperature Profiles
The coupled heat-flow-mass transport model was run for a total of 23 simulation cases by varying heat output (3 cases: 50, 500, and 1500 W), crystalline rock permeability (2 cases: and ) and disposal depth (3 cases: 0.5, 1.0, and 3.0 km). The nineteenth and twentieth cases involved the extended heat output (EHO) assumption to test if simplifications in the heat-decay model for the other three heat loads were acceptable; the EHO was only applied for the 0.5 km disposal depth (for the two values of permeabilities). The sensitivity of the model output to the effective diffusion coefficient was tested by increasing the reference value of 10-10 m2/s five times in the upper 500 meters of the crystalline rock with and for the 50 W, 1500 W, and EHO cases (21-23).
The run time for each case was approximately 8 hours for a simulation time of – years (depending on the case). Results include temporal and spatial concentrations (TOUGHREACT uses mole/L as concentration unit which was converted to equivalent radioactivity concentration (in Becquerel/m3 according to Equation (9)) for purposes of display) of all the radionuclides, pressure, temperature, and advection velocity (as a result of fluid density variations caused by the heat-generating waste).
Temperature evolution for heat generating wastes (50, 500, and 1500 Watt/canister) is shown at the middle of the 200 m long source zone (at depths of -400, -900, and -2900 m for the 500, 1000, and 3000 m deep borehole, respectively) at the interface between waste zone and rock (at radial distance ) and within the rock (at and 5 m) (Figure 6, rows 1-3). For all cases, the temperature transient does not extend beyond year 200 as a result of the assumed decay heat evolution (see Equations (7) and (8)) which only considers the effect of short-lived radionuclides. Figure 6 (row 4, third column) also shows the temperature evolution for the EHO case at and and 5 m. As expected from the heat output evolution for the EHO case (Figure 3), the effect of heat output on the local temperature increase lasts beyond 10,000 years.
The maximum recorded temperature for the 50-Watt simulation case occurs around year 10 and at the interface of waste zone/rock (Figure 6). The peak temperature is 68, 35, and 26°C for the disposal depths of 3 km (observation point at ), 1 km (observation point at ), and 0.5 km (observation point at ), respectively. The temperature at these three depths is the sum of the temperature increase from the waste (approximately 5°C, independent of depth) and the background temperature using a geothermal gradient of 20°C/km: 21°C at 0.4 km depth, 29°C at 0.9 km depth, and 63°C at 2.9 km depth. At the radial distance of 5 m (i.e., 4.7 m into the crystalline rock), peak temperatures occur at year 30 and have dropped to 23, 32, and 66°C at 0.4, 0.9, and 2.9 km depth, respectively. This is only a few degrees above the background temperature, meaning that the temperature effect on the rock for a low heat output of 50 Watt is also limited in space, in addition to being limited in time (up to 200 years).
Peak temperatures for the 500-Watt heat output case at the waste zone/rock interface are 75, 85, and 118°C for the 0.4, 0.9, and 2.9 km depth, respectively. About 20 years later peak temperatures have reached 4.7 m into the crystalline rock and approach 45, 50, and 85°C for the 0.4, 0.9, and 2.9 km depth, respectively. Maximum temperature increase in the host rock at these depths is thus 24, 21, and 22°C, which is more than an order of magnitude higher than for the 50-Watt case.
For the 1500 W simulation case, the overall maximum recorded temperature occurs around year 10 at the interface of waste zone/rock (Figure 6). The peak temperate is 230, 200, and 190°C for the deepest, medium depth, and shallowest disposal depth, respectively. At the radial distance of 5 m, peak temperatures at year 30 have dropped to 90, 100, and 120°C at 0.4, 0.9, and 2.9 km depth, respectively. This is still 69, 71, and 57°C above the background temperature. The effect of the heat production becomes negligible after around 200 years, when the maximum temperature has become only 10% above the background temperature. In other words, the transient temperature phase during which buoyancy could play a role is limited to less than 200 years for the 1500 Watt/canister heat load. As can be expected, the transient phase for the EHO case lasted much longer (at least for 10,000 years, see Figure 6).
The radial evolution of temperature versus time (1500 Watt/canister, observation depth at ) shows that the temperature increase in the rock mass does not extend beyond approximately (Figure 6; bottom left). Similar observations were recorded for other cases (not shown) and illustrate that a relatively small volume of rock would experience elevated temperatures.
The highest recorded temperatures and temperature increases relative to the background temperature at a radial distance for different heat outputs are summarized in Table 4. This additional observation point was chosen as it would also serve for outputting of radionuclide concentrations (see further). These maxima are observed at 10-15 years, irrespective of depth.
When temperatures exceed 100°C at pore-water pressures larger than atmospheric, boiling conditions are not necessarily met or exceeded. Typically, boiling temperatures increase with increasing pressure. At 500 m depth, boiling temperature is about 260°C, while at 3000 m, it is about 400 °C. The maximum temperature recorded in Table 4 was 211°C for 1500-Watt canisters. In other words, this is still far below the boiling temperature.
Nonetheless, the engineered barriers and rocks which surround waste packages may change their properties due to heating. For example, at temperatures above 100°C, bentonite clay may show mineralogical alterations [26, 81, 82] that may lead to a less performing barrier. Also, at temperatures above the boiling point of water, phase changes occur which would invoke several coupled multiphase processes that are difficult to characterize and simulate. Such conditions should therefore be avoided [83].
Under some extreme temperature conditions, rock can even melt. This can occur when the concept of rock melting is applied [11, 84], where the heat production from the waste is sufficiently large to melt the rock (for granite approximately 950°C [85]). The corresponding required thermal power for melting is approximately 16 kWatt/m3. This is much larger than heat productions considered here.
Effective diffusion coefficients, thermal conductivity, and heat capacity are known to be temperature dependent. Experimentally determined effective diffusion coefficients at 80°C for Callovo–Oxfordian claystones showed an increase by a factor of 3 for tritiated water and a factor of 2 for caesium compared to values measured at 21°C [86]. This is consistent with an increase by a factor of 3.2 at 85°C compared to 25°C, assuming diffusion can be scaled using the temperature/viscosity ratio [87]. For Beishan granite, effective diffusion coefficients for TcO4- increased from at 25°C to at 55°C [88]. Temperature dependence of thermal conductivity and heat capacity was experimentally demonstrated for crystalline rocks by Vosteen and Schellschmidt [70], with the former decreasing with increasing temperature while the latter increases as temperature goes up. From 50 to 200°C, conductivity decreased from 2.3 to 2.2 Watt/m K, while the heat capacity increased over this range from 800 to 925 (J/kg K).
Such temperature effects on petrophysical properties (e.g., diffusion coefficient, thermal conductivity, and heat capacity) of the rock were not considered here. For the current temperature range (Figure 6), variations in thermal rock properties (thermal conductivity and heat capacity) are considered to have a small effect on the temperature evolution, while temperature dependency of the diffusion coefficient is considered unimportant to radionuclide transport given the limited radial extent of the rock zone experiencing any significant temperature increase (Figure 6 bottom row). Also, in the current version of TOUGHREACT, thermal conductivity, heat capacity, and diffusion coefficients for the nongaseous phase are constants.
4.2. Convection-Driven Flowlines
The early-stage patterns of convection-driven flowlines and radionuclide activity concentrations in the vicinity of the combined heat and radionuclide source for the 0.5 km disposal depth are displayed in Figure 7 (streamlines and velocity magnitude). Figure 8 also shows the early-stage temporal variation of the water density at . In Figure 7, flow streamlines at and 1000 years are shown for an area between and and for depths between and . At , the greater heat output (1500 Watt) results in a greater area of the rock being impacted by the buoyancy effect as evidenced by the (1) convection streamlines extending further in the radial direction compared to the 50-Watt case and (2) the greater vertical distance over which the velocity increases. For instance, the furthest streamline for the 50-Watt case reaches as far as and then travels further upwards. At 1500 Watt, the furthest streamline goes as far as (). Comparing velocity magnitudes at for 50 and 1500 Watt shows values of ~10-13 m/s at for the former and 10-12 m/s for the latter. In other words, the buoyancy effect has not reached for 50 Watt whereas for 1500 Watt, it is still large enough to increase the velocity an order of magnitude above its background value (due to the background temperature profile, see background value of 10-13 m/s at years when temperature has returned to its background value). The effect of higher permeability on the area of host rock with a noticeable convective velocity is most pronounced for the greatest heat output. At 1500 Watt, the streamlines travel beyond for (not shown in the scale provided), while they reach only up to R = 35 m for k =10-17 m2. The buoyancy effect (most pronounced for the 1500-Watt case, see also Figure 8) will drive more radionuclides upwards per unit of time compared to diffusive transport, effectively causing higher concentrations at observation points above the waste zone.
The second column in Figure 7 shows the convective velocity magnitude () in an area above the waste zone to in year 50. At 1500 Watt, the magnitude of convection (note the changes in the colour scale) reaches a maximum of and (for m2 and 10-17 m2, respectively), at least one order of magnitude higher than at 50 Watt. This larger convective velocity, mainly oriented in the upward direction, will drive solute transport more in the upward direction rather than in the horizontal direction. In other words, more solutes will migrate upwards in the 1500-Watt case compared to the 50-Watt case (Figure 7).
Columns 3 and 4 of Figure 7 show the flow streamlines and the velocity magnitude after 1000 years. Clearly the effect of heat production has vanished, and the system returns to its natural state with convective velocities of about 10-12 m/s (for ) and 10-13 m/s (for ). Importantly, the buoyancy effect only lasts for several hundred years and has returned to its natural state within 1000 years, as also reported elsewhere [51, 89]. The simulations show that the maximum calculated convective velocity magnitude is (~5 mm/year), indicating that for long time frames, advection (compared to diffusion) is negligible for the assumed disposal depth.
The velocity magnitude for the 50-Watt cases reaches the background velocity at or around the elevation of . For the 1500-Watt cases, the velocity magnitude at year 50 remains 2 times above the background velocity at (not shown in Figure 7).
For the boreholes with disposal zones at deeper depths (2000 and 3000 m), results would be similar as similar zones of the host rock would experience similar changes in velocity for similar durations. These depths are therefore not further discussed.
4.3. Radionuclide Concentrations
Figure 9 shows the total radioactivity (the sum of the activities of all five radionuclides) at a radial distance from the centre of the source, or 0.7 m into the rock, for three heat outputs (50, 500, and 1500 Watt/canister). The maximum activity concentrations occur around year 40, regardless of the heat output and disposal depth. This represents a delayed transport of solute mass compared to the heat pulse from the waste, as was also observed in the dimensional analysis in Figure 1 (also see Figure 6, the temperature profiles at , where the maximum temperature occurs between 10 and 15 years). The contribution of the five radionuclides in the peak values (year 40) are approximately 1, 2, 41, 48, and 8 percent for 137Cs, 239Pu, 79Se, 99Tc, and 126Sn, respectively. The activity concentrations half-way the disposal zone (-direction) and for for all heat sources increase slightly with increasing depth of the disposal boreholes. This is due to accelerated upward radionuclide transport based on buoyancy-driven flow (see discussion in Figures 7–8). As higher temperatures decrease the viscosity, deeper boreholes and higher heat outputs will both contribute to decreasing viscosity (see Figure 5: the viscosity at is almost 2.5 times smaller than the counterpart value at ; versus 10-3 Pa.s). A decrease in viscosity will increase the hydraulic conductivity, which will further increase buoyancy-driven flow.
Figure 10 shows the activity concentrations of 79Se and 99Tc (at , or 50 m above the top of the radionuclide source zone) in the radial direction at early time following the release of radionuclides, i.e., at years. Among the five radionuclides selected for demonstration purpose, 79Se and 99Tc have among the longest half-life and are the least sorbed on the crystalline rock (Table 2), therefore are of greater interest from an impact assessment point of view (the other three had negligible (physically meaningless) concentration values). For both radionuclides, higher concentrations exist for the 1500 Watt compared to the 50-Watt case. Of further note is that the peak concentration for both heat load cases occur almost at the same distance (despite a slight shift to the right in the 1500 W cases); therefore, it is inferred that the shape of the radionuclide plume in the radial direction is somehow independent of the heat load. This seems in contradiction to the convective field at 1500 Watt being spread out over a greater radial distance (Figure 7, first column). However, if only the first 5 m are considered, then flowlines between 50 and 1500 Watt are similar (Figure 7). As a result, concentration plumes and breakthrough curves (along the -axis) are nearly identical, except for the magnitude of the radionuclide concentrations. Higher concentrations for the 1500-Watt case are the result of the greater upward convective velocity causing more radionuclides per unit of time to be vertically transported and accumulated at the observation depth ().
The effect of a 10 times higher permeability (10-16 compared to 10-17 m2) on activity concentrations (for the same heat production and radionuclide) is observed for 79Se and 99Tc, with peak concentrations increasing 1.5 and 5.9 times for 79Se at 10-16 m2 (for 50 and 1500 W, respectively) and 2.5 and 6.1 times for 99Tc at 10-16 m2 (for 50 and 1500 W, respectively). The higher convective velocity at will transfer more radionuclides per unit of time from the waste zone upwards resulting in higher concentrations. A simple dimensional analysis is helpful here. Considering Equations (4) and (6) (and approximating the effective diffusion coefficient equal to the free-water molecular diffusion coefficient times the porosity) for years, the diffusive length scale is 0.8 m. The convective velocity equivalent to this length scale over 500 years is then . With the velocity profiles given in Figure 7 (rapidly depleting after year 50), wherever the velocity magnitude is greater than , the effect of advective transport is dominant over diffusive transport. Effects of radioactive decay are negligible at very early stage (i.e., at 500 years, except for Cs-137); therefore, differences in transport velocity have little effect on a decrease in concentrations due to decay.
4.4. Dose Rates
Next, we discuss the maximum radionuclide activity concentrations and their corresponding maximum annual dose at a hypothetical domestic well which is screened at the depth of . The location in the radial direction is not fixed but is adjustable such that the well will always return the maximum activity concentration depending on where in space the maximum concentration occurs. The space-time evolution of the radionuclide plumes is a function of the coupled processes, i.e., buoyancy driven flow as the result of thermal gradient and heat generation (Equations (1)–(3)) and transport by molecular diffusion. The concentrations are further influenced by the rock properties (here adsorption and permeability) and radionuclide half-life. The magnitude of convective transport is driven by the buoyancy term (the body force in the momentum conservation equation) and proportional to the intrinsic permeability and inversely proportional to the kinematic viscosity (Equation (2)). Two different values of intrinsic permeability were considered for the shallowest borehole (0.5 km depth): and . For the 1000 and 3000 m deep disposal borehole, results are shown for the most permeable case (). The shallowest borehole was selected for the dose calculations as it has the shortest travel distance between source and receptor (screened well) and would produce the overall highest dose among the three disposal depths. This risk-based approach first considers the upper bound impact before proceeding, if required, with lower impacts.
Among the five radionuclides tested, 79Se and 99Tc returned the highest concentrations and corresponding doses; the discussion will therefore proceed based on these two radionuclides. The contributions from the other three radionuclides to the total dose will be discussed later. Activity concentrations and corresponding annual doses for 79Se and 99Tc are discussed based on longitudinal profiles from to at the depth of the water well () based on the shallowest disposal borehole (bottom hole at 0.5 km) (Figures 11 and 12). The annual dose was calculated as the activity concentration (Bq/m3) times the annual amount of drinking water intake per person (0.73 m3/year) [90] times the dose coefficient for ingestion ( Sv/Bq for 79Se and Sv/Bq for 99Tc) [91].
Inspection of the longitudinal profiles shows that for a given combination of heat output and permeability, the spatial locations of the centre of mass are not fixed (Figures 11 and 12). Indeed, at the time of the first profile (), the centre of mass at the observation point for both radionuclides initially occurs at approximately 40 m (for 79Se, Figure 11) and 60 m (for 99Tc, Figure 12) from the centre of the disposal borehole. From this, we conclude that at times up to , the borehole itself is not the main transport pathway for the conditions applied in these simulations. This is due to several reasons, including buoyancy taking place over a considerable fraction of the rock mass leading to smearing out of the radionuclide plume in radial direction (Figure 11), and the higher adsorption capacity of the borehole shaft (10 L/kg for selenium and 700 L/kg for technetium) relative to that of the rock matrix (0.2 L/kg for selenium and 0.1 L/kg for technetium, Table 2). The retardation factor influences both transport by diffusion and by advection, although the effect on the advection will be limited in time (for as long as the buoyancy effect exists, which is less than 1000 years for all but the EHO scenario) while the effect on diffusion is for the full duration of the simulation. For example, if the convective velocity was identical within the borehole and the surrounding rock, then the retarded velocity —with the retardation factor—in the borehole would be () compared to () in the crystalline rock for Selenium. For Technetium the difference would be even larger, i.e., () compared to () in the crystalline rock. Note that Figure 7 (second column) shows the highest convective velocities in a region including the disposal borehole; however, these are unretarded velocities and do not account for the radionuclide-specific retardation. The latter effect is clearly visible in Figure 11, where the peak concentration has moved away from the borehole owing to reduced sorption to the rock matrix (note that only 79Se has a in the borehole—at radial distance —as differences in retardation factor between rock and backfill are less than a factor 2; Tc consistently has near zero concentrations in the borehole, see Figure 12).
As time proceeds, 79Se and 99Tc start to behave differently. For 79Se, the centre of mass gradually moves towards the centre of the domain (i.e., the disposal borehole), which is especially visible at late time () (Figure 11). For 99Tc, the centre of mass remains more or less at the same distance from the borehole centre, i.e., at 60 m (Figure 12). Another important observation is that concentrations at the position of the borehole () gradually increase. This is true for both radionuclides, although more apparent for 79Se (Figure 11) than for 99Tc (Figure 12). This behaviour is due to lateral diffusion of the radionuclide plume back into the initially radionuclide-free disposal shaft as the plume appears in the upper part of the model. Any contribution to the increasing radionuclide concentrations in the disposal shaft from upward radionuclide movement within the shaft itself is negligible owing to the strong sorption to the cement backfill ( for 79Se and 700 L/kg for 99Tc) which retards the upward transport relative to the transport in the horst rock ( for 79Se and 0.2 L/kg for 99Tc). A second effect of the high sorption to cement is that pore-water concentrations will be much lower in the cement relative to the surrounding rock, as more radionuclides will be partitioned to the cementitious solid phase relative to the crystalline rock (Figure 13).
Until now, an effective diffusion coefficient was used for the upper 500 m of the crystalline rock. Across many studies, the diffusion coefficient has been observed to be correlated with rock porosity and therefore a function of depth (as porosity typically decreases with depth). Lower porosities usually exist at greater depth due to the greater consolidation [36, 92]. Studies for crystalline rock have assumed an increase in effective diffusion of one order of magnitude when moving from deep to shallow rock [87, 93]. For different argillaceous rocks including both indurated and plastic clay, was found to range between 10-11 and 10-12 m2/s [36]. Therefore, postclosure safety assessments for geological disposal have often applied an uncertainty factor of 10 in their radionuclide transport calculations [36, 94]. Therefore, when diffusion is the main transport mechanism over long time scales, uncertainties about the diffusion parameter need to be accounted for. To this end, additional simulation cases were considered based on a greater effective diffusion in the upper 500 m of the crystalline rock ().
Calculated annual total dose rates summed over all radionuclides are shown for three heat loads (50 and 1500 Watt, and EHO) and two diffusion parameters (Figure 14). The main contributors to the annual dose rates were 79Se and 99Tc (their combined contribution is 99% of total dose). For the smaller diffusion coefficient (10-10 m2/s), the radionuclide concentrations and corresponding annual dose rates (on the order of 10-6 mSv/y) are negligible and many orders of magnitude smaller than what the IAEA [95] considers an insignificant dose to humans (<0.01 mSv/y). Of further importance is the observation that the peak dose is rather insensitive to the heat load. Both the magnitude of the heat load and the duration of the heat production (short-term as per the 50 and 1500-Watt assumption or long-term as per the EHO data) have little influence on the dose and produce a nearly identical peak dose. This further shows that the initial assumption of simplifying the heat production based on decay of a short-lived radionuclides is valid, at least for the conditions of this study.
Use of a larger diffusion coefficient results in significantly larger dose values (Figure 14), with the maximum dose around . Such a peak dose is still considered very low, illustrating that the level of safety for the considered combination of borehole depth (500 m) and geologic conditions is very high. The very large sensitivity of the peak dose with respect to the magnitude of the effective diffusion coefficient requires that future rock characterization studies associated with deep geological disposal undertake appropriate testing to derive site-specific diffusion values and to test for any depth-dependency. A combination of laboratory tests on undisturbed cores [96], in situ testing in an underground research facility [97], or borehole research facility together with interpretation of large-scale environmental tracer data such as helium profiles [98] will likely provide the most representative data.
The annual insignificant dose to humans (<0.01 mSv/y) is then used as a threshold to visualize from what depth onwards (or travel distance from the source) the dose indeed becomes insignificant. Figure 15(a) shows that the total peak doses for a 500 m deep borehole have decreased to 0.01 mSv/y at depths between 100 and 150 m in the time interval 0.5- years. For times years, the 0.01 mSv/y threshold drops to greater depths owing to radioactive decay now becoming more effective. The effect of different assumptions about heat load is noticeable but relatively limited: at 106 years, the 1500-Watt case has its threshold of 0.01 mSv/y about 30 m higher in the profile compared to the EHO case, which itself is about 10 m higher than the 50-Watt case. Importantly, for all cases, annual dose rates remain very small with doses that would be observed at the hypothetical well even smaller (but not calculated because of insignificant values).
For the 800-1000 m disposal depth (Figure 15(b)), a somewhat similar behaviour of the peak dose is observed, although there are differences too owing to the use of a lower effective diffusion coefficient (10-11 m2/s) in the deeper rock zone (>500 m depth) compared to the shallower rock zone. The 0.01 mSv/y threshold is now reached after a travel distance of approximately 120 m above the top of the waste zone, after a travel time of approximately .
The travel distance at which the 0.01 mSv/y threshold is reached for the deepest borehole (waste emplacement between -2800 and -3000 m, Figure 15(c)) is almost identical to that for the previous depth, i.e., nearly 80 m above the top of the waste zone. Effects of different heat loads are again negligible, with the threshold being reached practically at the same depths and time. These shorter travel distances are due to the smaller effective diffusion coefficient in the deeper boreholes, compared to the 500-m deep borehole.
Figure 16 depicts from what depth onwards the dose for all the five radionuclides becomes insignificant (less than 0.01 mSv/y). As the short-lived 137Cs decays much more rapidly than the other radionuclides and has the overall highest adsorption (), it reaches the 0.01 mSv/y threshold after a much shorter travel distance (about 10 m and 200 years). The somewhat longer-lived 239Pu (24,100 y) reaches its threshold also at about 10 m from the top of the source, but remains at that level for a longer time than 137Cs due to its longer half-life and relatively strong adsorption ().
The behaviour of the 126Sn threshold is somewhat different: it increases very slowly during the first 105 y, the result of its relatively long half-life (230,000 y) in combination with adsorption coefficient .
This study evaluated what processes are first-order controls on radionuclide transport in the subsurface with the aim to reduce the CPU costs of the solution of these strongly nonlinear coupled equations commonly used in repository safety assessments, without loss of accuracy. The coupling of heat transport with radionuclide migration which involved accounting for buoyancy-driven transport was shown to have a limited importance, at least for the assumptions of this study. These findings suggest that simpler models can be used without loss of accuracy. We note that the importance of accounting for additional heat-mass transport dependencies is function of many factors; therefore, there is no single answer. Here, the final goal was to quantify radionuclide concentration and radiological dose to humans from potential exposure to groundwater containing several radionuclides. As was shown for different heat production scenarios, while also varying the rock permeability, the total dose from considering five radionuclides and different heat sources was virtually insensitive to the buoyancy process. From this result follows that accounting for temperature-dependent fluid viscosities was relatively unimportant. This conclusion should not be generalised, as it depended on many factors, such as negligible advective transport owing to the small permeability and hydraulic gradient, consideration of a limited number of radionuclides, a disposal zone between 300 and 500 m depth, and a rather low geothermal gradient of 20°C/km.
5. Conclusions
Coupled flow, heat, and radionuclide transport calculations were undertaken to assess to what degree heat produced from the decay of radioactive waste leads to buoyancy-driven migration of radionuclides from a deep disposal borehole. Radionuclide migration was uniquely by molecular diffusion in the absence of a hydraulic gradient; a convective component was added due to a transient heat source associated with the decaying radioactive waste. Of the four heat sources compared, three had a short-duration heat output for about 200 years (50, 500, 1500 Watt/canister), while the fourth had a long-term heat output for more than 10,000 year (2200 Watt/spent fuel assembly).
Regardless of the depth of the disposal borehole, a rock volume with a maximum radius of approximately 60 m surrounding the disposal borehole is influenced by the transient heat production for the 1500 Watt/canister with buoyancy creating upward convective velocities one order of magnitude larger than the velocities at 50 Watt/canister heat source. The transient nature of this buoyancy effect is limited to approximately 200 years, after which the heat source is exhausted. During this transient, radionuclide migration is affected by the temporarily increased upward velocities which causes slightly higher concentrations at observation points above the disposal zone for the highest heat source.
Even for the shallowest borehole (500 m), the radionuclide concentrations and annual dose rates are negligible, and many orders of magnitude smaller than what the IAEA [95] considers an insignificant dose to humans (0.01 mSv/y). This threshold total dose was reached after a travel distance between 150 and 200 m from the source. Increasing the diffusion coefficient by a factor 5 resulted in a much larger dose rate, but still low enough to be safe. This result showed the sensitivity of the dose rate to the diffusion parameter under the conditions of the model, and the importance to accurately measure influential parameters like the effective diffusion for site-specific applications. As a unique finding, the coupling of heat transport with radionuclide migration which involved accounting for buoyancy-driven transport was shown to have limited importance, at least for the assumptions of this study.
For the two deeper boreholes (1000 and 3000 m), the required travel distance to reach the 0.01 mSv/y total dose threshold was similar (about 80 m from the top of the waste zone, with 99% contribution from 99Tc and 79Se), but somewhat less than the shallower borehole owing to the lower diffusion assumed at greater depth (five times smaller). The shorter-lived 137Cs and 239Pu reached their 0.01 mSv/y threshold after 10 m travel distance, which underscores the containment capacity of the assumed host rock at 1000 and 3000 m depth. For 126Sn, the 0.01 mSv/y threshold was reached after 20 m travel distance, in between that for 137Cs/239Pu and 99Tc/79Se.
These results are preliminary in that they are based on very conservative assumptions about radionuclide release. Future work will include use of more realistic release models and the testing of scenarios that include seismic faults and poor sealing of the borehole as potential conduits for radionuclide transport.
Data Availability
PetraSim (a commercial software by Thunderhead Engineering Consultants) was used for preprocessing. It included TOUGHREACT, which is a member of TOUGH group of commercial codes, delivered by Lawrence Berkley National Laboratories (United States). The in-house pre- and postprocessing codes used in this study were developed in MATLAB, which is a commercial software delivered by MathWorks. Theses codes may be obtained by contacting Kaveh Sookhak Lari (Kaveh.sookhaklari@csiro) or Dirk Mallants (dirk.mallants@csiro).
Conflicts of Interest
The authors declare no conflict of interests for this submission.
Acknowledgments
This study is a part of the Safe Waste Disposal Project (SP-IE-Safe Waste disposal FY21/22) which was funded by the Land and Water Business Unit of CSIRO, Australia’s National Science Agency.