1. Introduction

After the disposure of radioactive waste in the deep geological disposal system, various thermal, hydraulic, mechanical, and chemical complex processes, called multi-physics processes, occur within the system. The multi-physics processes affect the physical properties of engineered barrier system (EBS) and natural barrier system (NBS) constituting the disposal system, and directly or indirectly affect the migration of radionuclides released from radioactive wastes [1]. Therefore, the essential of total system performance assessment for the disposal system is to predict the migration of radionuclides released from radioactive wastes considering multi-physics processes not only related to each other but also occurring simultaneously and hence a set of partial differential equations needs to be solved together.

However, it is hard to obtain analytic solutions for partial differential equations due to the complicated structure of the disposal system, complex phenomena and diversity of materials and physical properties and thus numerical solutions are required to be obtained using numerical analysis methods such as, finite difference method, finite volume method and finite element method (FEM). Among the numerical analysis methods, FEM has been widely used in several codes for the multi-physics simulation of subsurface environment such as OpenGeoSys [2], and iCP [3]. Korea Atomic Energy Research Institute (KAERI) also has developed APro [4] to perform the FEM based numerical analysis for total system performance assessment (TSPA).

The FEM divides the whole domain into a finite number of elements and defines local basis functions for each element. Then, the numerical solution is obtained with a linear combination of those local basis functions. In other words, it is a method of obtaining the numerical solution by using as many local basis functions as a finite number of elements and obtaining the coefficients for the linear combination. Therefore, as the number of elements increases, the size of linear system and computing time increases but the accuracy of numerical solution increases.

In this regard, in the FEM based analysis, the mesh element quality such as the mesh size, number of meshes, and mesh shape affects not only the accuracy of numerical solutions but also the memory requirement and computing time. Low mesh resolution induces inaccurate results or convergence issues, while excessively high mesh resolution leads to huge computational burden. In other words, using a highquality mesh may be helpful for reducing numerical errors and producing reference solutions but may cause waste of computing resources and time.

Therefore, in this study, in order to efficiently and accurately perform TSPA using APro, a mesh stability study is performed prior to full-scale analysis. The mesh stability study determines the appropriate mesh structures varying size of elements, which result in efficient simulations in aspects of calculation time and memory requirement with reasonable numerical error. In addition to that, based on the results of mesh stability study, an adaptive mesh refinement is proposed employing the finest mesh structures for the initial time-period only.

2. Modeling Specification for the Mesh Stability Study

DECOVALEX-2023 is an international collaboration research project focusing on the development of understanding, models and codes in complex multi-physics phenomena in geological and engineering applications. Specifically, the Task F [5] of the project involves the comparison of models and methods used in post-closure performance assessment of deep geological repositories. In this regard, a generic reference case for a mined repository in fractured crystalline rock is proposed for the simulation of stationary groundwater flow and time-dependent conservative transport of two tracers.

Although the geologic setting for the reference case covers not only the EBS but also the NBS, in this study, only the region of EBS and near-field rock region are considered. The repository comprises 50 deposition tunnels branching off two parallel access tunnels and 50 deposition holes are emplaced in each disposal tunnels, which results in a total of 2,500 deposition holes as shown in Fig. 1 [5].

Fig. 1

Cross section of deposition tunnel (left) and repository layout (right).

In this study, three assumptions are introduced for the mesh stability study. First of all, it is assumed that the tracer release mechanism occurs at the surface of deposition holes and the transport simulations of the tracer are modeled to be started from the surface of deposition holes. Based on this, the mesh stability study focuses on disposal tunnels and near-field rock regions, which occupy most of the repository area. Second, the time-dependent conservative transport of tracers is simulated up to 100 years although the benchmark proposes 100,000 years since a period of 100 years is enough to compare the accuracy of solutions and calculation efficiency for various mesh structures. At last, the geometric structure of two parallel access tunnels is removed and replaced by a near-field rock region.

For the physics setup, Darcy’s law is employed for the simulation of stationary groundwater flow as follows, which the hydraulic head boundary conditions are applied to the surfaces of near-field rock region and no flow conditions are applied to the surfaces of deposition hole [5]:

where ρ and μ are density and viscosity of water dependent on the temperature, respectively, and Q_m is the mass source, u is the Darcy velocity vector, κ is the permeability, p is the hydraulic pressure, and g is the acceleration of gravity.

Once u is obtained from the stationary groundwater flow equation, the time-dependent transport simulations are performed by solving the advection-diffusion equation as follows:

where D_i, R_i and c_i are the effective diffusion coefficient, reaction rate and concentration of the species i, respectively, and ε_p is the porosity. Since the transport simulations of tracers are modeled to be started from the surfaces of deposition hole, the flux boundary conditions are applied to the surfaces of deposition hole preserving the total inventory of tracers. At the surfaces of near-field rock region, outflow boundary conditions that the tracers would be transported by the fluid motion are applied [6]. The permeability and porosity of tunnel backfill were set to single values of 1.5×10⁻²⁰ m² and 0.46 [5], respectively. And, those of rock domain were obtained as the equivalent continuous porous medium (ECPM) model which was upscaled from the discrete fracture network (DFN) model of the stochastically generated fractures. The effective diffusion coefficients of the tunnel backfill and the rock were set to single values of 1.6×10⁻¹⁰ m²·s⁻¹ and 2.0×10⁻¹⁴ m²·s⁻¹ [5], respectively. As shown in Table 1, two tracers with different release mechanism and inventory per waste package are simulated. Initially, the domain is empty of tracers everywhere except in the waste packages. The instant release mechanism is modeled with Gaussian pulse function and the congruent release mechanism is modeled with a fractional rate.

Table 1

Tracer inventories and release mechanisms

For the mesh generation, the maximum element sizes for the disposal tunnels and surrounding bedrock are adjusted in APro. The geometric specifications of model domain and the mesh generation specifications are shown in Tables 2 and 3, respectively. In case of the disposal tunnels, the maximum element size is set as 5 m considering the structure of 4.8 m height and 4.2 m width as shown in Fig. 1. Also, in case of the near-field rock regions, the maximum element size is set as 50 m since the center-to-center distance between disposal tunnels are 40 m as shown in Fig. 1. Once the mesh structure is generated, the number of meshes per unit volume is used as a measure for the mesh generation criterion. It is empirically demonstrated that changes in the computing time and memory requirements highly follow changes in the number of meshes per unit volume rather than changes in the maximum element size or the number of meshes. As shown in Table 2, although the rock domain occupies most of the repository area, the number of meshes per unit volume of the rock domain is smaller because it has a relatively large mesh size compared to the disposal tunnel. In Table 3, ‘CASEID’ denotes various model cases which have different mesh structures, some of which are shown in Fig. 2. Note that ‘T’ denotes the mesh stability case for disposal tunnel and ‘R’ denotes the mesh stability case for near-field rock region.

Table 2

Geometric specifications of model domain

Table 3

Mesh generation specifications of the stability study for disposal tunnel and near-field rock region

Fig. 2

Mesh structures for the cases of ‘T0’ (left) and ‘R0’ (right).

Total height of the repository is 100 m from 500 m to 600 m on the z-axis, which the floor of the disposal tunnels is at z=550 m. In x- and y-directions, the dimensions of the repository are 1,000 m from x=520 m to x=1,520 m and 696 m from y=652 m to y=1,348 m respectively. In order to estimate the accuracy of analysis results, the points at equal intervals are selected in the three-dimensional repository area. The starting and ending points, the distance between each point and the number of points are given in Table 4. Also, four representative points considering the groundwater flow results are selected as shown in Fig. 3.

Table 4

Selection of estimation points at equal intervals

Fig. 3

Four representative estimation points considering the groundwater flow results.

For the accuracy, the relative discrepancy for absolute difference of concentration for each time step is measured where the reference value is set from the results of finest mesh structure as shown in Equation 3:

where c_m,i,t is the concentration of tracer i at time step t for the mesh stability study index m and c_m,0,t stands for the case employing the finest mesh structure.

For the calculation efficiency, computing time and memory requirement are compared.

3. Results and Discussion

3.1 Results of Mesh Stability Studies for Disposal Tunnel and Near-Field Rock Regions

As shown in from Fig. 4 to Fig. 7, the mesh stability study results at four representative points show larger difference of tracer concentration depending on the tracer release mechanism rather than the mesh structure of disposal tunnel or near-field rock region. In case of the instant release mechanism, the tracer concentrations varying mesh structures show a large difference compared to the results using from the finest mesh structures, denoted as ‘T0’ and ‘R0’. Contrary to that, in case of the congruent release mechanism, there is no significant difference in the tracer concentrations obtained from varying mesh structures. This tendency appears in all cases for the representative points, and it is clearly shown that a large numerical error in the tracer concentration is mainly induced from instant release mechanism.

Fig. 4

The concentration of tracer simulated with instant release mechanism at four representative points varying the mesh structure of disposal tunnels.

Fig. 5

The concentration of tracer simulated with congruent release mechanism at four representative points varying the mesh structure of disposal tunnels.

Fig. 6

The concentration of tracer simulated with congruent release mechanism at four representative points varying the mesh structure of near-field rock region.

Fig. 7

The concentration of tracer simulated with congruent release mechanism at four representative points varying the mesh structure of near-field rock region.

The average of relative difference of tracer concentrations for all the estimation points, i.e. 504 points in whole repository area, is shown in Figs. 8 and 9. In the figures, a large difference in the results from different tracer release mechanisms is clearly shown. The numerical solutions are unstable in the early stage of the simulation time-period because the initial concentration of tracers is zero and the amount of release is fairly small compared to the size of the repository. As a result, the numerical errors for both instant and congruent release mechanisms fluctuate to large values. Also, for the instant release mechanism, despite the Gaussian pulse function is used for relaxing the step change, there is a large numerical error of tracer concentrations in the early stage of the simulation time-period more than 4% and it is maintained for the remaining time-period. In case of the congruent release mechanism, the initial numerical error is not maintained over time and is reduced to less than 2%.

Fig. 8

Average of relative difference for all the estimation points in case of instant (left) and congruent (right) release mechanism according to the mesh stability study for disposal tunnel.

Fig. 9

Average of relative difference for all the estimation points in case of instant (left) and congruent (right) release mechanism according to the mesh stability study for near-field rock region.

Tables 5 and 6 present the comparison of computing time and memory requirement according to the mesh stability study for disposal tunnel and near-field rock region respectively. Considering the case of congruent release mechanism only, it might be recommended to use the coarsest mesh structure because the coarser mesh structures are applied, the less computing time and memory requirement are, and the accuracy of numerical solution is maintained. However, in case of the instant release mechanism, the use of coarser mesh structure reduces the computing time and memory requirement, but the initial numerical errors tend to remain unreduced during the simulation time-period. Therefore, when considering two mechanisms together, it might be prudent to use coarser mesh structure.

Table 5

Comparison of computing time and memory requirement according to the mesh stability study for disposal tunnel

Table 6

Comparison of computing time and memory requirement according to the mesh stability study for near-field rock region

3.2 Enhancement of Accuracy Employing Adaptive Mesh Refinement

The discrepancy in the early stage of the simulation time-period is maintained for remaining time-period and its effect is non-negligible especially for the cases of instant release mechanism. Therefore, so-called the adaptive mesh refinement is proposed in order to reduce the discrepancy in the early time-period. The adaptive mesh refinement employs the finest mesh structure in the early time-period and coarser mesh structure for the remaining time-period. In this study, among the 100-year period, a time point after 20 years when the discrepancy starts to be maintained are set as a dividing point.

Note that solid and dotted lines in Fig. 10 represent the results without and with adaptive mesh refinement respectively. Since the adaptive mesh refinement cases identically employ the finest mesh structure in the early stage of the time-period, the relative differences for all cases up to 20 years are identical and plotted in gray dotted line. As shown in Fig. 10, in case of the instant release mechanism, the relative differences for the adaptive mesh refinement cases are significantly reduced from more than 4% to less than 2%. For the congruent release mechanism, the adaptive mesh refinement slightly improves the accuracy especially for the cases employing coarser mesh structures, denoted as T4-AM and T5-AM, but the amount is not significant.

Fig. 10

Significant reduction of relative difference with the adaptive mesh refinement.

3.3 Comparison of Computational Efficiency With Fixed Time Intervals

In this study, the backward differential formula (BDF) is employed as a temporal discretization method and optimized time intervals for the BDF method are taken by COMSOL [6]. For the simulation cases without adaptive mesh refinement, the optimized time intervals are taken for entire simulation time-period. On the other hand, when the adaptive mesh refinement is employed, the optimized time intervals are taken for the period up to 20 years and for the period from 20 years to 100 years, respectively. Therefore, the time intervals used in the simulation cases without and with adaptive mesh refinement are different, which can affect the estimation of computational efficiency. In this regard, the time intervals used in the case employing adaptive mesh refinement are manually set as same as the time intervals used in the cases employing the finest mesh structure. Therefore, only the difference in computational efficiency due to the difference in mesh structure can be compared.

Fig. 11 shows the relative differences of tracer concentrations for the cases using adaptive mesh refinement with different time intervals, and Table 7 summarizes computing time for the cases. Note that ‘AM’ and ‘FT’ stand for adaptive mesh refinement and fixed time intervals respectively. For the cases from ‘T1-AM-FT’ to ‘T5-AM-FT’, the computational efficiency according to the mesh structure can be clearly compared and the computing time is significantly reduced to under 50%.

Fig. 11

Comparison of relative difference of tracer concentration according to the introduction of fixed time intervals while the adaptive mesh refinement is employed.

Table 7

Comparison of computing time according to the adaptive mesh refinement and fixed time intervals

In aspect of accuracy, for both instant and congruent release mechanisms, the effect of applying the fixed time intervals on the relative difference results is not significant. Although there is a slight increase in relative difference results, the value is still less than 2% and 1% for instant and congruent release mechanisms, respectively.

4. Conclusions

In this work, the mesh stability study for disposal tunnel and near-field rock regions is performed using DECOVALEX- 2023 Task F reference case which models the stationary groundwater flow and the time-dependent conservative transport of two tracers in the deep geological repository. The stability analysis results are confined to the hydraulic and conservative transport, but are expected to be helpful in setting up an efficient mesh in the model development for the TSPA of a deep geological repository that addresses the THMC complex processes. In order to estimate the accuracy varying mesh structure, the tracer concentrations at 504 points with the same interval and 4 specific points considering the groundwater flow feature are selected. And the relative difference of tracer concentrations, derived from the varying mesh structure and the finest mesh structure is compared. For the calculation efficiency, computing time and memory requirement are compared. As a result of the mesh stability study, it is shown that the computing time and memory requirement are reduced as the coarser mesh structures are used, but the initial numerical error is nonnegligible and is maintained for remaining simulation timeperiod in case of the instant release mechanism. In order to resolve this problem, the adaptive mesh refinement which employs the finest mesh structure in the early time-period and coarser mesh structures for the remaining time-period is proposed. With the adaptive mesh refinement for the instant release mechanism, the relative difference of tracer concentrations significantly reduced from about 8% to less than 2% and the computing time is reduced to less than 50% based on same time intervals.