1. Introduction

Nuclear technology is used in a wide range of fields, including power generation, medicine, industry, and research. In medicine, radioisotopes are used to diagnose and treat disease. In industry, radiation technology plays a vital role in non-destructive testing, sterilization and process control. In research, it is used for material structure analysis, radiocarbon dating and other scientific studies. As a vital source of energy in modern society, nuclear power generation offers two significant advantages. First, it provides a stable supply of electricity regardless of climate or seasonal changes. Second, it produces minimal greenhouse gas emissions during operation, making it an effective energy source for combating climate change. However, the production and utilization of nuclear energy inevitably generate radioactive waste. The safe disposal of radioactive waste is essential for the long-term use of nuclear energy.

The classification of radioactive waste is based on the standards of the International Atomic Energy Agency (IAEA), but is adapted to the specific waste generation characteristics and regulatory requirements of each country. In South Korea, radioactive waste is managed under a fourtier system categorized by radioactivity levels and half-life: high-level, intermediate-level, low-level, and very low-level waste. High-level radioactive waste primarily consists of spent nuclear fuel, which requires long-term management against decay heat. Intermediate-level waste includes various devices and components generated during reactor operations. Low-level and very low-level radioactive waste encompass materials such as work clothing, decontamination residues, and wires produced during the operation of nuclear power plants and the use of radioactive isotopes. This four-typed classification system enables tailored treatment and disposal methods based on the risk and characteristics of the waste, contributing to the establishment of an efficient and safe management framework.

According to the International Atomic Energy Agency standards [1], radioactive waste should be isolated from the ecosystem through a multi-barrier system composed of engineered barriers, such as concrete walls, and the host rock of deep geological formations. The multi-barrier system consists of disposal containers, engineered barriers, and natural bedrock. This triple-layered barrier is designed to prevent the leakage of radioactive materials and to ensure safety without human intervention.

The Wolsong Low and Intermediate-Level Waste (LILW) Disposal Centre is in operation to manage domesticallygenerated intermediate- and low-level radioactive wastes. The phase 1 underground disposal facility, which began its full operation in 2015, is designed for the safe disposal of intermediate-level and lower-level radioactive waste. It consists of six silos, each 50 meters high, located at elevations between E.L. −80 m and −130 m. Within the same site, a near-surface disposal facility is under construction for the disposal of low-level and very low-level radioactive waste [2]. This facility is situated near the surface and is designed to prevent the release of radioactive materials through the use of natural or engineered barriers. This approach is widely adopted, with similar facilities operating in the United Kingdom, the United States, France, and Japan.

When radioactive waste comes into contact with water, numerous chemical and physical reactions can occur, potentially leading to the leakage of radioactive materials into the surrounding environment. Through this process, water acts both as a physical solvent and as a mediator for chemical reactions, potentially facilitating the transport of radioactive substances within the waste. Therefore, it is crucial to thoroughly analyze interactions with water and implement measures to control these interactions during the design of disposal facilities to ensure safety. For the safe disposal of radioactive waste, a multi-layer cover system (MLCS) should be executed right after the closure of nearsurface disposal facilities. The cover system consists of several layers, each with a distinct function. The vegetation layer on top prevents erosion, the drainage layer diverts surface water, the impermeable layer blocks infiltration, and the buffer layer protects the underlying disposal facility to ensure stability.

This study evaluated infiltration rates under various scenarios to assess the long-term performance of MLCS. Environmental factors such as rainfall intensity, topographic changes, soil characteristics, and climate change were considered for each scenario. Particular focus was given to analyzing the impacts of extreme weather conditions, natural disasters like earthquakes, and long-term topographic changes on the performance of the cover system. The study aimed to verify the suitability of the cover system design and identify input parameters for the safety assessment of disposal facilities to evaluate long-term radiological impacts. This systematic approach is expected to improve the reliability of near-surface disposal facilities and minimize the environmental impacts of radioactive waste.

2. Methodology

2.1 Conceptual Model

The near-surface disposal facility is constructed on an unsaturated layer of rock overlying a sedimentary bedrock. MLCS, as shown in Fig. 1, is installed to prevent water infiltration after closure of the disposal. This system consists of topsoil and side slopes. The topsoil layer supports vegetation and provides protection against erosion, animal and plant encroachment, and prevent inadvertent human access. Below the topsoil, several layers are designed to act as capillary and permeability barriers based on the permeability coefficients of the soils. They work together to ensure the effective containment of radioactive materials.

Fig. 1

Design of multi-layer cover system for near-surface disposal facility.

Precipitation is lost to the atmosphere through evapotranspiration, as illustrated in Fig. 2. A portion of the remaining precipitation manifests as surface runoff, which flows along the topsoil layer of the disposal cover and discharges at the surface. The infiltration water, which penetrates below the surface, represents another portion of the precipitation and can be expressed using the following equation:

Fig. 2

Three-dimensional conceptual model and water balance of disposal cover system.

Here, I represents the infiltration amount, P the precipitation amount, ET the evapotranspiration amount, and SR the surface runoff amount. This equation is utilized to quantitatively evaluate the amount of water transitioning from precipitation to infiltration, serving as a critical factor in assessing the hydrological performance of the disposal cover.

The 2D cross-sectional conceptual model of the disposal cover is shown in Fig. 3(a). The near-surface disposal facility is located directly above the prepared ground level, assuming that the infiltration beneath the surface exhibits a vertical and symmetrical pattern in the 2D section. In the model, the red dashed line represents the hydraulic divide. Under conditions where the hydraulic properties of each layer are homogeneous and vertical flow is dominant, horizontal infiltration behavior across the hydraulic divide is assumed not to occur. This assumption simplifies the analysis of infiltration behavior and evaluation of the hydraulic properties of the disposal cover, playing a crucial role in assessing the balance and stability of flow centered around the hydraulic divide.

Fig. 3

Two-dimensional conceptual model and water balance of disposal cover system.

The 2D longitudinal cross-sectional conceptual model of the disposal cover is illustrated in Fig. 3(b). When the disposal cover functions as intended, precipitation passes through the topsoil layer, the gravelly sand layer, and the pea gravel layer to reach the first drainage layer. Infiltrated water reaching the first drainage layer, composed of sand, primarily flows horizontally due to the low permeability of the underlying clay layer, with only a minimal amount infiltrating into the clay. Water entering the second sand layer comes from infiltration through the upper clay layer. Similar to the first drainage layer, most water in this layer also moves horizontally, with a small portion penetrating the lower clay layer. This multi-layer structure significantly limits the amount of infiltrated water reaching the disposal facility via the sand layer beneath the second clay layer. As a result, the quantity of water infiltrating into the disposal facility is expected to be extremely minimal.

2.2 Governing Equations

Buckingham extended Darcy’s law [3], originally formulated for saturated soils, to unsaturated soils by proposing the Darcy-Buckingham equation as follows:

Here, q represents the Darcy velocity (cm∙s⁻¹), K(θ) is the unsaturated hydraulic conductivity (cm∙s⁻¹), h is the pressure head (cm), and z is the elevation head (cm). This equation indicates that flow occurs according to the gradient of the pressure head and the elevation head.

Richards applied the law of mass conservation to the Darcy-Buckingham equation to derive the governing equation for unsaturated flow [4]. The law of mass conservation is expressed as:

Here, θ represents the soil’s volumetric water contents (cm^³·cm^−³), t denotes time (s), and ∇∙q refers to the divergence of the water flux (cm^³·s⁻¹·cm^−³). Equation (3) indicates that the rate of change in water content per unit volume over time is equal to the divergence of water flux passing through that volume.

Following equation represents the divergence form of the Darcy-Buckingham equation, which describes the flow of water in unsaturated soils:

This expression demonstrates how the divergence of flow velocity is influenced by the spatial variations in pressure head and the unsaturated hydraulic conductivity. To further expand, the divergence form can also include additional terms accounting for gradients of both pressure and gravitational effects. Thus, Equation (4) can be reformulated as follows:

In this form, the first term, −K(θ) ∇²h, represents the diffusion-like flow caused by the curvature of the pressure head distribution. The second term, −K(θ) ∙ (∇h + e_z), includes contributions from both the gradient of the pressure head and the gravitational head (z), where e_z is the unit vector in the vertical direction.

Substituting Equation (5) into the mass conservation law (Equation (3)) gives the Richards equation, which governs the flow in unsaturated soils:

This formulation shows that changes in water content (∂ θ ∂ t) are driven by the combined effects of hydraulic conductivity, pressure head gradients, and gravitational forces. When K(θ) is assumed to be independent of depth (z), the third term vanishes, simplifying the equation.

The Richards equation, as shown in Equation (6), integrates the Darcy-Buckingham equation with the principle of mass conservation to model unsaturated water flow. It is widely used in soil physics and hydrology to describe the transient behavior of water in the vadose zone, where saturation and conductivity vary dynamically.

Van Genuchten (1980) proposed a widely used relationship to describe the hydraulic properties of unsaturated soils as follows [6]:

Here, S_e is the effective saturation, θ_r is the residual water content (cm^³·cm^−³), θ_s is the saturated water content (cm^³·cm^−³), α is the sensitivity parameter for the pressure head (cm^⁻¹), n is a dimensionless constant that characterizes the hydraulic properties of the unsaturated medium, m is defined as 1−1/n.

The hydraulic conductivity function for unsaturated soils is given as follows:

Here, K(S_e) is the unsaturated hydraulic conductivity (cm∙s⁻¹), K_s is the saturated hydraulic conductivity (cm∙s⁻¹), l is an empirical constant, typically assigned a value of 0.5.

Combining the Richards equation with the Van Genuchten model results in the following equation:

Equation (9) serves as the governing equation for describing water movement in unsaturated soils over time and space. It is used to predict the distribution of soil water content and pressure head.

2.3 Numerical Modelling

2.3.1 Geometric Model

The near-surface disposal facility involves installing concrete disposal vaults on stable bedrock. Preliminary site investigations reveal that the unsaturated bedrock beneath the planned location has a thickness of approximately 70 meters. For modeling purposes, the unsaturated bedrock was set to extend from the ground level EL. 107 m to EL. 37 m, with the lowest boundary of the modeling domain located 10 meters below the groundwater table at EL. 27 m, as shown in Fig. 4(a).

Fig. 4

Numerical modeling for multi-layer disposal cover infiltration.

To ensure compatibility with MLCS design and numerical stability, quadrangle mesh elements ranging from a minimum of 0.134 m^² to a maximum of 10.09 m^² were generated, as depicted in Fig. 4(b). The program used for this analysis was FEFLOW, a widely recognized tool for simulating flow and transport processes in porous media.

2.3.2 Boundary Conditions

At the bottom boundary of the model in the saturated zone, a fixed head condition (Type 1 or Dirichlet boundary) was applied. At the lateral boundaries, a no-flow condition was implemented to represent impermeable barriers. At the top boundary, a hydraulic head constraint was applied to simulate rainfall infiltration. The applied boundary conditions are shown in Fig. 5(a).

Fig. 5

Numerical model for infiltration through multi-layer disposal cover.

2.3.3 Barrier hydraulic properties

The material properties for each layer in the engineered barrier system, including the disposal cover and disposal cells, were applied based on values referenced from NUREG/ CR-6114 [5], as shown in Table 1, while Fig. 5(b) illustrates the stratigraphic distribution of the analysis model, where the waste within the disposal cells was treated with properties equivalent to those of crushed stone. The heterogeneous hydraulic properties of the bedrock were derived using Frac- Man, a specialized tool for modeling fracture networks and their influence on hydrological characteristics.

2.3.4 Precipitation condition

To estimate the infiltration rate through the multi-layer cover of the near-surface disposal facility, rainfall conditions at the top boundary of the model were set using meteorological observation data from the Ulsan region. Climate normals (30-year averages) were calculated based on rainfall records spanning 1985 to 2014, and recharge rates were derived and applied as boundary conditions. As shown in Fig. 6, actual daily rainfall, monthly averages, and annual averages were input in the model to assess and compare the resulting infiltration rates.

Fig. 6

Rainfall conditions based on climatological normals at Ulsan weather station (1985−2014).

2.3.5 Modelling scenarios

Once the site is closed, a disposal cover will be installed. During the 300-year post-closure management period, systematic maintenance will be carried out to ensure the performance of the cover system. Considering that the hydraulic properties of the materials within MLCS may gradually degrade over time, the maintenance program is designed to minimize this degradation and maintain its critical functions. For numerical analysis, it has been reasonably assumed that the materials will maintain their initial hydraulic properties throughout the management period under this maintenance regime.

To enhance computational efficiency, only half of the symmetric structure was activated for the analysis. Sensitivity analysis of infiltration flow was conducted by varying input parameters (rainfall patterns and bedrock properties) under the following four cases:

2.3.6 Analysis points for rainfall infiltration

To analyze the rainfall infiltration rate through the disposal cover system, the flow characteristics of infiltrating water were considered, and boundary interfaces of the area of interest were defined in each model. A water budget was calculated for them, representing the balance of water inflows and outflows within the system. The infiltration rate was then determined by dividing the water budget by the effective area, yielding the amount of water infiltrating into the area per unit time.

To quantitatively evaluate the amount of infiltrated water through the disposal cover system during the postclosure management period, three boundary interfaces were defined, as shown in Fig. 7. Firstly, at the topsoil upper boundary interface, the amount of rainfall infiltration recharged through the disposal cover system was assessed. Secondly, at the lower boundary of the second clay layer, the performance of MLCS was evaluated to verify whether it meets design specifications. Thirdly, at the disposal cell inner wall boundary interface, the quantity of infiltrated water entering the radioactive waste storage area was assessed.

Fig. 7

Boundary conditions for infiltration water budget analysis.

3. Modelling Results

3.1 CASE #1

Fig. 8 illustrates the time-series variation of infiltration rates per unit area, calculated based on the water budget for each boundary interface defined in Fig. 7 under CASE #1. The hydraulic properties of the bedrock were assumed to be equivalent to those of sand, and the annual average rainfall infiltrating through the topsoil layer of the disposal cover was set to 339.2 mm∙year⁻¹.

Fig. 8

Zonal infiltration rates during closure period (CASE #1).

Rainfall infiltrates the topsoil and percolates through six layers to the upper clay. The infiltration rate initially increases but converges to an average of 5.9 mm∙year⁻¹ (18.3% of the design criterion of 32 mm∙year⁻¹) 40 years after the facility’s closure. The infiltration rate into the concrete vault stabilizes approximately 70 years after closure at an average rate of 0.03 mm∙year⁻¹, confirming the longterm effectiveness of the cover.

The saturation and pathway analysis results in Fig. 9 illustrate water movement dynamics within the disposal system. The left panel shows the pathlines of infiltration water (blue lines), while the right one displays the saturation distribution, ranging from fully dry (brown) to fully saturated (blue). Most rainfall infiltration passing through the topsoil flows along the sand layers, which are originally designated for drainage, leading to rapid surface runoff at the sloped edges of the multi-layer cover.

Fig. 9

Analysis of flow pathlines and saturation distribution at different post-closure periods (CASE #1).

The figure captures a moment in the continuous process of water migration under gravity, particularly within the vadose zone, where water initially fills the void over time before moving downward due to gravitational forces. The saturation levels within the clay layers gradually increased, while no significant changes were observed inside the concrete disposal vaults. After 300 years, a hydraulic head formed between the disposal vaults, but the internal changes in infiltration remained negligible. These findings confirm the effectiveness of the multi-layer disposal cover in redirecting infiltrated water and minimizing its impact on the interior of the disposal facility over the long term.

3.2 CASE #2

Fig. 10 depicts the time-series variation of infiltration rates per unit area, calculated based on the water budget at each boundary interface in Fig. 7 under CASE #2. In this scenario, the hydraulic properties of the bedrock were assumed to be equivalent to those of sand.

Fig. 10

Zonal infiltration rates during closure period (CASE #2).

Rainfall infiltration into the topsoil of the disposal cover was modelled using monthly average rainfall values, which accounted for rainfall frequency. The rainfall conditions ranged from 0.0 mm∙day⁻¹ on sunny days to a maximum of 10.6 mm∙day⁻¹ on the rainy days.

Rainfall infiltrates the topsoil and percolates through six layers to the overlying clay. After closure of the facility, the infiltration rate under the lower clay initially increases, but stabilizes at an average of 6.8 mm∙year⁻¹ (21.25% of the design criterion, 32 mm∙year⁻¹) after 30 years. The infiltration rate into the concrete vault stabilizes at an average of 0.03 mm∙year⁻¹ (maximum 0.04 mm∙year⁻¹) approximately 50 years after closure, confirming the long-term effectiveness of the cover.

The saturation and flow path analysis results shown in Fig. 11 indicate that drainage proceeded smoothly, similar to CASE #1. The increase in saturation levels within the clay layers was minimal. After 300 years, a hydraulic head was observed to have formed between the disposal vaults. However, changes in internal infiltration remained negligible. These results represent a characteristic snapshot of unsaturated flow patterns, where water movement is primarily governed by gravitational forces and capillary action in the vadose zone.

Fig. 11

Analysis of flow pathlines and saturation distribution at different post-closure periods (CASE #2).

3.3 CASE #3

Fig. 12 shows the time-series variation of infiltration rates per unit area, calculated based on the water budget at each boundary interface in Fig. 7 under CASE #3.

Fig. 12

Zonal infiltration rates during closure period (CASE #3).

In this scenario, the hydraulic properties of the bedrock were assumed to be equivalent to those of sand. Rainfall infiltration into the topsoil was modeled using daily rainfall data recorded at the Ulsan meteorological station from 1985 to 2014. These observations were repeated in 30-year cycles to reflect long-term variability in daily rainfall patterns.

Rainfall infiltrates the topsoil and percolates through six layers to the upper clay. During initial operation, the infiltration rate increases below the lower clay, but stabilizes at an average of 6.6 mm∙year⁻¹ (20.63% of the design criterion of 32 mm∙year⁻¹) within 30 years. The infiltration rate into the concrete vault stabilizes at an average of 0.03 mm∙year⁻¹ (maximum 0.04 mm∙year⁻¹) approximately 60 years after closure, confirming the long-term effectiveness of the cover.

The saturation and flow path analysis results in Fig. 13 show that, similar to CASE #1 and CASE #2, infiltrated water was primarily discharged along the sand layers. Over the long term, the saturation level of the clay layers in the multi-layer cover gradually increased. However, the interior of the concrete disposal vaults maintained its initial state with no significant changes observed. After 300 years, a hydraulic head developed between the disposal vaults. Nevertheless, changes in internal infiltration remained negligible. The figure provides a snapshot of ongoing water dynamics and saturation changes within the disposal system.

Fig. 13

Analysis of flow pathlines and saturation distribution at different post-closure periods (CASE #3).

3.4 CASE #4

Fig. 14 illustrates the time-series variation of infiltration rates per unit area, calculated based on the water budget at each boundary interface in Fig. 7 under CASE #4. Unlike previous conditions, CASE #4 incorporates the Discrete Fracture Network (DFN) concept to represent bedrock heterogeneity. The DFN model was constructed based on field-investigated fracture characteristics such as orientation, frequency, and size distribution. The equivalent hydraulic properties from this model were applied to the numerical simulation. Rainfall conditions were represented using monthly average data with precipitation frequency.

Fig. 14

Zonal infiltration rates during closure period (CASE #4).

Rainfall was observed to pass through the six layers before reaching the upper clay. After closure, the infiltration rate under the lower clay initially increases but stabilizes at an average of 5.1 mm∙year⁻¹ (15.94% of the design criterion, 32 mm∙year⁻¹) within 30 years. Internal infiltration stabilizes at an average of 0.03 mm∙year⁻¹ (maximum 0.04 mm∙year⁻¹) approximately 50 years, confirming the effectiveness of the cover.

The saturation and flow path analysis results in Fig. 15 indicate that MLCS satisfied its design performance criteria. Even with fracture characteristics incorporated into the model, the overall behavior of infiltrated water was similar to the patterns observed in homogeneous media (CASE #1 to #3). This suggests that the performance of the disposal cover is not significantly sensitive to changes in the hydraulic properties of the underlying bedrock. Fig. 15(a) shows variations in saturation levels along fault structures, indicating potential differences in permeability compared to the surrounding media. With the Discrete Fracture Network (DFN) model assigning varying hydraulic conductivity to all elements, selective flow pathways could develop along regions of high permeability. While these fracture-influenced flow patterns do not significantly impact the cover system performance, they may affect water infiltration and radionuclide transport from the vadose zone to the saturated zone, requiring further investigation.

Fig. 15

Analysis of flow pathlines and saturation distribution at different post-closure periods (CASE #4).

4. Findings

This study quantitatively evaluated the infiltration rates through MLCS of a near-surface disposal facility over a post-closure management period of 300 years. The evaluation considered the hydraulic properties of the layers and various rainfall conditions, under the assumption that the initial hydraulic properties and geometry of MLCS remain unchanged.

The analysis demonstrated that MLCS effectively fulfilled its designed drainage and impermeability functions under all cases, showing a clear reduction in infiltration between the rainfall entering the top soil and the infiltration discharged at the bottom of the lower clay. The infiltration delay increased progressively from the top to bottom layers, and infiltration rates exhibited an initial increase followed by stabilization over time. By incorporating DFN characteristics (such as fractures and faults) into the bedrock, the study identified potential preferential flow paths that could develop due to bedrock heterogeneity in realistic scenarios. Nevertheless, the distribution and variation of overall infiltration rates were similar across scenarios, indicating that the heterogeneous hydraulic properties of the bedrock had a limited impact on infiltration rate evaluation. Furthermore, the sensitivity differences among rainfall conditions (annual average, monthly average, and daily average) were minimal, suggesting that infiltration evaluations could be reliably conducted under CASE #2 (sand-like hydraulic properties with a monthly average rainfall pattern).

It is assumed that the initial geometry and material properties of the cover system remain unchanged. However, by focusing on different rainfall patterns and bedrock characteristics, this study may differ from actual long-term performance as it excludes factors like erosion, intense rainfall events, chemical degradation (e.g., changes in porosity and permeability due to cement deterioration), and structural stability. Therefore, additional design considerations, such as drainage facilities, should be investigated in follow-up studies that take these factors into account.

Additionally, it is essential to verify the long-term performance of the system through full-scale experimental facilities, trial execution, and the establishment of a robust long-term monitoring framework. To address uncertainties, studies involving probabilistic distributions of input parameters using methods like Monte Carlo simulations should be conducted to quantify uncertainties effectively.

These additional studies aim to provide stronger scientific evidence for the long-term safety of disposal facilities and serve as foundational data for demonstrating technical feasibility during the licensing and regulatory approval process.