1. Introduction

Among studies regarding the geological disposal of high-level radioactive waste, it is necessary to evaluate the performance of natural barriers that must maintain stability during long-term geological evolution. As a consequence of the difficulty in predicting the future geological environment, various research and construction methods for reliable evolution have been conducted globally. The modeling research must be conducted to specify the inferred geological evolutionary event and apply it to a model that simulates the disposal site. This modeling technique is crucial because it can predict the critical behavior of natural barriers caused by evolutionary events in the geological environment. To develop modeling technology for the simulation of each scenario, related research is being conducted targeting specific research sites. Research from Finland, Sweden, and Canada suggest a numerical method with a discrete element or finite element method and report evaluation results of the compossible behavior of the disposal site subsequent to the specific evolution process [1,2,3]. The numerical model was evaluated as a hydromechanically coupled model. Most models are constructed with the possibility of the re-activation of fracture zones and changes in the stress field in response to geostructural events such as glacial loads and earthquakes. Therefore, it is very important to construct an optimized coupled model to accurately evaluate the interactions and impacts of geological environmental changes. Various geological evolutionary processes and scenarios can be listed over time and designated as FEPs (features, events, and processes). It is essential to consider all features because each factor interacts with the others in terms of mechanics, hydraulics, tectonics, and bedrock heat content. However, it may also be problematic to integrate all the features into a quantitative model that predicts the evolutionary process. Fig. 1 is a systematic table showing the geological, rock-mechanical, hydraulic, geochemical, and thermal mechanism interactions in the disposal environment [4]. The interactions offer various physical phenomena and are complex. When the performance of the rock as a natural barrier to protect engineering is evaluated, the priority consideration is the hydraulic condition, with the mechanical state representing long-term geological environment changes. As shown in Fig. 1, the hydro-mechanical coupled process is a preferential concern for safety investigations, and the physical changes interact independently of geophysical features. Therefore, the evaluation model varies according to the rock type and geological scale. This study reviewed and systematized the principal features of long-term geological stability. In addition, as a case study to develop a reasonable model for the geological characteristics of Korea, various modeling methods with specific cases are summarized. Modeling techniques corresponding to the tectonic history and conditions of the Korean Peninsula were categorized and identified.

2. Hydro-Mechanical Evolution of the Natural Barriers

2.1 Categorization of FEPs and Configuration of an Interaction System

Since 2000, the Nuclear Energy Agency (NEA) has categorized FEPs that can be utilized for deep geological disposal, and has proposed a structure and systematic format. The FEPs list known as the NEA International FEP (IFEP) was recently revised to version 3.0, in 2019 [5]. Fig. 2 shows a schematic diagram of the relationship between the external FEPs proposed by the NEA and the internal FEPs corresponding to the disposal system components. Among 268 FEPs offered by NEA, the external FEPs were selected only for a related list of natural barrier evolution, and internal systemic FEPs, which may occur with the interaction between the physical components of the rock and external FEPs. As such, the FEP list arrangement and the formulation of each discipline are significant and essential initial stages for investigating the long-term safety of disposal systems [6].

Finland’s Posiva and Sweden’s SKB also classify and systematize FEPs, and propose a process that can occur because of the interaction of each feature in the natural barrier system. Based on the proposed interrelation matrix from Posiva, the processes and features related to the hydro-mechanical evolution of natural barriers included in the scope of this study were selected and are presented in Table 1. In particular, the relationship with the FEPs presented by NEA in Fig. 2 is indicated. Consequently, it was determined that all geosphere processes suggested by the NEA correspond to the groundwater flow characteristics and the disposal system stress state. The relationship between the properties and physical processes, newly marked in red in the table, is also characterized in this study. It was determined that the pore pressure effect could be related to the change in the rock mass stress state or the re-activation of discontinuous fractures based on the effective stress principle [7,8].

This study reviewed the mechanisms and interactions between rock characteristics and hydro-mechanical conditions. In the section below, the specific mechanism responsible for the evolution of hydro-mechanical behavior caused by a series of long-term evolutionary processes of natural barriers is reviewed in detail.

2.2 Mechanism of a Natural Barrier Evolution Process

The rock stress state is an important mechanism and a significant factor for determining the evolution and stability evaluation of natural barrier rocks. Geological structural characteristics and various pressure conditions may modify the stress state. Geological structures representing largescale discontinuities may be natural factors that induce stress-field disturbances. The representative evolution in pressure conditions is the result of tectonic changes (uplift or subsidence) that occur on long-term and large scales. Earthquakes occur on a local scale [9]. The main mechanisms driving the four hydro-mechanical evolutionary processes of natural barriers shown in Table 1 are briefly described in the following section.

2.2.1 In-situ Stress & Stress Redistribution

Although there is no internationally accepted common terminology to describe the stress states in rocks, Zang and Stephansson [10] presented different types of stress that can exist in the Earth’s crust (Fig. 3). The first hierarchial level of stress (A‒D) can be classified into two pairs: in-situ/ perturbed in-situ stresses and structural/perturbed structural stresses. In-situ stress means the stress state in the rock mass before any man-made disturbance, and the term structural stress is used instead if the rock mass is anisotropic or heterogeneous and therefore the stress distribution is affected by such natural geological structures. The term ‘perturbed’ is added when the in-situ or structural stress is changed by artificial disturbance such as excavation, explosion, pumping, etc.

On the second hierarchial level, each stress on the first level is subdivided into four types by its origin: gravitational, tectonic, residual/remnant, and terrestrial stresses (A1‒A4). For general applications in rock engineering, gravitational and tectonic stresses are dominant factors of the stress states in the rock mass. Gravitational stress increases with depth in the Earth’s crust, and is largely governed by the density distribution in the rock mass and affected by surface topography. Tectonic stress is largely driven by tectonic movements in the Earth’s crust, and can be subdivided according to the contributing factors at more detailed level, such as tectonic plates, plate bending & isostacy, and individual faults/fracture zones. These various rock stresses occur as the result of highly complicated processes and mechanisms. Because such geological and geomechanical processes and corresponding geostructures cannot be precisely determined across the wide temporal and spatial domains of interest, it is impossible to accurately compute the stress state from its driving forces. Therefore, in general, rock stress is determined as the best ‘estimated’ stress based on available information, such as field test data and results through various observation methods [11].

The stress state can be expressed as the magnitudes and directions of the maximum, intermediate, and minimum principal stresses. If stress disturbing factors such as the effects of large geological structures or surface topology are expected to be negligible in a rock volume of interest, vertical stress can be typically assumed to be one of the principal stresses; thus, in such cases, the maximum and minimum horizontal stresses are also considered as the principal stresses. As shown in Table 2, the stress regimes are classified as normal faulting (NF), strike-slip faulting (SSF), and reverse faulting (RF) regimes according to the magnitude orders of the three principal stresses [12].

Using the Mohr-Coulomb (M-C) failure criterion, which represents the critical level of stress state for shear failure of rock with maximum and minimum principal stresses, the magnitude range of the maximum and minimum principal stresses in which failure does not occur at a specific depth (i.e., a specific vertical stress) can be expressed mathematically according to the three stress regimes. The Mohr-Coulomb criterion equation is as follows:

where τ is the shear stress applied on the discontinuity plane which is judged for the criticality of shear slip, σ_n is the normal stress applied on the same plane, P_p is the pore pressure, σ ′ n is the effective normal stress, and c, μ, ϕ are the cohesion, friction coefficient, and friction angle, respectively, of the discontinuity. The Mohr-Coulomb criterion then can be applied and transformed to address the relationship between the maximum and minimum principal stresses which can be maintained without shear failure.

The stress regimes bounded by the above relation are shown on the coordinate axis for arbitrary vertical stress (e.g., 70 MPa) in Fig. 4. This type of graphical tool, representing the possible range of horizontal stresses without failure when the vertical stress is given, is called as the stress polygon. A stress polygon is divided into three triangular regions corresponding to the three stress regimes of NF, SSF, and RF.

A few key characteristics of the stress polygon are worth to be noted for its expression and application. First, the stress state for a point lying inside a stress polygon is most stable at the center where the three principal stress magnitudes equalize (S_v = S_HMax = Sh_min = 70 MPa in Fig. 4), whereas it becomes more unstable (closer to failure) as the point moves farther from the center to the upper and left boundaries of the polygon (1, 2, and 3 in Fig. 4). In terms of the size of a stress polygon, as one may expect from the M-C criterion defining the boundaries, the size increases or decreases when the M-C criterion parameters such as the cohesion, friction coefficient, and applied pore pressure, are changed. For example, the size of a stress polygon decreases as the pore pressure increases, and increases as the friction coefficient increases. However, because a stress polygon represents the possible range of horizontal stresses for a fixed vertical stress, the region of interest is defined by only one vertical stress. If the targeted rock volume has considerable depth range of its own, it is difficult to determine a possible stress state corresponding to the entire depth range and thus varying vertical stress on a single stress polygon. In this case, a normalized stress polygon can be used by dividing the minimum and maximum horizontal stress magnitude axes by the vertical stress magnitude [14]. Fig. 5(a) shows an example of a normalized stress polygon and how its size varies for each different friction coefficient.

Stress polygon indicates a range of possible stress states using easily assumed parameters, such as rock mass density, depth, and friction coefficient, without prior information such as in situ stress measurement tests. Therefore, a stress polygon is used as a basis for integrating the stress indicating information and determining the best estimated stress model. Fig. 5(b) shows a specific case of determining the stress model by comparing the stress indicators such as borehole breakout observations, stress inversion of seismicity focal mechanisms, and hydraulic shearing by fluid injection; although each of these indicators alone cannot determine the stress state, the best estimated stress model can be determined by comparing the stress indicators and finding the stress condtion that best fulfills and explains them on a normalized stress polygon [14].

2.2.2 Re-activation of Existing Discontinuities

The fracture zone behavior depends on the applied pore pressure in the discontinuity plane and the mechanical properties. On a discontinuous surface, closing and opening phenomena may occur owing to shear displacement and vertical displacement; additionally, shear dilation occurring along the rough surface may be induced. The processes that cause displacement include large-scale plate movement, excavation, thermal stress changes, and seismic and glacial processes. This process eventually results in changes in the magnitudes of the shear and normal stresses applied to the discontinuity. The Mohr circle was used to visually characterize the stress state for a discontinuous surface. Mohr circles can be used to examine fracture formation and reactivation potential [16,17]. The 3D-Mohr circle can be divided into three circles based on the principal stress magnitude. A large-diameter circle representing the maximum and minimum principal stress difference was drawn, and the other two circles were drawn inside the large circle using the medium principal stress magnitude (Fig. 6). Currently, the magnitude of the normal and shear stresses acting on the discontinuous surface can be expressed inside the Mohr circle. The normal and shear stresses can be derived using the slope and direction information of the discontinuous plane, respectively, as shown in the following equation:

where dip and dip-direction information β₁, β₂, and β₃ represent the angles formed by the maximum, minimum, and intermediate principal stresses, respectively, with respect to the normal direction of the discontinuity.

Fig. 7 shows principal stress and direction information that include results from Eq. (3). The direction of principal stress can be determined through the stress regime of NF, RF, and SSF. Using this relationship, the slip tendency (T_s) of a fracture zone can be analyzed by the ratio of the shear stress and the normal stress applied to the discontinuity plane [18].

where σ ′ n is effective stress (σ_n − P_p) considering the pore pressure.

The above formula is applicable when the cohesion (c) of the fracture zone is negligible. Moreover, slippage occurs when the slip tendency T_s is greater than the friction coefficient (μ) of the fracture zone. Hawkes et al. [19] reported the shear tendency, including the dip angle and physical properties of the fracture zone, as follows:

where δ is the dip angle, ϕ_fault is the friction angle of the fracture zone, and P_p is pore pressure.

The re-activation of the existing fracture zone can be characterized using information regarding the direction of the fracture zone and the state of local accompanying stress. In the case of a non-planar fracture zone, the slip tendency can be examined using directional information regarding the fracture zone center or by subdividing the fracture plane and including all directional information at each point [11]. Using the above-mentioned three-dimensional Mohr circle and the equations for slip tendency, it was confirmed that the re-activation of the fracture zone can be investigated by considering all specific factors (groundwater, fracture geometry, and fracture properties).

2.2.3 Deep Depth-induced Stress Failure: Spalling

Before the construction of the highly radioactive disposal repository, the target rock mass was placed in a specific stress state as a result of global-scale tectonic conditions and its own weight. In the general case, the stability of an inhomogeneous rock mass in a low-stress field is focused on the rock block behavior near the tunnel and the growth and re-activation of the fracture. However, in the deep depth range, where high stress is applied, the behavior of the fracture is suppressed by a high-stress magnitude. In this case, spalling failure caused by damage with overpressure or deviation stress may be the main problem in the behavior of the rock around the cavity [20]. In underground cavities, high local and tangential stresses are concentrated near excavation walls. Mechanical spalling occurs at a high stress that exceeds the rock strength yield stress. This causes brittle failure, in which rock fragments fall out with a plate impact parallel to the excavation surface. Therefore, it is important to understand the mechanism for V-shaped and plate-shaped failure at the borehole wall resulting from the construction of a large diameter, deep-depth disposal hole. The Atomic Energy of Canada Limited (AECL) conducted a series of experiments and proposed evaluation models regarding the brittle failure mechanism in deepdepth rock mass [21]. Because of the destructive characteristics of brittle rock, the state of the excavation hole after failure becomes stable, and the depth of fracture and extent of the damage may be the principal factors for the stability analysis rather than the magnitude of the applied stress at the wall. Fractures are generated when the deviation stress at the borehole wall exceeds a certain critical value. The critical value refers to the strength of the rock mass, and is reported as 30‒50% of the uniaxial strength. During the fractures in the critical direction, the brittle failure process resulted in developing a V-shaped notch. As such, it is important to analyze the failure strength of the in-situ rock, which indicates the point at which spalling occurs with crack-initiation stress. Several studies have reported that the stress magnitude when spalling occurs is consistent with the lateral strain onset strength, which indicates that the longitudinal extension strain is measured in a uniaxial compression test in the laboratory [22]. Simultaneously, the depth of the notch in the spalling phenomenon is important for the evaluation of borehole stability. In general, a thin plate rock fragment is formed in the notch-formed area, and the thickness is observed to be as small as 1–5 cm. Martin et al. [21] found that brittle failure occurs when the intrinsic cohesion of the rock mass is diminished while focusing on the parameter m, in the Hoek-Brown equation, which is assumed to be negligible. Based on this assumption, the fracture propagation equation based on the cohesion and maximum and minimum principal stresses is presented as follows:

where, in the case of brittle rock, s usually is in value of s = 1 / 3. where σ₁ and σ₃ are the maximum and minimum principal stresses at the borehole wall, respectively, when a failure occurs.

The spalling notch depth was estimated using the proposed empirical failure equation (Fig. 8(a)). The proposed model predicts the maximum damage depth within a reasonable deviation stress (σ₁ − σ₃) range. Subsequently, the critical state of crack initiation, which is the first step in the spalling process, is defined. Martin et al. [22] derived a spalling-depth estimation equation for a specific excavation radius through eight case studies of induced spalling failure (Fig. 8(b)).

The upper linear correlation equation predicts the depth of failure (R_f), which uses the ratio between the maximum tangential stress and uniaxial compressive strength of the rock mass, such as σ_max = 3σ₁ – σ₃. The brittle fracture yielding point is σ_max / σ_c ≈ 0.4 ± 0.1.

2.2.4 Time-dependent Stress Failure: Creep

To evaluate the long-term process and stability of a repository, it is necessary to evaluate the long-term bedrock behavior. Long-term stability is important, including slowly disturbing stress states between fracture zones that can trigger earthquakes. In addition, the creep phenomenon of the bedrock can be discussed as an example of a representative time-dependent deformation. The creep phenomenon is a state in which the deformation of the rock continues, despite the magnitude of the stress being constant. When a constant stress is applied to the rock for a prolonged duration, the creep behavior can be estimated using a temporal graph of strain [23]. The three stages of creep behavior are expressed by a quadratic function. Each step has a different rate and strain ratio as a response to time (Fig. 9). The primary creep velocity exhibited faster and nonlinear growth as the loading stress increased. The duration of the primary creep is short, usually less than 1 day, and strain recovery occurs along the line ‘P-Q-R’ (as shown in Fig. 9) unless stresses exceeding the primary stage are applied. In the line ‘P-Q-R’, strain recovery consists of elastic recovery with an immediate response and creep recovery with a time-dependent response. In the steady-state creep stage, the strain increases linearly. A secondary creep increase occurred over a relatively prolonged period. The period varied according to the magnitude of the applied load. In the secondary creep stage, strain recovery occurs along the line ‘T-U-V’ (as shown in Fig. 9), and irreversible strain occurs, which indicates permanent deformation. Finally, if the loading condition continued, the creep strain accelerated, and the rock mass failed. The duration of the tertiary creep stage is short and proportional to the n-sequence of time. The step-by-step strain characteristics of a series of creep behaviors can be expressed in terms of time as follows [23]:

where ε_i is the primary elastic strain, ε₁(t) is the timedependent strain, V_t is the stabilized and linearly increasing strain, ε₃(t) is the accelerated creep strain immediately prior to failure. The time-dependent strain, ε₁(t), has been proposed in various forms, and Jaeger et al. [23] summarized the results of several previous studies as follows:

where A, n, and a included in each formula are constants for the dimensional analysis. It is suitable for predicting the strain at a specific time during which transient creep occurs.

3. Numerical Method for Simulation of Hydro-Mechanical Behavior of Rock

To reasonably simulate hydraulic-mechanical rock behavior, the coupled behavior mechanism should be reviewed first. It is necessary to determine the numerical technique to use based on the relative size and geometry of the fractures in the rock. The model should be conducted as either a continuum model or a discontinuous model, depending on the stress state and characteristics of the fracture sets [24]. The continuum approach can be applied when small-scale fractures are present in the rock mass or when the distinct behavior of each block is not a significant factor [25]. As a continuum model application method, rock media with many small-scale fractures are represented as an equivalent continuum model using a porous medium with low permeability. An equivalent continuum model can contain largescale fracture zones using a ubiquitous fracture model [21]. The discontinuous method is suitable when the number of fractures is significant, and fractures are a major factor for analyzing the rock mass behavior. After suitable domain construction based on the fracture geometry and characteristics, a numerical technique for analyzing hydro-mechanical coupled behavior should be adopted and achieved. As a method for coupled analysis, a monolithic approach derives the results through one governing equation in a single code. Moreover, a sequential analysis method uses two or more principles in an interactive-related process [26]. Therefore, the construction of a numerical model for coupled behavior constitutes a solution by coupling the differential equations representing each physical mechanism. Thus, in the hydromechanically coupled process, the effective stress in the mechanical matrix is calculated using the derived fluid pressure. The porosity and transmissivity are updated as a function of the stress or strain. The following section introduces commercial programs that are generally used for continuum and discontinuity analyses. The coupled analysis process and the characteristics of each numerical technique were reviewed.

3.1 Hydro-Mechanical Modeling in Distinct Element Method

Among the various numerical methods that can be applied to geotechnical engineering processes, the discrete element method (DEM) is a representative numerical modeling technique that directly simulates the behavior of fractures in a rock mass. As mentioned in Section 2.2.2, the rock mass contains various discontinuities, including fractures, joints, and fracture zones, created by geological evolution. These discontinuities influence the hydraulicmechanical behavior of the rock mass because of the low strength and high permeability compared to intact rock. The distinct element method is distinguished from a continuous analysis based on the existence of an interface or contact between discrete blocks. The entire domain system will be made up from the discrete block and the behavior of the domain will be analyzed from the motion and interaction between blocks. The assigned constitutive law in the contact and interface between the blocks describes the model behavior [25]. Representative commercial programs based on the distinct element method include Itasca’s UDEC and 3DEC software [27]. In addition, there are many studies on hydro-mechanical analysis [28]. The principal points of the distinct element analysis are the geometric data and the properties of the fracture interface. The geometric data are related to the implementation of the fracture or fracture zone information mentioned in Section 2.2.2, and the physical properties correspond to the stiffness and strength values between blocks given to the constitutive equations. The UDEC and 3DEC analyze the disturbance phenomenon depending on the geometrical and physical properties of the discontinuous interface in the block system. The block motion is invoked by the body force and changeable contact force to derive the new nodal velocities and displacements.

The motion may be resolved into normal and shear displacement components along the direction of the applied stress. The contact stiffness must be considered to calculate the displacement between blocks in the elastic deformation stage.

where Δσ_n and Δτ_s are the effective normal and shear stress increments, k_n and k_s are the normal and shear stiffness of the discontinuities, Δu_n and Δu_s are the normal and shear displacements, respectively. The dilation of apertures occurs only when the discontinuities are at slip, considering the roughness characteristics. An irreversible increase in the discontinuity transmissivity can be induced by dilation slip. Generally, the dilation angle (ϕ_d), which represents the strength and roughness characteristics of the discontinuities, can be adopted for the simulation of the aperture opening.

where Δ u n s is the normal displacement according to the dilation induced by shear displacement. Discontinuity interface hydraulic characteristics are calculated by aperture changes along with the motion of the blocks. Therefore, the aperture size, which was modified for each defined mechanical interpretation cycle, was applied to the hydraulic analysis cycle. The aperture is occasionally defined as a separate mechanical and hydraulic aperture. It is difficult to measure the microscopic value of a small hydraulic aperture compared with a mechanical aperture [29]. Therefore, the hydraulic aperture is a function of the mechanical aperture and is defined as the average distance of the contact interface between the blocks. In the 3DEC simulation, the fluid flow behavior in the discontinuity plane was examined by calculating the hydraulic aperture change with respect to the effective stress (Fig. 10).

where, e_h0 is the initial aperture at zero normal stress and Δe_n is the aperture change according to the normal displacement. In Fig. 10, e_res is the residual aperture size, which is the minimum aperture during the simulation. This size has no effect on the transmissivity modification. In addition, the maximum hydraulic aperture size could be inputted, 3DEC assumes that the hydraulic aperture and the mechanical aperture are identical until the aperture size reaches the maximum hydraulic aperture size owing to normal displacement. The mechanical aperture size that is increased by normal displacement can be larger than the hydraulic aperture size. However, the hydraulic aperture size is adopted in the cubic law, which calculates the flow rate in the discontinuities. Therefore, the exceeded mechanical aperture was not considered a transmissivity change.

Therefore, in the 3DEC simulation, the flow rate between the block interface is calculated using the above equation, considering the hydraulic aperture size changes due to the external stress effect. The fluid pressure and volume change in the flow knot, which indicate the fluid flow path between the blocks, depend upon the derived flow rate. Under this series of processes, a hydraulic-mechanical coupled analysis was sequentially performed. As in the equations listed above, a discontinuous aperture representing a fluid flow path in the rock mass influences the hydraulic properties of the entire rock mass. In the case of most crystalline rock simulations, the analysis was performed because the medium of the intact rock was considered unsaturated.

3.2 Hydro-Mechanical Modeling by the Finite Element Method

The finite element method (FEM) is widely used in most engineering fields, including rock engineering. It is commonly used to analyze the applied loads and deformations in the rock mass. FEM requires a modeling process called discretization, which divides the analysis object into small units (finite elements) and configures a system by connecting those elements with nodes or boundaries [25]. Because FEM can apply various boundary conditions and solve problems with different equations for each element, it is convenient to model a medium composed of different phases, as in hydraulic-mechanical analysis. In the FEM analysis, the fractured rock mass was hydromechanically simulated by composing the rock mass and fractures into porous-elastic domains. It is based on a coupling analysis of the flow rate and fluid pressure by examining the permeability or storativity, which varies depending on the stress. COMSOL Multiphysics is a widely used FEM program that applies the constitutive equations for fluid flow in a rock mass as follows.

where ρ_w is the density of water, ϕ is the porosity (ϕ_f and ϕ_m are the porosities of the fracture and rock mass, respectively), t is time, q is the flow rate, U_m is the source term for the mass, ε_v is the volumetric strain, and α is the Biot constant. The flow rate of the fluid was determined using the Darcy equation.

In COMSOL, the permeabilities of the fracture and rock mass can be distinguished as k_f and k_m, respectively. The storativity, which is the water storage capacity in the fracture and rock mass, can be expressed by the following equation:

where S is the storativity (S_f and S_m are the storativities of the fracture and the rock mass, respectively). The equations above simulate the fluid flow in the rock mass, as shown in the following equation:

The fluid flow in the fracture can be simulated by adding a fracture flow interface inside the domain, which is defined as an impermeable or porous rock mass. The fluid flow characteristics are determined based on properties such as porosity, permeability, and aperture of the fractures. If the flow characteristics follow the cubic law, the flow rate per unit length (Q) of the fracture is expressed through the aperture (e), as shown below:

where f_f is a parameter related to the fracture roughness and has the following relationship with the fracture permeability (k_f) and aperture (e):

As in Eq. (16), considering the storativity of the rock mass, the governing equation considering the storativity of the fracture is expressed in Eq. (20), where the storativity of fracture (S_f) exhibits the relationship shown in Eq. (21).

where C_w is the compressibility of water, and K_n is the normal stiffness of the fracture.

As mentioned above, there is a case of analyzing the hydro-mechanical behsavior by FEM analysis by constructing rock masses and fractures as different porous media or fluid conduit domains (Fig. 11).

4. Case Studies to Review the Hydro- Mechanical Modeling Method

4.1 Case Study 1: Boundary Condition for the Long-term Evolution of Tectonic Structure

4.1.1 Critical Issue of the Model

Numerical modeling for medium/large-scale rock masses is generally performed by applying an in situ stress state (direction and magnitude) as a model boundary condition and subsequently simulating the rock mass behavior within the model domain. Hakala et al. [31] conducted a study to analyze the spatial variability of the strata and changes in the distribution of stress fields at the Forsmark site, which is being discussed as a possible disposal site in Sweden, through a regional-scale model applying the movement of the thrust fault. The geological structure of the brittle deformation zone at the Forsmark site was first modelled, and a regional-scale model was constructed considering the geometric information and mechanical properties of the brittle deformation zone. The thrust fault (compressible) velocity and in-situ stress state were applied as boundary conditions.

4.1.2 Analysis Method

3DEC, a three-dimensional distinct element method program, was used to characterize the geological information from 110 brittle deformation zones in the regional model (15 × 11 × 2.1 km) for the Forsmark site (Fig. 12). The analysis consisted of the following two phases.

4.1.3 Key Results

Considering the deformation zone’s geometry, the shear displacement magnitude was greater for the undulating surface than the planar surface. In addition, the distribution range is more significant for the former. The magnitude of the shear displacement not only varied depending on the depth, but was also affected by the horizontal stress acting on the model and the friction angle of the deformation zone. In particular, the difference in the displacement according to the magnitude and direction of the stress increased in the brittle deformation zone with low shear strength. The mechanical properties of the bedrock other than the brittle deformation zone did not significantly affect the variation in the stress state (Fig. 13).

The modeling of Phase 1 has limitations regarding its stress state estimates. The stress state may vary by the depth and may be disturbed by the deformation zone shear displacement by applying the in situ stress measured in the field as a boundary condition. Therefore, an analysis technique was developed for Phase 2, which involved thrust fault velocity as a boundary condition. The analysis results for Phase 2 differ from those of Phase 1, specifically regarding the variation in depth of the stress magnitude and direction distribution. As the shear displacement of the brittle deformation zone increases, the deviation in the magnitude and direction of the principal stress increases. However, as the depth increased (to less than 550 m), the magnitude and direction of the principal stress converged within a specific range. As a result, the hydro-mechanical model with thrust fault velocity as the boundary condition was reasonable for long-term behavior analysis because it can effectively simulate the temporal and spatial variation of the site stress.

4.2 Case Study 2: Upscaling Problem in Fractured Rock Simulation

4.2.1 Critical Issue of the Model

For analyzing the coupled behavior of rock mass, particularly fractured rock mass, upscaling is essential for expanding the analysis domain to a large-scale disposal site and evaluating the site performance. This is because fractures, although small, can influence the hydro-mechanical characteristics of the rock mass. The focus of studies on upscaling in the hydro-mechanical coupled behavior model is to determine the representativeness of the observable coupled behavior pattern, which can be used to evaluate the disposal site’s overall performance. In this regard, calculating the equivalent continuum properties by determining the representative size of the fracture system that can be applied to the disposal site performance evaluation is a priority. The representative size is generally referred to as the representative element volume (REV), and it is accompanied by an operation that determines its existence. The REV is determined by the size at which the hydraulic and mechanical properties of the fracture system converge uniformly, regardless of the model size. Consequently, upscaling in the fractured rock mass was performed using the equivalent continuum approach, which assumes that a greater rock volume than the REV is an equivalent continuum medium [32]. The results of previous studies on upscaling may help to construct a reasonable numerical model for the performance evaluation of natural barriers and verify the validity of the model. Benchmark Test 2 of DECOVALEX III addressed upscaling issues related to a hydro-mechanical analysis by analyzing groundwater flow in the case where the disposal system is located within a fractured rock mass [33].

4.2.2 Analysis Method

The benchmark problem applied boundary conditions and geological characteristics based on data surveyed near Sellafield, United Kingdom. The analysis techniques for the benchmark problem were reviewed by comparing the results of studies conducted in various countries such as France, Japan, and Sweden. The upscaling processes applied by each research team were as follows.

Each analysis technique constituted a two-dimensional or three-dimensional model depending on computer performance and software. The domain sizes for determining the REV in each analysis technique varied. Most of the research teams used the distinct element method program (3DEC/UDEC) on a small scale. After scale expansion, the analysis was performed using the finite difference method (FLAC/TOUGH2) or the finite element method (ROCMAS) program. The crack tensor theory was mainly applied using discontinuum analysis. The upscaling process was performed by assuming the rock mass as a porous medium in the continuum analysis. For example, KTH in Sweden and INERIS in France performed an upscaling process for the hydraulic properties using UDEC and 3DEC [34,35]. Min et al. [34] performed analyses considering the boundary effect by extracting the size of 0.25 to 10 m from the two-dimensional 300 × 300 m domain. The hydraulic aperture of the model was determined based on the experimental results, and the permeability of the fracture system, which could be direction-dependent, was analyzed by comparing the permeability tensors of the rotated models. Thoraval and Renaud [35] performed analyses to determine the REV by varying the model size from 2 m to 5 m using 3DEC. The results analyzed the effect of geometrical properties, such as fracture length, regarding the REV determination and the fracture system’s permeability. The two research teams derived representative equivalent properties by examining the compliance tensor for the mechanical properties of the fracture system using the twodimensional and three-dimensional models constructed, respectively.

4.2.3 Key Results

The geometrical properties of the fracture (length, aperture, density) determined by statistical methods could significantly affect the REV. Therefore, in constructing the DFN model, research on the applicability of statistical methods based on field data should be conducted. However, the correlation between the size (length) and aperture of the fracture assumed during the modeling process has not yet been verified. The uncertainty caused by this assumption is still present. In addition, the three-dimensional analysis was limited in domain size owing to the computational cost [35]. Although the two-dimensional analysis has limitations regarding model reliability, it was possible to examine the REV by constructing domains of various sizes compared with the three-dimensional analysis (Fig. 14).

In the hydraulic-mechanical coupled model, the maximum mechanical aperture and minimum hydraulic aperture had a significant role. Therefore, it is necessary to analyze the changes in the aperture and fluid flow path according to the applied stress and their correlation through laboratory tests. In particular, the differential stress (maximum principal stress–minimum principal stress) was emphasized with respect to the increase in rock mass permeability and fracture connectivity. The fluid flow channels are expected to be identified by reflecting these correlations in the upscaled model.

The hydraulic behavior of the upscaled model was not significantly affected by thermal factors. It was primarily dependent upon the external stress and pore pressure (fluid pressure) within the fracture aperture. In addition, in the upscaled model, the fracture aperture was maintained at the residual size because high in situ stresses were applied according to the depth. Therefore, the variation in the hydraulic properties due to the applied stress was marginal. These results indicate that the hydraulic properties of the rock mass measured at an appropriate depth can be represented within a depth interval.

As mentioned above, the correlation between the applied stress and hydraulic aperture of the fracture may have a significant role in constructing a hydro-mechanical coupled model. Therefore, the appropriate fracture density that can represent the hydraulic properties of a rock mass must be established by analyzing the minimum hydraulic aperture so that fluid flow, such as solute transport, in a fractured rock mass, can be reasonably simulated.

4.3 Case Study 3: Stress Field Distribution and Disturbance Problem

4.3.1 Critical Issue of the Model

The rock mechanics site descriptive model (RMSDM) is a geological structure model constructed to evaluate the mechanical stability of the disposal site or to couple the hydro-mechanical characteristics in consideration of the in-situ stress state variation with depth. The site descriptive model aims to evaluate the current state and long-term geological evolution (past-present and present-future) of the site by appropriately simulating the rock mass properties and topographical characteristics observed in field tests. Related studies have primarily discussed the model construction methodology and have suggested practical evaluation and application methods for the data acquired during the site investigation stage [36]. In these studies, the site descriptive model has a role in empirical, theoretical, and numerical simulations of the properties of the rock mass and fracture zones distributed on the site and in evaluating the in-situ stress distribution considering the fracture zone. Establishing a descriptive site model is essential because it can be used to predict and review the geotechnical parameters required for the design and stability evaluation of a disposal site.

4.3.2 Analysis Method

The benchmark model for constructing the site descriptive model was the Äspö Hard Rock Laboratory in Sweden. A study regarding rock mechanical evaluation and the establishment of a geological structure model was conducted on a three-step scale for the model site (Fig. 15(a)). The construction of the site descriptive model could be developed by acquiring sufficient field data through site investigation and applying these data. In addition to constructing the model, it is necessary to establish an application procedure for other research sites in the future by supporting a systematic explanation of the applied theory and field data analysis method. The priority in the model establishment process is to consider the rock mechanical factors required for site design and stability evaluation. In particular, the model should adequately represent in situ stresses and mechanical properties of the rock mass and fractures. Fig. 15(b) shows the input data and flowchart for the rock mechanics site descriptive model presented by SKB, which is an essential component of the stress field evaluation process using the site characterization model.

Determining the properties input to the model is important to ensure the representativeness of the site descriptive model. Properties related to site characterization can be determined through the selection and analysis of input data, which are primarily limited to field data provided at the surface or borehole scale [36]. Therefore, the primary consideration in site characterization is the applicable range of rock-mass properties determined from field data. Methods regarding rock mass classification, such as the Q-system and RMR analyzed in borehole surveys, were widely used to determine rock mass properties. Estimating the representative properties using the kriging technique to obtain a three-dimensional statistical distribution was used to extrapolate the information measured in the borehole into the model [37].

Because the stress state of the site can be affected locally or regionally by a large-scale fracture zone. It is also likely to be changed by the weight of the rock mass distributed at the site and the resulting geostatic stresses. Therefore, it is necessary to understand how multiple fractures of various sizes interact with the stress field acting on-site. Consequently, a method for constructing a numerical model with suitable boundary conditions should be discussed to determine the stress state for geological change simulations. In addition, by performing a sensitivity analysis on the model variables, it is possible to estimate the model’s stress field variation and examine alternative models (changes in the geometrical properties of the fracture zone, mechanical properties, and loading conditions). A site descriptive model was constructed using 3DEC for the Äspö Hard Rock Laboratory site (Fig. 16), and the model’s applicability was reported for the Forsmark site in Sweden [39]. As a model construction step, the fracture zone was first simulated on a regional scale centered on the site. The principal stresses that could be induced by the geological structure were analyzed. Next, the rock mass behavior of the site was evaluated by applying the stress boundary condition acting on the local scale model, which was included as a part of the regional scale model. Additionally, a sensitivity analysis of the stress field variation to the input properties of the rock mass and fracture zone, and the geometry of the brittle deformation zone was performed [40].

4.3.3 Key Results

When inputting properties into the site descriptive model, it is essential to understand and evaluate geological evolution in terms of rock mechanics. The method for determining rock mass properties can be validated through the geostatistical theory of the spatial distribution of properties. In addition, it was confirmed that applying an empirical index, such as the Q-system or RMR, to determine the properties of the site descriptive model is only useful at the corresponding site. Therefore, when determining the site descriptive model properties, methods based on empirical relationships can be applied only to geologically similar sites. In other words, both empirical and theoretical approaches for determining the properties of site descriptive models contain uncertainties; however, these uncertainties can be reduced by targeting a specific site having a fresh rock mass. A disposal site with optimal geological conditions should be selected based on these considerations.

In the stress field analysis using the site descriptive model, the size of the analysis model domain and the geometry of the simulated fracture zone had a more significant effect on the input mechanical properties. The comparison results of the analyses indicated that the reliability of the insitu stress measurement and borehole investigation results should be supported to ensure the validity of the model. In addition, the stress state analyzed in the regional-scale model is meaningful for determining the boundary conditions required for the local-scale analysis within the model. However, it should be noted that the presence of large-scale fracture zones in the model or the inhomogeneity caused by the intersection of rock blocks may result in significant variability in the in situ stresses at local locations.

5. Conclusions

For the deep geological disposal of high-level waste, studies evaluating the performance of natural barriers that must be stably maintained during long-term geological evolution are essential. Studies examining prediction scenarios for future disposal sites and applying geological evolutionary events to models that simulate disposal sites are currently being conducted. Such modeling techniques are essential because they can predict the actual behavior of natural barriers caused by evolutionary events in the geological environment by considering the characteristics of the target disposal site. The geological evolution of the natural barriers was evaluated using hydromechanically coupled models. Most models have been constructed considering the possibility of fault re-activation and stress conditions in response to the evolving geological environment, such as glacial loads and earthquakes. Therefore, it is essential to construct an optimized hydro-mechanical model to reliably evaluate the interaction and influence of geological properties. In this study, an optimal modeling technique that can simulate the evolution of the geological environment and structure in Korea was reviewed by examining the coupled modeling techniques developed and utilized by leading countries in nuclear research, such as Sweden, Switzerland, and Finland.

First, processes corresponding to long-term geological changes were analyzed and identified by reviewing the defined FEP and suggested by foreign organizations. Next, the hydro-mechanical behavior mechanisms involved in each process were summarized. In addition, factors related to the long-term behavior of natural barriers were identified, and significant evolutionary mechanisms of natural barriers were analyzed through the selected FEP. Long-term behaviors such as stress states, fault re-activation, spalling, and creep in the rock mass were considered the principal mechanisms.

Consequently, various analytical techniques have been reviewed to reflect these geotechnical mechanisms. In addition, case studies were conducted regarding the stability of the disposal system in underground research laboratories or disposal candidate sites using each analysis technique. A new approach to boundary conditions applied to numerical models was investigated by analyzing the stress field distribution that evolves as a result of the movement of the strata. In addition, studies on upscaling and equivalent continuum approaches in fractured rock masses were reviewed to discuss the extension of observable fractures’ hydraulic properties to regional site scales. Finally, studies on the utilization of field data and reasonable stress field analysis were summarized by examining the construction method of the site descriptive model.

This study reviewed modeling techniques that can effectively examine the long-term performance of natural barriers. It is helpful as a preliminary study to develop an optimal modeling technique that can simulate the in situ stress distribution and geological evolution scenarios in Korea. Based on the results reviewed in this study, future studies are required for developing technologies related to the long-term stability of natural barriers.

The remaining issues are summarized as follows. These include the construction of a hydro-mechanical site descriptive model and analyzing stress field variation using field data. Specifically, this includes:

In addition, applying long-term evolution scenarios in coupled models and analyzing site stability for disposal candidates requires additional consideration.

The techniques and case studies of the coupled model summarized in this study are expected to be used as preliminary data for selecting candidate sites for disposal and evaluating the site stability in the future.