Thermoelastic analysis of a geological repository with distributed decay heat sources by the image method in combination with a numerical integration scheme

: The disposal of nuclear waste involves thermo-mechanical reaction of the host rock to the buried waste – a distributed heat source that decays. To solve the problem within the half infinite space confined by the ground surface, an image method is developed. Specifically, a negative image of the heat source with the ground surface as the mirror, i.e. a mirrored heat sink is utilized so that the normal traction generated by the heat source can be counterbalanced, and a numerical scheme of integration of the classical Cerruti solution is developed to include the effect of tangential shear traction on the ground surface caused by the heat sources and their mirrored sinks. For a conceptual repository model, large thermal shear stress, tensile stress, and deformation occur at the corner, between adjacent drifts, and at the boundary of the repository area, respectively. For a prescribed thermal loading, it is more efficient to mitigate the thermo-mechanical effects through enlarging the pit spacing than increasing the drift spacing.


Introduction
Due to the radioactivity and toxicity of nuclear waste, a final disposal solution is required for a sustainable nuclear power program [1]. Constructing geological repositories is currently recommended as a feasible and reliable option to safely isolate radionuclides from the biosphere for very long time. For example, the Swedish Nuclear Fuel and Waste Management Company has submitted permit applications for the construction of a final repository in crystalline bedrock [2].
In the design of a geologic repository, the thermomechanical interaction induced by heat emitted from nuclear waste is a major factor. The heat can change the temperature field and create thermal stress in surrounding rocks [3]. The thermal stress may cause rock fractures increasing the possibility of groundwater intrusion into storage area [4]. The salt in groundwater may corrode the cylindrical container made up of copper or cast iron. Meanwhile, the high temperature environment can promote this corrosion process [5]. Once containers are damaged, radionuclides may quickly reach the ecosystem by groundwater flowing through fractures. To evaluate the safety of a potential repository, the coupled thermo-mechanical process is of interest.
Crystalline, as one kind of candidate rocks, has very low permeability. Thus heat propagates mainly by conduction in a geologic repository. Several appropriate local-scale models are developed to investigate the heat transfer or thermo-mechanical interaction around the central canister [6][7][8][9][10][11][12][13][14][15][16]. Although multiple barriers consisting of the container, buff, backfill and rock can be considered, these local-scale models are only available for the symmetrical position at any time scale or asymmetrical region at short-term scale. Designing a repository also needs some large-scale models which are capable of describing the variations of temperature and thermal stress at any location in the repository during the long-term disposal period. Mapping the canister on a rectangular plane, point or line, authors proposed several analytical or semi-analytical models for the repositoryscale heat transfer analysis [17][18][19][20]. Analytical [21] and regional numerical methods [22] are also adopted to analyze the repository-scale thermo-mechanical interaction based on the plane model. But, large errors of thermal stress-displacement exist near the disposal area on account to the significant geometric difference between the cylindrical canister and the rectangular plane. Therefore, a more efficient repository-scale model is of requirement for the thermo-mechanical assessment.
In the present paper, the thermal loading of a nuclear waste repsoitory is simplified into distrbuted and decaying line heat sources. To satisfy the constant temperature and stress-free bounday conditions at the ground surface, an image method in cobination with a numerical integration scheme is developed. The method is first comapred with existent lierature, and then applied to some hypothetical situations to examine the thermomechnical effects. Fig.1 illustrates a schematic of the KBS-3V repository [23]. In this study, the crystalline is assumed as a semiinfinite, isothermal, homogeneous, elastic medium.

Mathematical equations
Temperature-independent rock properties, uniform initial temperature field, constant ground surface temperature and rectangular disposal region are postulated as well.

Energy conversation equation
The heat propagation is dominated by conduction due to the low permeability of crystalline. The energy conversation equation is where θ is the temperature increment; ρ rock density; λ rock thermal conductivity; c rock specific heat; i= x, y, z; the subscript comma indicates partial differential; the superscript point represents the derivative of time t.

Initial and boundary conditions
The following initial and boundary conditions are used.
where the superscript x= (x, y, z).

Solution for an infinite space
Due to an instantaneous unit point heat source in the infinite region, Carslaw et al. [25] gave the temperature solution which met the Eqs.(1) and (5) except for the temperature boundary condition. Using the temperature solution, Dargush et al. [26] derived the thermoelastic solution which satisfied the Eqs.(2), (3), (4), (6) and (7) except for the stress boundary condition. Integrating the thermoelastic solution with respect to time, we obtain the thermal displacement u ch and stress σ ch induced by a constant unit point heat source at point x' and time t 0 .
With regard to a time varying point heat source, we divided the period t 0 -t into intervals of N. The heat source strength is assumed to be constant and equal to Q(t n ) during the n th time lag. Based on the solutions in Eq.(8), the thermal displacement u̅ and stress σ̅ are Except for σ yz and σ zx at the ground surface, other boundary conditions are satisfied by adding the effect of a heat sink with equal magnitude at point (xʹ, yʹ, -zʹ). However, these two shear stresses caused by the heat source are the same as ones resulted by the heat sink, respectively. Therefore, the solution of the following stress boundary problem should be superimposed.

Solution for a semi-infinite space
For the tangential force on a half-space sufrace, the theoretical solutions were found by Cerruti. The generalized integrals of the Cerruti solutions related to the load in x-direction are where the superscript denotes the load direction; ξ and η denote the boundary coordinates, and are parallel with xand y-axis, respectively; ( ) ( ) The improper integrals of the Cerruti solutions, related to the y-direction force can be calculated by using the abbove approach. The ground surface is discretized by rectangles of S. Assuming uniform load over each element, combining Eqs. (10) and (11) yields where s denotes the s th element; f c sx and f c sy the values of f x and f y at the element center, respectively; ξ s1 and ξ s2 the region limits of the element in x-direction; η s1 and η s2 the region limits of the rectangle in y-direction.
The abbove results are applied to derive the solutions for multiple line heat sources. Among heat sources of M, the m th one with a strength of Q m (t) is L m high and evenly divided into K m segments. The k th segment approximates to a point with the Simpson's weight of w k . The Simpson's weights are [28] [ ] The thermoelastic induced by all heat sources are The heat released by per canister is [18] where the unit of t is year. Some parameters are listed in Table 1[18].

Model comparison
The present model is compared with the one [21] which applies the planar heat source model and neglects the stress-free condition at the ground surface. The line L̅ ( z 1 =472m, z 2 =497m) shown in Fig.2 is selected for examination of the comparison.  Fig.3 illustrates the thermal effects along line L̅ after 10 years of closure. The compressive stress is defined positive while tensile stress negative. Good agreements can be found between the two solutions when the distance from the examination point to the central canister center exceeds 25m. But as the distance decreases, large differences of the two solutions can be observed due to the different heat source model used. Ref. [21] (σ xx ) Ref. [21] (σ yy ) Ref. [21] (σ zz ) Present (σ xx ) Present (σ yy ) Present (σ zz ) (a) 3 8

Thermal effects at the examination points
Six points in Fig.2 are chosen as examination points. Their coordinates and thermal effects are listed in Table  2 and shown in Fig.4, respectively.  Fig.4a shows the variation of horizontal stress σ xx with time. Although the stresses at points B 2 and C 2 are tensile before the 80 th year, their values are below 0.1 MPa. By contrast, the stresses at other four points are compressive. Among these points, the peak of compressive stress at point A 1 is the largest and exists during the 40 th ~ 70 th year with a value of 15.2MPa. Fig.4b illustrates the evolution of horizontal stress σ yy with time. Before the 100 th year, the stresses at points A 2 , B 2 and C 2 are tensile and less than 0.5MPa. In contrast, the stresses at other three points are compressive. The maximum of compressive stress peaks at these points is 21.5MPa occurring in the 20 th year at point A 1 . compressive after the 50 th year. Before the 70 th year, the stress at point C 1 is tensile and less than 0.5MPa. By comparison, the stresses at other points are compressive. Among these points, the compressive stress at point A 1 is the largest having two peaks. The first peak is 9.4MPa in the 3 th year and the other is 5.3MPa in the 700 th year.
The shear stresses at points A 1 , A 2 , B 1 and B 2 are nearly zero because these points close the symmetrical location in the repository. Fig.4d only plots the development of shear stresses with time at points C 1 and C 2 . The shear stress at point C 1 increases to the maximum of 2.1MPa in the 50 th year. The tangential stress at point C 2 is lower than 1.0MPa during the sealing period.  Fig.4. Thermal stresses and displacements at interest points.

Thermal effects on the examination plane
The present model is employed to investigate the thermal effects on the quarter of the disposal region plane(Plane  Fig.2. The thermal stress and displacement after 10 years of sealing are illustrated in Fig.5. The spatial differences of the thermo-mechanical effects are significant. In the vicinity of the central pit, the compressive thermal stress of 20 MPa occurs. Howerver, the tensile thermal stress of 1 MPa exists between adjacent drifts. Large shear stress of 2 MPa mainly occurs at the repository corner. There are large deformation near the repository boundary.

Influences of pit/drift spacing on the thermal effects
For constant heat power per unit disposal area, the values of d×D in three schemes are set as 4 m × 60 m (S 1 ), 5 m × 48 m (S 2 ) and 6 m × 40 m (S 3 ), respectively. Fig.8 illustrates the thermal stress of point A 1 in these scenarios. The thermal stress at point A 1 increases with the enlargement of pit spacing and decrease of drift spacing. Time t (year) S1 (σ xx ) S2 (σ xx ) S3 (σ xx ) S1 (σ yy ) S2 (σ yy ) S3 (σ yy ) S1 (σ zz ) S2 (σ zz ) S3 (σ zz ) Fig. 8. Thermal stress at point A1 of the three senarios.

Summary
To analyze the thermo-mechanical response of a repersitory, the thermal loadings are characterized as distributed decaying line heat sources, an image method in combination with a numerical integration scheme was developed. The results from several hypothetical calculation examples demonstrate that (1) The presented method is efficient to calculate the effects of the thermal loading pattern and the burial depth of a geological repository on the thermal stressdisplacement in the host rock.
(2) Due to the decay of the heat from the nuclear waste, the thermal stresses and displacements in the host rock increase at first and then decrease with time.
(3) The spatial differences of the thermo-mechanical effects are significant. Large thermal shear stress, tensile stress, and deformation occur at the repository corner, between adjacent drifts, and at the repository boundary, respectively.
(4) For a prescribed thermal loading, it is more efficient to mitigate the thermo-mechanical effects through enlarging the pit spacing than increasing the drift spacing. E3S 143 201430 ARFEE 2019 1 16 (2020) 6