A Multimodel Framework for Quantifying Flow and Advective Transport Controlled by Earthquake-Induced Canister Failures in a Reference Case for Radioactive Waste Geological Disposal

Abstract

Characterizing flow and transport for earthquake-induced shear canister failure is critical for the performance and safety assessment of radioactive waste geological disposal. The study presents a modeling framework that integrates multiple models to account for fractures produced by shear displacements, evaluate canister failures, and simulate flow and advective transport in a conceptual repository site based on a selected reference case in an offshore island in western Taiwan. The typical KBS-3 disposal concept associated with 500 realizations of the shear-induced fracture properties is employed to quantify the uncertainty of flow and advective transport in the geological disposal site. The radionuclides in canisters are assumed to migrate through the shear-induced fractures surrounding the deposition holes. The results from 500 realizations show that two types of fractures produce a high potential to destroy canisters induced by the shear displacements. The earliest canister failure time influenced by possible shear movements is 0.23 million years for the reference case. The modeling framework identifies five canisters and the associated shear-induced fractures for flow and advective transport simulations. Based on the results of the density-dependent flow fields, the particle tracking algorithm enables the calculations of flow and transport parameters, including equivalent initial flux, equivalent flow rate, path length, travel time, and flow-related transport resistance for the identified five canisters. These parameters are critical for the performance and safety assessments of buffer erosion and canister corrosion near the disposal repository and the far field of the radioactive waste disposal site.

1. Introduction

High-level radioactive waste and spent nuclear fuel are critical issues worldwide. Many investigations have shown that deep geological final disposal is currently the best method for managing radioactive waste. Svensk Kärnbränslehantering AB (SKB) proposed the popular KBS-3 disposal concept based on the primary purposes of isolation, containment, and retardation [1]. The multi-barrier system in the KBS-3 disposal concept comprises an engineering barrier system and a natural barrier. The natural barrier is the host rock, while the engineering barrier system consists of the canister, buffer, and backfill. The current design is based on the original concept presented in the KBS-3 study [2]. The canister has three safety functions for containment in the KBS-3 disposal concept, including (1) providing a corrosion barrier, (2) withstanding isostatic loads, and (3) withstanding shear loads [3]. Canister integrity may generally be threatened chemically or mechanically [4]. For mechanical load, in the Nordic countries, the canister must withstand an isostatic load of 45–50 MPa, which is the sum of the maximum swelling pressure and maximum groundwater pressure caused by ice sheets with 3–4 km width in the glaciation period. In addition, the canister should remain intact after a 5 cm shear displacement caused by intersecting fractures [5]. However, the isostatic load is more critical than the shear load in Sweden and Finland since the magnitude and frequency of earthquakes are relatively low [3,4,5].

The territory of Taiwan is located on the boundary of two tectonic plates (the Eurasian Plate and the Philippine Sea Plate). It belongs to the Pacific Ring of Fire, so the geological environment, faulting, and seismicity are more frequent than in other countries. Earthquakes frequently occur in Taiwan, making them the critical mechanisms for destroying the canister induced by the shear movement of fractures intersecting the deposition hole (DH). In general, the shear displacement of fractures induced by earthquakes is due to the deformation of host rocks and the wide-spreading seismic waves. For the current design requirement, canisters should be able to withstand 5 cm shear displacement without failing. Fractures accumulating shear displacement larger than 5 cm may yield intersected canisters [3,4,5]. To evaluate the risks of canister failure caused by shear displacement, SKB uses a numerical program called 3DEC, based on the three-dimensional (3D) distinct element method, to calculate the possible shear displacement on fractures induced by earthquakes [3,6,7,8]. The simulations indicated that the shear displacement of fractures with various sizes could be evaluated and identified.

If an earthquake-induced shearing fracture intersects the repository and damages the canister, groundwater will come in contact with the radionuclides inside. The radionuclides will determine the rate at which they dissolve in the groundwater, migrate across the near field (repository area) and far-field (host rock), and potentially reach the biosphere [5,9]. For conservative purposes, most studies assume that the buffer loses all of its mass in the near field and that the radionuclides directly migrate to the biosphere based on the calculations of release paths in the flow and transport modeling. The assessments rely on determining the smallest trace length or highest transport velocity in the far field [10,11,12,13]. This release path modeling might be too simplistic to serve as the fundamental parameters for the following safety assessment in the far field. In addition, such flow and advective transport models neglect the earthquake-induced shearing fractures and the detailed flow mechanisms of this type of fracture. Therefore, the safety assessment results for withstanding shear load may be overly conservative and impractical for the cases in Taiwan.

According to the “Nuclear Materials and Radioactive Waste Management Act”, Taiwan Power Company (TPC), the radioactive waste producer, is responsible for managing, storing, and finally disposing of all radioactive materials, and one of them is radioactive waste [14]. Recently, TPC has adopted the KBS-3 disposal concept to develop geological disposal capability and technology and selected the crystalline rock as one of the potential host rocks. To improve the assessment technologies, TPC generated a reference case that consists of the geological unit and thermal, hydraulic, mechanical, and chemical (THMC) properties based on the experience from the Japanese H12 report [9] and the site-descriptive model (SDM) concept from Swedish SKB [15]. Additionally, TPC designed a conceptual repository based on the reference design of the KBS-3 concept [3,4,5]. These reference case data and the conceptual repository are open to researchers with the same basis of information to iterate, improve, and modify the available technologies for assessing radioactive waste geological disposal sites in Taiwan [16].

The study aims to develop a modeling framework considering earthquake-induced shearing fractures for groundwater flow and advective transport simulations. A reference case in Taiwan is selected for presenting the developed framework. The study first focuses on calculating the geometry and locations of the high probability earthquake-induced shearing fractures and their associated magnitudes of shear displacements based on the discrete fracture network dataset (DFN recipe) and a synthetic reference case in Taiwan. The earliest canister failure time due to shear movements is then calculated based on the specified failure criteria. Finally, the identified earthquake-induced shearing fractures will be employed in the density-dependent groundwater flow and advective transport processes. The site-specific flow and transport parameters for performance and safety assessments are then estimated based on the failure canisters induced by the earthquakes.

2. Materials and Methods

2.1. Workflow for the Study

There are three major tasks in this study. The first task involves calculating fracture shear displacement induced by earthquakes. We use earthquake source data from Taiwan and the DFN recipes to calculate the shear displacement of different fracture sets. The second task is to evaluate the potential failure of canisters due to the shear displacement of fractures. Based on the results of the first task and the geometry of the conceptual repository, we could calculate the time and location of canister failure due to shear displacement. In addition, information on specific shear fractures, including location, size, and attitude, is employed for groundwater flow and transport modeling. The last task involves simulating the groundwater flow field based on the hydrogeological conceptual model. Specifically, we calculate the critical flow and transport parameters, or named as performance measures in the KBS-3 disposal concept [3,5,10,11,12,13], including U, Q_eq, L, t_w, and F, for the DH destroyed by the shear displacement.

There are several steps involved in the workflow. Figure 1 shows the main steps for numerical modeling. First, we examine input data, including the geological conceptual model, earthquake simulation model, and layout of the conceptual repository, and generate the associated objects for other steps. Next, we create the numerical domains for shear displacement and potential canister failure and conduct groundwater flow and advective transport simulations. The simulation domain of shear displacement is created to calculate the maximum and cumulative displacements of five fracture sets in the R# region. The associated fracture displacements are then used to calculate potential canister failure caused by shear displacement. In this step, the geometry of the conceptual model is considered to determine the types of shearing fractures and their locations. The critical shearing fractures are then identified and incorporated into the numerical domain for the flow model. The flow model comprises the hydrogeological conceptual model and an example of DFN realization based on the DFN recipes. The example DFN comprises millions of fractures. The critical shearing fractures are incorporated into the example DFN, and the whole fracture systems (i.e., the critical shearing fractures and an example DFN) are to be the foundation for flow and advective simulations. Finally, the critical shearing fractures causing potential canister failure are integrated into the DFN and used in the following flow simulations. In this study, pressure and salinity distributions are obtained in the flow simulations. The particle tracking algorithm is employed for the release path simulations from the conceptual repository to the biosphere. We consider the Q1 path based on the KBS-3 disposal concept and release particles at the intersection between fractures and DHs. Specifically, we consider the Q1 paths released from potential failure canisters to the biosphere and calculate the performance measures for those canisters destroyed by shear displacements. The trajectory of the Q1 path and the associated performance measures are demonstrated for the safety assessment.

2.2. Conceptual Repository and the Reference Case for the Study

Kinmen Island is one of the offshore islands selected for developing technologies for geological investigations. The site is also a simple SDM based on limited information to improve performance and safety assessment technologies. The SNFD2017 report submitted by TPC demonstrated that Taiwan’s geological conditions are suitable for geological disposal [16]. In this study, we generate a conceptual repository and use the same geological data as the reference case shown in the SNFD2017 report for simulation purposes. In the following sections, we present the key information for calculating the shear displacement induced by earthquakes, the potential failure canister destroyed by shearing fracture, and the hydrogeological conceptual model for flow and advective transport simulations in this study.

2.2.1. Layout and Hydraulic Properties of the Conceptual Repository

Figure 2 shows a reference design layout for a conceptual repository with a depth of 500 m. Figure 2a presents the top view of the conceptual repository composed of the main tunnel (MT), disposal tunnel (DT), DH, and excavation damaged zone (EDZ) beneath MT and DT. The design concept is based on the KBS-3 system, an underground geological disposal concept developed by SKB over several decades [1,2,3,4,5]. In the study, the MT has a height of 7 m and a width of 10 m (Figure 2b). DH has a diameter of 1.75 m and a height of 8.155 m. DT has a height of 4.8 m and a width of 4.2 m (Figure 2c) [16]. The first DH position is at least 20.6 m from the entrance of the DT, and the last DH position is located 10 m from the end of the DT. There are two panels spaced 150 m apart, each containing 52 DTs. The DTs in the western panel have a length of 250 m and contain 25 DHs, while those in the eastern panel have a length of 300 m and a capacity of 30 DHs. The center-to-center spacing for the DTs is 40 m, and the center-to-center spacing for the DHs is 9 m. In summary, the conceptual repository comprises 2 MTs, and 104 DTs and has a gross capacity of 2860 DHs.

The study incorporates the repository into our mechanical and hydrogeological modeling.

Table 1. The flow properties of components of the conceptual repository in this study.

Units	Material	Hydraulic Conductivity (m/s)	Porosity (–)
MT	Backfill material	1.0 × 10−10	4.0 × 10−1
DT	Backfill material	1.0 × 10−10	4.0 × 10−1
DH	Buffer material	1.0 × 10−12	4.0 × 10−1
EDZ	Granitic gneiss	3.3 × 10−8 *	1.0 × 10−4

* The transmissivity of EDZ is measured to be 1.0 × 10⁻⁸ m²/s [18]. Considering the width of the EDZ, which is approximately 30 cm, we calculate the hydraulic conductivity by dividing the transmissivity by the width.

shows the flow properties of each component of a conceptual repository [16,17,18]. The hydraulic conductivity for MT and DT is 1.0 × 10⁻¹⁰ m/s, representing the flow properties of backfill. The hydraulic conductivity is 1.0 × 10⁻¹² m/s inside DH, representing the property of buffer [17]. According to the formal KBS-3 disposal concept, a canister should be deposited at the center of the DH. However, we exclude the canister and assume that each DH is entirely filled with buffer for simulation purposes. Additionally, previous studies have indicated that the width of the EDZ is approximately 30 cm. The transmissivity ranges from 1.0 × 10⁻⁹ to 1.0 × 10⁻⁷ m²/s [18]. Therefore, the recommended transmissivity for EDZ is 1.0 × 10⁻⁸ m²/s, and the porosity is 1.0 × 10⁻⁴ (–) [18]. Consequently, we calculate the hydraulic conductivity of the EDZ by dividing the transmissivity by the width of the EDZ. Based on this calculation, the hydraulic conductivity of the EDZ is determined to be 3.3 × 10⁻⁸ m/s.

2.2.2. Reference Case and the Relative Properties of the Geological Conceptual Model

The range, geometry, and classification of geological units are the basis of the reference case. Figure 3 shows the geological conceptual model of the reference case. The reference case is primarily composed of granitic rock, with lower fracture intensity in the host rock and higher intensity in the thinner sedimentary layer on the top. The geological conceptual model integrates in situ geological investigations and data from six additional boreholes and their logging. All information derived from field data, in situ hydraulic experiments, and borehole logging data are recognized as the DFN recipe for fracture generations. The generated fractures are for detecting intersections between fractures and repository, shear movement induced by earthquakes, hydrological flow, and transport, among others [16,17,18,19,20,21,22,23]. There are three geological units in the domain, including regolith (R0), host rock (R#), and major water-conductive fault/fracture (F#). The depth of R0 is 70 m based on the geological investigation and six borehole logging data. R0 comprises fractures and sedimentations, and the hydraulic conductivity of sedimentation is specified as 1.0 × 10⁻⁵ m/s.

The fracture information within R0 will be presented in the table list of DFN recipes later. The host rock is also composed of fractures and rock matrix, and the hydraulic conductivity of the rock matrix is 1.0 × 10⁻¹⁰ m/s. The fracture information within R# will also be presented in the table list of DFN recipes. Two major water-conductive structures are inside the domain, i.e., fault (F1) and fracture zone (F2). The attitude of F1 is N64E/70N, and the thickness is 200 m, while the attitude of F2 is N80W/50S, and the thickness is 20 m. The F1 and F2 penetrate the R0 and R# and have the same hydraulic conductivity as 5.0 × 10⁻⁶ m/s. The geological units and their hydraulic properties are listed in

Table 2. The hydraulic properties of hydrogeological units for the conceptual model in this study.

Units	Lithology/Material	Distributions/Attitude and Width	Hydraulic Conductivity (m/s)	Porosity (–)
R0	Regolith	70 m thickness of the topmost of domain	1.0 × 10−5	1.0 × 10−3
R#	Granitic gneiss	-	1.0 × 10−10	5.4 × 10−3
F1	Fault	N64E/70N, 200 m width	5.0 × 10−6	1.0 × 10−2
F2	Fracture zone	N80E/50S, 20 m width	5.0 × 10−6	1.5 × 10−2

.

Table 3. The DFN recipe for R0 and R# (modified from Yu et al., (2021) [22]).

Fracture Domain		R0	R#
Elevation		Depth below surface < 70 m	Depth below surface > 70 m
Fracture clusters (Pole trend, pole plunge, κ, P32,rel)	Cluster 1	(198, 18, 18, 26%)	(65, 17, 20, 15%)
	Cluster 2	(155, 4, 15, 24%)	(344, 38, 18, 24%)
	Cluster 3	(264, 23, 16, 18%)	(281, 29, 16, 30%)
	Cluster 4	(98, 81, 11, 32%)	(174, 22, 17, 10%)
	Cluster 5	-	(175, 75, 19, 21%)
Fracture intensity (P32)		2.4	0.3
Fracture size		Power law: ${\mathrm{k}}_{\mathrm{r}}=2.6$, ${\mathrm{r}}_0=0.1 \mathrm{m},$ ${\mathrm{r}}_{\min }=4.5 \mathrm{m},$ ${\mathrm{r}}_{\max }=564 \mathrm{m}$	Power law:${ \mathrm{k}}_{\mathrm{r}}=2.6,$ ${\mathrm{r}}_0=0.1 \mathrm{m},$ ${\mathrm{r}}_{\min }=4.5 \mathrm{m},$ ${\mathrm{r}}_{\max }=564 \mathrm{m}$
Fracture location		Stationary random (Poisson) process	Stationary random (Poisson) process
Fracture transmissivity (T, m2⁄s)		$\mathrm{T}=1.51\times {10}^{-7}\times \left({\mathrm{L}}^{0.7}\right)$; $\mathrm{L}=\sqrt{\left({\mathsf{\pi }\mathrm{r}}^2\right)}$ $\mathrm{L}$ is the equivalent size (m) of a square fracture.	$\mathrm{T}=3.98\times {10}^{-10}\times \left({\mathrm{L}}^{0.5}\right)$; $\mathrm{L}=\sqrt{\left({\mathsf{\pi }\mathrm{r}}^2\right)}$ $\mathrm{L}$ is the equivalent size (m) of a square fracture.
Fracture Aperture (e, m)		$\mathrm{e}=0.5\sqrt{\mathrm{T}}$	$\mathrm{e}=0.5\sqrt{\mathrm{T}}$

presents the DFN recipe within R0 and R#. We focus on DFN generation, shear displacement calculation, and flow and advective transport simulation in the numerical domain composed of the fractured rock and a conceptual repository in this study. For detailed information on establishing the DFN recipe, one could refer to previous studies for details [16,17,18,19,20,21,22,23]. The fractures are generated using the Fisher distribution. The fracture cluster is defined by pole trend, pole plunge, and concentration. Based on the analysis of borehole fracture orientation, R0 consists of four clusters, while R# consists of five clusters. The fracture intensity (P₃₂) is assumed to be 2.4 for R0 and 0.3 for R#. A power law shaper factor (${\mathrm{k}}_{\mathrm{r}}$) of 2.6 is assigned for the fracture size distribution. The minimum radius location parameter (r₀) is set to 0.1 m, and the minimum fracture radius (r_min) and maximum fracture radius (r_max) are set to 4.5 and 564 m, respectively. The spatial arrangement of the fracture center is assumed to follow a stationary random (Poisson) process. The transmissivity of fracture is determined by the fracture size using an empirical function derived from the Swedish Forsmark site. Additionally, the fracture aperture is determined by the fracture transmissivity using an empirical function obtained from Forsmark due to the lack of site-specific data [15].

2.2.3. Hydrogeological Conceptual Model

Figure 4 presents the hydrogeological conceptual model in this study. The length of the modeling domain is 10.5 and 13.8 km on the x- and y-axis, respectively (in Figure 4a), and the bottom boundary is at a depth of −2000 m based on the geological conceptual model (in Figure 3). The domain boundary is stored in a text file (XY file), and the surface topography is stored in another text file (XYZ file). Based on the geological conceptual model, the reference case is mainly composed of granitic gneiss, a kind of fractured rock. Two main water-conductive structures exist, including the F1 and fracture zone F2 within the area. The geometry of F1 and F2 are stored as STL (Stereo Lithography) files to facilitate the generation of R0 and R#. The topmost 70 m thickness is assigned as R0, while the remaining part is assigned as R#. Figure 4b presents the location and layout of a conceptual repository in the numerical domain. According to the layout of the conceptual repository, there are MTs (in blue color), DTs (in red color), DHs, and EDZs. The geometry and location of four components are stored as four STL files, respectively. Figure 4c shows the conceptual repository in a 3D view, and Figure 4d is the close view of parts of an MT, DTs, and DHs (in green color).

Once the geometry and distribution of the conceptual repository and hydrogeological units are generated, assigning the hydrogeological properties (Table 1 and Table 2) becomes feasible. R0 and R# are fractured rock, i.e., composed of the rock matrix and fractures, so there should be a fracture system within R0 and R# based on the DFN recipe (in Table 3). The more detailed DFN generation results will be demonstrated in the following section.

Based on the geological conceptual model, the numerical domain only includes land. However, a specific salinity of 3.2% is assigned at the lateral and bottom boundary to reveal the character of the offshore island [16,19,20,21]. The density of seawater is specified as 1.03 g/cm³. We assume the temperature is uniform, so there is no buoyancy due to temperature variations. The lateral boundary condition is assigned as a specified hydrostatic pressure. However, the westernmost lateral boundary (i.e., the red dotted line in Figure 4a) is assigned a no-flow boundary condition, as the reference case is located in the eastern part of the island. Thick sedimentary layers are observed between the eastern and western parts of the island, so a no-flow boundary condition is a reasonable assumption for the westernmost boundary of the reference case [16]. The bottom boundary is given as a no-flow boundary. The top boundary is assigned a specified recharge rate with fresh water of 66.8 mm/year [16]. The initial salinity condition is assigned as fresh water for the entire domain, which is an initial guess. The model is run for a period enough to obtain a steady-state flow field, and fresh and saltwater interfaces are available based on the model simulations.

2.3. Canister Failure Induced by Fracture Shear Movement

Fracture shear displacement across a DH is one of the canister failure modes over the safety assessment timescale [3,5]. We estimated the magnitude of shear movement induced by earthquakes based on seismicity data, the layout of a conceptual repository, and a fracture system generated using the DFN recipe from the reference case. If a fracture intersects with the canister and the resulting shearing force exceeds the design requirement for the canister’s shear resistance, its integrity could be compromised and seem a failure [3,5]. In this study, we utilized 3DEC for DFN generation and to evaluate canister failure due to shear movements. We describe the numerical algorithm in greater detail as follows.

2.3.1. Mathematical Formulation for Mechanical Modeling

The Coulomb slip model, which includes calculations of shear and tensile failure as well as shear dilation, is a built-in constitutive model in 3DEC [24]. As for the weak plane that has not failed, the upper limit of tensile stress (${\mathrm{T}}_{\max }$) gives

$${\mathrm{T}}_{\max }=-{\mathrm{TA}}_{\mathrm{c}}$$

(1)

where $\mathrm{T}$ and ${\mathrm{A}}_{\mathrm{c}}$ is the tensile strength and area of the weak plane.

The upper limit of shear stress (${\mathrm{F}}_{\max }^{\mathrm{s}}$) can be expressed as

$${\mathrm{F}}_{\max }^{\mathrm{s}}={\mathrm{cA}}_{\mathrm{c}}+{\mathrm{F}}^{\mathrm{n}}\tan \varnothing$$

(2)

where c and ∅ is the cohesion and friction angle of the weak plane. ${\mathrm{F}}^{\mathrm{n}}$ is the normal force acting on a weak plane.

When a weak plane fails in the 3DEC model, regardless of shear failure or tensile failure, the tensile strength and cohesion can then drop to zero. Equations (1) and (2) can then be rewritten to:

$${\mathrm{T}}_{\max }=0$$

(3)

$${\mathrm{F}}_{\max }^{\mathrm{s}}={\mathrm{F}}^{\mathrm{n}}\tan \varnothing$$

(4)

Strength reduction on the weak plane is similar to the displacement–weakening behavior on fracture. The stress acting on the fracture can be expressed as

$${\mathrm{F}}^{\mathrm{n}}<{\mathrm{T}}_{\max },{ \mathrm{F}}^{\mathrm{n}}=0{ \mathrm{and} \mathrm{F}}_{\mathrm{i}}^{\mathrm{s}}=0$$

(5)

where ${\mathrm{F}}_{\mathrm{I}}^{\mathrm{S}}$ is the shear force on the fracture for tensile failure. As for shear failure, the equations are given by:

$${\mathrm{F}}^{\mathrm{s}}>{\mathrm{F}}_{\max }^{\mathrm{s}},{ \mathrm{F}}_{\mathrm{i}}^{\mathrm{s}}:={\mathrm{F}}_{\mathrm{i}}^{\mathrm{s}}\frac{{\mathrm{F}}_{\max }^{\mathrm{s}}}{{\mathrm{F}}^{\mathrm{s}}}$$

(6)

$${\mathrm{F}}^{\mathrm{s}}=\sqrt{{\mathrm{F}}_{\mathrm{i}}^{\mathrm{s}}{\mathrm{F}}_{\mathrm{i}}^{\mathrm{s}}}$$

(7)

Shear dilation can occur as shear displacement and is accumulated on fractures. When the shear displacement is expressed as:

$$\Delta {\mathrm{U}}^{\mathrm{s}}=\sqrt{\left({\mathrm{U}}_{\mathrm{i}}^{\mathrm{s}}{ \mathrm{U}}_{\mathrm{i}}^{\mathrm{s}}\right)}$$

(8)

The shear dilation is given by

$$\Delta {\mathrm{U}}^{\mathrm{n}}\left(\mathrm{dil}\right)=\Delta {\mathrm{U}}^{\mathrm{s}}\tan \mathsf{\phi }$$

(9)

where $\mathsf{\phi }$ is the dilation angle.

In general, the relationship between shear displacement, shear stress, normal stress, shear dilation, and normal displacement is shown in Figure 5. Before the weak plane fails, the shear displacement is proportional to shear stress. As normal stress increases, the critical shear displacement and shear stress required for failure increase while the shear dilation decreases.

2.3.2. DFN Generation for Canister Failure Calculation

Figure 6 illustrates the fracture domain of the DFN model, where fracture centers can be generated within a volume of 1770 m × 1544 m × 500 m. Figure 6 also shows an example of the DFN realization based on the DFN recipes provided in Table 3 and the layout of the repository. The dimensions refer to the repository layout, including a 250 m margin for the maximum fracture radius setting. The repository model is located in the center of the fracture domain. In the study, the minimum fracture radius was set to 2.88 m, and the size is equivalent to the tunnel radius. The shape of the fractures is assumed to be an infinitely thin and circular disc. Furthermore, we extended the fracture outside the domain instead of truncating and recalculating the equivalent radius to maintain a consistent radius for those fractures near domain boundaries.

2.3.3. Deposition Hole Rejection Criteria

The geometrical criteria for rejecting intersected deposition holes developed by SKB were considered to represent the failure canisters caused by intersecting fractures. The objective of rejection criteria is to avoid deposition holes with large fractures intersecting them because such fractures may be subject to movement during earthquakes with subsequent damage to the canister in a deposition hole. Once the geometrical criteria reject a deposition hole, there will be no canister deposited into it. Figure 7 shows the conceptual profile of the deposition tunnel and the associated deposition holes and canisters. During the repository design stage, the space of the host rock must be assessed by considering the number of intersected deposition holes. This study uses a 3DEC numerical model to assess the deposition holes possibly intersected by fractures [24]. The fracture is assumed to be an extremely thin circular plane in the assessment, and the intersected logic is assumed to be within the timescale of the safety assessment, which is one million years. Even if shear displacement occurs on the fracture, the radius of the fracture will not increase [25,26,27].

In the study, we consider two criteria for the rejected deposition holes [28]: (1) The full perimeter criterion (FPC), which considers that if a fracture intersects the wall of the disposal tunnel and penetrates the tunnel perimeter completely and the hypothetical linear extension of the fracture intersects the canister, that deposition hole is to be rejected (see the left deposition hole in red color shown in Figure 7). The canister in orange color is not intersected by the hypothetical linear extension of the fracture, so it will not be rejected. (2) The extended full perimeter criterion (EFPC), which considers that five deposition positions are intersected by connective fractures, states that these five or more deposition holes are to be rejected (see five deposition holes in a row on the right of Figure 7).

2.3.4. Calculation of Fracture Shear Movement

This study evaluates the shear displacement induced by earthquakes using the 3DEC model [24]. The fracture and earthquake simulation models calculate the shear displacement of fractures induced by earthquakes. For the earthquake simulations in 3DEC, the source data of earthquakes in Taiwan is collected and classified into fault source and diffuse seismicity. According to previous geological evolution investigations, the tectonic settings in the reference case will remain invariant within 1,000,000 years from now [16]. Three major parameters for seismic sources are geometry, maximum magnitude, and seismicity rate. These parameters are derived from current geological and seismic data, and their uncertainties are considered through a logic tree shown in Figure 8. The geometry, maximum magnitude, and seismicity rate for different types of seismic sources in the reference case are:

• Fault source: The active fault near the reference case is the Binhai fault (Figure 9a). Based on data from previous studies and workshops, the Binhai fault can be divided into two rupture models (Figure 8).
  • Model 1: The length and dip of the fault are 450 km and 60 degrees toward the west, respectively. Three possible seismogenic depths are 10, 15, and 20 km [16]. The maximum magnitudes evaluated from in situ stress are Mw 7.93, 8.27, and 8.51. Slip rates are 0.02, 0.2, and 0.5 mm/yr;
  • Model 2: The length of the fault is 71 km, and the dip of the fault plane is 54 degrees toward the east. The seismogenic depth in this region is around 25 km. The maximum magnitude of the Binhai fault can reach Mw 7.3 [29,30]. Slip rates are 0.02, 0.2, and 0.5 mm/yr.
• Diffuse seismicity: Diffuse seismicity is also called area sources in probabilistic seismic hazard analysis (PSHA). There is no definite geometry for the rupture plane of diffuse seismicity. A boundary with earthquakes shares similar focal mechanisms for diffuse seismicity could be defined through geophysical and geological surveys, narrowing down the uncertainties of rupture planes. Three boundaries of diffuse seismicity referring to previous studies for the reference case are listed as follows (see Figure 9b):
  • A circle with a 200 km radius from the reference case: The length of the radius is determined by the distance from the reference case to the deformation front in the Taiwan region. The edge of the circle also coincided with the front of Peikang High and the seismicity distribution in Taiwan [16,31,32];
  • AS_K01 and DS_K01: AK_K01 and DS_K01 are the area sources that cover the Taiwan Strait region. Their boundaries are modified from areal sources proposed in previous seismic hazard analyses for nuclear power plants in Taiwan and the design earthquakes in the Taiwan Strait region;
  • The upper and lower depths of diffuse seismicity are 2 km and 35 km, referring to the depths of areal sources in the past PSHA studies in Taiwan. Since this diffuse seismicity cannot be linked to any precise geological structure, a non-surface rupture model is assumed in the study. Based on the assumption and the study of Shimazaki (1986) [33], the maximum magnitude for the diffuse seismicity was specified as 6.5. This maximum magnitude of 6.5 is consistent with the observed seismic data within this region after eliminating events related to the Binhai fault [34]. Additionally, the seismicity rate of diffuse seismicity is derived based on the truncated exponential model [35]. The truncated exponential model employs Gutenberg–Richter’s law, which represents the counts of cumulative annual numbers of earthquakes of a certain magnitude [36].

The fractures with diameters ranging from 240 to 360 m were used to simplify the simulations for the fault source model. However, the diffuse seismicity model ranges from 280 to 320 m in the study. Because of the limited rock mechanic parameters of the fractures, we used the value of fracture cohesion, friction angle, dilation angle, residual cohesion, and residual friction angle proposed in this study of SKB [37].

2.3.5. Evaluation of Potential Canister Failure Due to Shear Movement of Fractures

The generated DFN and the conceptual repository model allow for analyzing the intersections between the two models. This study assumes that the fracture size remains constant throughout the safety assessment timescale. The number of intersections between fractures and canisters is determined by applying the deposition hole rejection criteria. For the remaining intersections that do not fulfill the criteria between fractures and canisters, the cumulative fracture shear movement induced by earthquakes is then calculated. In addition, we assume that each shear direction induced by an earthquake is the same. If the cumulative fracture shear movement exceeds 5 cm, we consider the canister failed.

Section 2.2.3 has presented the concept of fracture shear displacement (${\mathrm{u}}_{\mathrm{sim}}$) induced by an earthquake. The fracture size is denoted as ${\mathrm{r}}_{\mathrm{sim}}$. Based on the analytical solution of fracture in a linearly elastic medium [38], we could calculate the shear displacement $\left({\mathrm{u}}_{\mathrm{r}}\right)$ at different sizes ($\mathrm{r}$) or the shear displacement (${\mathrm{u}}_{{\mathrm{r}}^{\prime }}$) at various distances (${\mathrm{r}}^{\prime }$) from the fracture center:

$$\frac{{\mathrm{u}}_{\mathrm{r}}}{{\mathrm{u}}_{\mathrm{sim}}}=\frac{\mathrm{r}}{{\mathrm{r}}_{\mathrm{sim}}}$$

(10)

$$\frac{{\mathrm{u}}_{{\mathrm{r}}^{\prime }}}{{\mathrm{u}}_{\mathrm{r}}}=\sqrt{1-{(\frac{{\mathrm{r}}^{\prime }}{\mathrm{r}})}^2}$$

(11)

In most cases, the actual number of intersections between fractures and canisters can only be determined during the construction of an underground facility. Additionally, it is difficult to investigate the size of fractures accurately. Therefore, the study calculates the canister shear failure based on the stochastic method and the associated DFN realizations. The canister shear failure ratio $({\mathrm{p}}^{\prime }$) is defined as:

$${\mathrm{p}}^{\prime }=\frac{{\mathrm{N}}_{\mathrm{failure}}}{{\mathrm{N}}_{\mathrm{r}}{\mathrm{N}}_{\mathrm{hole}}}$$

(12)

where the notation ${\mathrm{N}}_{\mathrm{failure}}$ is the total number of canister failures in ${\mathrm{N}}_{\mathrm{r}}$ realizations, ${\mathrm{N}}_{\mathrm{r}}$ is the number of realizations (i.e., 500 realizations in this study) and ${\mathrm{N}}_{\mathrm{hole}}$ is the number of deposition holes in the simulation model. Note that 2860 DH were considered in the study.

2.4. Groundwater Flow and Particle Tracking

This study employs DarcyTools numerical model to simulate groundwater flow and advective transport. DarcyTools is a non-commercial code developed by SKB and is primarily used for groundwater flow and advective transport modeling in sparsely fractured rock at a deep geological disposal facility [39,40,41]. DarcyTools has previously been employed to calculate groundwater flow and advective transport for two potential sites in Sweden, namely Forsmark and Laxemar. As DarcyTools is based on the finite volume method (FVM), the properties of the generated fracture system must be transformed into effective properties in each computational cell, which belongs to the equivalent continuum porous medium (ECPM) model [42,43,44,45]. Previous studies have focused on comparing flow and transport results with different approaches for disposal facility analysis, such as FracMan (discrete approach) and hybrid approach. The acceptable upscaling concept has been validated based on quantitative comparisons [22,23,46,47]. The more detailed applications of the DarcyTools model can refer to various reports and previous studies [42,43,44,45]. Here we briefly describe the mathematical and numerical algorithms for the numerical model.

2.4.1. Mathematical Formulation for Flow

DarcyTools utilizes the mass conservation equation, expressed as follows for a two-component system (water and solid):

$$\frac{\partial \mathsf{\rho }\mathsf{\phi }}{\partial \mathrm{t}}+\frac{\partial }{\partial \mathrm{x}}\left(\mathsf{\rho }\mathrm{u}\right)+\frac{\partial }{\partial \mathrm{y}}\left(\mathsf{\rho }\mathrm{v}\right)+\frac{\partial }{\partial \mathrm{z}}\left(\mathsf{\rho }\mathrm{w}\right)=\mathrm{Q}$$

(13)

where $\mathsf{\rho }$ represents density, $\mathsf{\phi }$ denotes porosity, $\mathrm{t}$ indicates time, and $\mathrm{u}$, $\mathrm{v},$ and $\mathrm{w}$ represent the directional components of the Darcy velocity. Additionally, $\mathrm{x}$, $\mathrm{y},$ and $\mathrm{z}$ represent the location in a Cartesian coordinate system. $\mathrm{Q}$ represents the source/sink term per unit volume of fluid mass. Based on Darcy’s law, the Darcy velocity in both the mass conservation equation and the mass transport equation can be expressed as follows:

$$\mathsf{\rho }\mathrm{u}=-\frac{{\mathrm{K}}_{\mathrm{x}}}{\mathrm{g}}\frac{\partial \mathrm{P}}{\partial \mathrm{x}}$$

(14)

$$\mathsf{\rho }\mathrm{v}=-\frac{{\mathrm{K}}_{\mathrm{y}}}{\mathrm{g}}\frac{\partial \mathrm{P}}{\partial \mathrm{y}}$$

(15)

$$\mathsf{\rho }\mathrm{w}=-\frac{{\mathrm{K}}_{\mathrm{z}}}{\mathrm{g}}\frac{\partial \mathrm{P}}{\partial \mathrm{z}}-{\mathrm{K}}_{\mathrm{z}}\left(\mathsf{\rho }-{\mathsf{\rho }}_0\right)$$

(16)

where ${\mathrm{K}}_{\mathrm{x}}$, ${\mathrm{K}}_{\mathrm{y}},$ and ${\mathrm{K}}_{\mathrm{z}}$ is the hydraulic conductivity and ${\mathsf{\rho }}_0$ is the reference fluid density. Note that the dynamic pressure relative to the reference hydrostatic pressure is described as:

$$\mathrm{P}=\mathrm{p}+{\mathrm{P}}_0$$

(17)

where $\mathrm{p}$ is the gauge pressure. The fluid properties, such as dynamic viscosity and density, at an isothermal condition, are given by the state laws:

$$\mathsf{\mu }={\mathsf{\mu }}_0$$

(18)

$$\mathsf{\rho }={\mathsf{\rho }}_0\left[1+\mathsf{\alpha }\mathrm{S}\right]$$

(19)

where ${\mathsf{\mu }}_0$ is the reference viscosity, $\mathsf{\alpha }$ is a constant value, and $\mathrm{S}$ is the salinity.

2.4.2. Particle Tracking

DarcyTools employs the traditional particle tracking algorithm to simulate the transport dynamics in the numerical domain [42,43,44,45]. Here, we only consider the advective process and ignore sorption, degradation, and decay during transport. Particle movement is simulated by monitoring the particle location at each time step. We use the Euler method to calculate the particle movement and its location at each time step once the flow velocity is available in the simulation domain [48]:

$$\mathsf{{\rm X}}\left(\mathrm{t}+\Delta \mathrm{t}\right)=\mathsf{{\rm X}}\left(\mathrm{t}\right)+\mathrm{V}\left(\mathrm{x},\mathrm{t}\right)\Delta \mathrm{t}$$

(20)

where X represents the particle location in the space, t indicates the time, $\Delta \mathrm{t}$ is the specified time step, and $\mathrm{V}$ is the seepage velocity at a specific location. Particle tracking is conducted to calculate the potential radionuclide releases from the repository through different release pathways.

2.4.3. The Concept for Radionuclides Released from Canisters and Flow and Transport Properties

Isolation, containment, and retardation are the main safety functions for geological disposal in the KBS-3 disposal concept [3,5]. For each engineering barrier (i.e., canister, buffer, backfill) and natural barrier (i.e., host rock), secondary safety functions detail how containment and retardation are maintained. One of the secondary safety functions of the natural barrier is to provide a suitable hydrogeological environment, and two critical issues are to control transport resistance at the buffer/rock interface and to have high resistance in the geosphere to limit the transport of solutes. Expressing these functions as measurable or calculable quantities is desirable to evaluate the safety quantitatively. Two indicators related to the critical issues are equivalent flow rate (Q_eq) (m³/s) and flow-related transport resistance (F) (year/m) [3,5,10,11,12,13].

Radionuclides might be released from the canisters and migrate into the buffer and host rock once the canisters are damaged for the KBS-3 disposal concept. Among various release patterns, the migration of radionuclides into the mobile water in fractures surrounding the deposition hole (referred to as the Q1 path) is particularly significant [10,11,12,13]. Specifically, a shear fracture, which intersects a DH and destroys the canister inside the DH by shear movement, could be the major pathway for radionuclide transport [3,5]. Figure 10 shows the schematic view of a conceptual repository and the Q1 pathway for radionuclides to leave from the DH. The Q_eq represents the concentration of the compartments in contact with flowing water in a fracture. This is a fictitious water flow rate with a concentration equal to that at the interface between the buffer and host rock [49,50]. Q_eq depends on the geometry of the contact area, water flux, and diffusivity. In the study, the Q_eq for a fracture intersecting a DH is given by:

$${\mathrm{Q}}_{\mathrm{eq}}=2\mathrm{UW}\sqrt{{4\mathrm{D}}_{\mathrm{w}}{\mathrm{t}}_{\mathrm{c}}/\mathsf{\pi }}$$

(21)

where $\mathrm{U}$ represents an equivalent initial flux of a fracture averaged over the rock volume adjacent to the DH. $\mathrm{W}$ indicates the diameter of DH, D_w is the diffusion coefficient in free water, and t_c is the contact travel time, which is the time a water parcel is in contact with the compartment.

For examining the Q1 path released from the DH to the biosphere, DarcyTools utilizes the particle tracking method to simulate the transport dynamic [42,43,44]. In this study, a particle is released for the DH, which is intersected by a shear fracture, and the path length (L), advective travel time (t_w), and F are calculated accordingly. Note that the particle is only released at the fracture intersecting the DH with the highest flow rate. L is obtained by adding the path length of each time step, which could be described as:

$$\mathrm{L}=\sum \Delta \mathrm{l}$$

(22)

where Δl is the path length of a time step. t_w is obtained by monitoring the time needed for a particle transport from the start location to the end location, which is written as:

$${\mathrm{t}}_{\mathrm{w}}=\sum \frac{\mathsf{\varphi }\Delta \mathrm{l}}{\mathrm{q}}$$

(23)

where $\mathsf{\varphi }$ is kinematic porosity, and q is Darcy flux. Notation F is a measurement of the potential for retention and retardation of the radionuclide transport along the trajectory, which is given by:

$$\mathrm{F}={\displaystyle \sum }_{\mathrm{i}}{\left(\frac{{\mathrm{a}}_{\mathrm{r}}\mathrm{L}}{\mathrm{q}}\right)}_{\mathrm{i}}$$

(24)

where ${\mathrm{a}}_{\mathrm{r}}$ is flow-wetted surface per volume of rock, $\mathrm{L}$ represents the length of the cell, $\mathrm{q}$ indicates Darcy flux through the cell, and i is the fragment numbers of a flow path. The total F relies on the summation of piecewise F along the flow path of a particle.

In this study, U and Q_eq can be estimated based on the groundwater flow solution, and the L, t_w, and F could be calculated using the particle tracking method.

3. Results and Discussions

3.1. Earthquake Induces Potential Failure Canister Destroyed by Fracture Shear Movement

3.1.1. Calculation of Fracture Shear Movements for the Reference Case

Figure 11 illustrates the simulation results of fracture displacements induced by a single earthquake. In the figure, the x-axis represents the time in seconds, and the y-axis represents the shear displacements of fractures in millimeters. The coseismic shear displacement of each logic tree branch is evaluated by multiplying the fracture displacement induced by a single earthquake with the associated weightings. This result is further multiplied by the number of earthquakes to calculate the cumulative coseismic shear displacement. The results from each branch in the logic tree are summed together to determine the total cumulative coseismic shear displacement during the safety assessment period. The calculations are based on several assumptions, including (1) fractures do not propagate with the accumulation of shear displacements; (2) fractures accumulate coseismic shear displacement only, meaning that interseismic shear displacement is not considered in the study; (3) the strength of fractures remains consistent; (4) directions of shear displacement triggered by every source model are consistent; (5) the shear displacement of fractures is regarded as the permanent displacement; and (6) only earthquakes that coincide with the source model in the logic tree will occur.

Table 4. The maximum and average accumulation shear displacements for five different fracture clusters in the entire assessment period.

Fracture Cluster	Accumulation of Shear Displacement	Fault Source	Diffuse Seismicity	Total
1	Max. (mm)	65.88	122.16	188.04
	Avg. (mm)	23.68	71.26	94.94
2	Max. (mm)	45.40	178.15	223.55
	Avg. (mm)	9.00	96.71	105.71
3	Max. (mm)	36.11	47.85	83.96
	Avg. (mm)	12.00	5.09	17.09
4	Max. (mm)	51.87	122.16	174.03
	Avg. (mm)	11.07	90.09	101.16
5	Max. (mm)	13.44	45.30	58.74
	Avg. (mm)	2.54	22.40	24.94

lists the results of the maximum and average accumulation of earthquake-induced shear displacements of the fractures. These results demonstrate that the accumulation of shear displacements of fractures for the entire assessment period can exceed up to 5 cm of the design requirement of canisters, indicating the possibility of causing failure of intersected canisters. Therefore, further assumptions and calculations should be made to evaluate the probabilities of the canister failure.

3.1.2. The Information of Shear Fracture-Destroying Potential Canisters

The canister shear failure ratio is evaluated by applying 500 DFN realizations over the safety assessment timescale (see Figure 12). The earliest possible time for canister failure due to shear displacement (when the accumulated shear displacement exceeds 5 cm) is estimated to be around 172 thousand years after closure, with a canister shear failure ratio of approximately 1 × 10⁻⁶. As time progresses, fractures accumulate shear displacement, resulting in an increased canister shear failure ratio of 3 × 10⁻⁴ by the end of the assessment period. The changing rate of canister shear failure over time is attributed to several factors. Firstly, the magnitude of fracture shear movement induced by earthquakes depends on the sizes of the fractures in the model. Thus, larger fractures result in displacement exceeding 5 cm earlier. Secondly, the distribution of shear displacements for a fracture depends on the distance from the intersecting position to the fracture center. Therefore, the shorter the distance, the larger the displacement can be. Thirdly, the number of fractures with varying radii behaves differently. A larger quantity is often observed for a smaller radius. Nonetheless, the displacement of fractures tends to accrue with time, leading to a gradual reduction in the radius of fractures that could cause damage to the canisters.

Figure 13 shows the histogram depicting the distribution of canister shear failures across different deposition holes. The horizontal axis represents the DH numbers ranging from 1 to 30 based on the layout of DT, and the vertical axis shows the cumulative numbers of canister failures based on the 500 realizations. The analysis reveals that DHs 1 and 25 have a relatively high cumulative canister shear failure. In contrast, a lower number of failures is observed at DHs 26 to 30 compared to DH number 25, which can be attributed to the variation in tunnel length designs, resulting in no DHs numbered 26 to 30 assigned in the 250 m DT. This finding indicates that the DH at the end of the tunnel is more likely to escape the rejection criteria and cause more canister failures. Thus, DH numbers 1 and 25 in the 250 m DT and numbers 1 and 30 in the 300 m DT have a higher probability of canister failure induced by fracture displacement. Therefore, in the test cases, we choose three representative fractures leading to the canister failure by fracture displacement.

Table 5. Parameters of the three selected fractures leading to canister failures.

Fracture Case	Dip Direction	Dip Angle	Fracture Radius (m)	x	y	z	Failure Time (Year)
1	249.6	38.6	249.4	−1161.5	593.8	−599.7	234,000
2	9.6			101.1	111.8	−599.7
3				600.7	978.0	−599.7

lists the parameters of the three shear fractures. The relative locations of the fractures are shown in Figure 14. The assessment shows that the locations of shear fractures are close to the end of DT based on the calculation in Figure 13. The geometry, attitude, and location are input data for flow and advective transport modeling.

3.2. Density-Dependent Groundwater Flow, Advective Transport, and the Relevant Quantitative Parameters

3.2.1. The Groundwater Flow and the Advective Transport for Potential Failure Canisters

Figure 15 presents a DFN realization using the DFN recipes (in Table 3), three shear fractures (in Figure 14), and the layout of the conceptual repository. Figure 15a shows the top view of the DFN intersecting the conceptual repository, while Figure 15b presents a 3D view of the DFN and the conceptual repository. Three fracture cases have been merged into the DFN and are marked in blue (Case 1), green (Case 2), and red (Case 3). The calculations of the intersected numbers of DHs rely on the proposed DH rejection criteria. In the particular DFN realization (comprising three shear fractures), the total number of DHs intersected by fracture is 298. There are 148 DHs that will be rejected after applying the rejection criteria. Therefore, the remaining number of DHs is 150 (i.e., 298 minus 148). The DHs are the targets for Q1 path simulations and performance measure calculations. Details regarding the generation of computational grids, fracture analysis, and up-scaling of effective hydraulic parameters can be found in various reports and previous studies [22,23,39,40,41,42,43,44,45]. In this study, our primary focus is to demonstrate the workflow to quantify flow, advective transport, and performance measures of potential failure canisters destroyed by shear displacement.

Figure 16 shows the results of pressure and salinity distributions and the associated Q1 paths for the reference case. Figure 16a presents the top view of the pressure distribution and the 150 Q1 paths, while Figure 16b shows the pressure distribution and the 150 Q1 paths on the repository level and in two vertical planes across the conceptual repository for presentation purposes. Figure 16c presents the top view of the salinity distribution and the 150 Q1 paths, while Figure 16d illustrates the salinity distribution and the 150 Q1 paths on the repository level and in two vertical planes across the conceptual repository. The pressure results indicate higher pressure in the central area and gradually decrease from the center toward the coastline. This result suggests a radial flow of groundwater from the central area to the coastline. The salinity results show a distinct interface between seawater and freshwater beneath the island, which would affect flow paths accordingly. The number of released particles depends on the repository layout, fracture system, and geometrical rejection criteria. Based on these parameters, we have detected 150 potential release locations for the Q1 paths. The results reveal that the groundwater pressure field and salinity field significantly influence Q1 paths. Moreover, the particles will likely move towards the north, northwest, and northeast due to the hydraulic gradient controlled by the relatively high surface elevation in the central area. The presence of seawater and freshwater interfaces also influences flow paths near the coastline, where downward flow paths turn upward near the interface, resulting in a majority of the particles through the coastline.

Figure 17 shows the Q1 paths of the five potential failure canisters destroyed by the shear displacement, while

Table 6. The case numbers, DH numbers, and the start locations of the Q1 paths.

Case Number	DH Number	Center of DH	Start Location of the Particle
Case 1	DH-2524	(190,892.64, 2,710,417.50, −504.08)	(190,893.40, 2,710,418.00, −505.25)
	DH-2525	(190,884.84, 2,710,422.00, −504.08)	(190,885.10, 2,710,423.00, −502.75)
Case 2	DH-671	(191,922.34, 2,709,736.00, −504.08)	(191,922.90, 2,709,737.00, −506.75)
Case 3	DH-2833	(192,414.22, 2,710,606.75, −504.08)	(192,415.10, 2,710,607.00, −508.25)
	DH-2843	(192,422.02, 2,710,602.25, −504.08)	(192,422.90, 2,710,603.00, −505.75)

lists the fracture case number, DH numbers, the center of the DHs, and the start locations of the particles. The start locations of the particles heavily rely on the calculations of the intersections between the fractures and DHs. The highest flux value on an intersection is calculated to determine the start location of a Q1 path. Ultimately, three fracture cases could lead to the failures of five canisters by the shear displacement. Figure 17a shows the top view of the five Q1 paths from the conceptual repository to their end locations in the entire domain, while Figure 17b is the 3D view of the results. Figure 17c is the close view of three fracture cases, the conceptual repository, and the five Q1 flow traces. The patterns reveal that the Q1 paths released from DH-2524 and DH-2525 tend to move toward the northwestern direction, and the end location appears to be on the surface of the domain. On the other hand, the rest of the Q1 paths tend to migrate northward, and their end locations are on the coastline. Figure 17d is the close view of Q1 paths of DH-2524 and DH 2525 for fracture case 1, while Figure 17e is the translucent view of Figure 17d for presentation purposes. Figure 17f is the close view of the Q1 path of DH-671 and fracture case 2, while Figure 17g is the translucent view of Figure 17f. Figure 17h is the close view of the Q1 paths of DH-2833 and DH-2843 and fracture case 3, while Figure 17i is the translucent view of Figure 17h. The start locations of the Q1 paths are exactly at the intersection between the fractures and DHs. Thus, the z-value locations of the five DHs differ from one another.

3.2.2. Flow and Transport Parameters for Potential Failure Canisters Induced by the Shear Displacement

Table 7. Flow and transport parameters of the three fractures induced by the shear displacement and leading to the canister failures.

DH Number	U (m/Year)	Qeq (m/Year)	L (m)	tw (Year)	F (Year/m)	End Location of Particle	End Location Type
DH-2524	2.06 × 10−5	7.78 × 10−5	2.44 × 103	4.34 × 102	1.19 × 107	(189,213.70, 2,711,510.00, −1.00)	Shallow
DH-2525	2.14 × 10−5	7.92 × 10−5	2.43 × 103	4.39 × 102	1.18 × 107	(189,214.90, 2,711,050.00, −1.00)	Shallow
DH-671	1.28 × 10−5	6.13 × 10−5	2.31 × 103	6.25 × 102	2.50 × 107	(192,447.20, 2,711,510.00, −47.95)	Lateral boundary
DH-2833	1.66 × 10−5	6.98 × 10−5	1.21 × 103	1.26 × 102	6.29 × 106	(192,884.30, 2,711,382.00, −25.68)	Lateral boundary
DH-2843	1.74 × 10−5	7.15 × 10−5	1.21 × 103	1.24 × 102	5.93 × 106	(192,880.50, 2,711,382.00, −25.48)	Lateral boundary

presents the flow and transport parameters for five potential failure canisters induced by the shear displacement. The performance measure calculations employ the U and Q_eq in the near field and use the L, t_w, and F in the far field. The results show that DH-671 is destroyed by fracture case 2 and has the lowest U and Q_eq. The DH-2524 and DH-2525 are destroyed by fracture case 1 and have the highest U and Q_eq. Two near-field parameters are essential input data for buffer erosion and canister corrosion calculations. DH-2524 and DH-2525 have high L values because their trajectories move toward the northwest and yield longer migration distances. However, DH-2833 and DH-2843 have low L and t_w and high F values. The study also presents the results of the end locations of particles and the end location types. Since the reference case is an offshore island, the end locations and the associated types are important for biosphere calculations. The results show that the end locations of the Q1 paths for DH-2524 and DH-2525 are very close to the domain surface. Here, we assume that this type of end location is called shallow, representing that the radionuclides might be on the shallow formation or aquifer. Other end location types of the Q1 paths are on the lateral boundary, showing that the radionuclides could migrate to the sea area. The exact end locations of the Q1 paths could be used in biosphere calculations, such as the biosphere dose conversion factor (BDCF) [16].

4. Conclusions

The study has developed a workflow for simulating and calculating performance measures of potential canister failures induced by earthquake-induced shear movements. Specifically, a conceptual repository and geological conceptual model of a reference case were used as the hydrogeological conceptual model for flow and advective transport simulations. The study used 3DEC to model fracture shear movements induced by earthquakes based on the Coulomb slip model. Two main numerical models, including the fracture generation and the earthquake simulation models, have been integrated into the workflow. For earthquake simulations in 3DEC, the source data of earthquakes in Taiwan were collected and classified into fault sources and diffuse seismicities. Three major parameters, geometry, maximum magnitude, and seismicity rate of seismic sources, were derived from current geological and seismic data. The uncertainties of these parameters were considered through a logic tree. Results showed that the accumulation of shear displacement of fractures for the entire assessment period could exceed the 5 cm design criterion of the canisters, indicating the possibility of causing the failures of intersected canisters. This information is the fundamental input data for calculating the probabilities of canister failures.

This study has calculated the intersection pattern between the fracture system and the conceptual repository based on information derived from the DFN recipes and the layout of the conceptual repository. Two deposition hole rejection criteria (i.e., FPC and EFPC) have been proposed to evaluate the canister failures in the conceptual repository. The study assumed that each shear direction induced by an earthquake was the same. A canister was considered to have failed if the cumulative fracture shear movement exceeded 5 cm. In addition, the study used the 3DEC model to calculate the canister shear failure ratio and the cumulative numbers of the failure canisters. The results showed three possible canister failure cases at the end of the tunnel based on 500 DFN realizations over the safety assessment timescale. The geometries, attitudes, and locations of the three fracture cases are the input data for the flow and advective transport simulation.

This study used the DarcyTools model to simulate the flow and advective transport for a radionuclide waste repository in an offshore reference case in Taiwan. Specifically, three potential canister failure cases were considered to accurately calculate the performance measures of failed canisters induced by fracture shear movement. The study considered the Q1 paths released from the canisters based on the KBS-3 disposal concept. The results of pressure in the central area of the island were relatively high, with a gradual decrease towards the coastline. The salinity results showed an evident interface between seawater and freshwater beneath the island, which could influence the flow paths. The particle tracking algorithm was employed for the Q1 path simulations. Results demonstrated a strong influence of the groundwater pressure and salinity fields on the flow paths. The particles were found to move predominantly towards the north, northeast, and northwest, following the gradient caused by the central area. At the same time, the seawater and freshwater interface control the flow paths near the coastline. The performance measures calculations for five potential failed canisters (DH-671, DH-2524, DH-2525, DH-2833, and DH-2843) have successfully determined the critical parameters of buffer erosion, canister corrosion in the near field, and radioactive migration in the far field, namely U, Q_eq, L, t_w, and F. Additionally, the end locations and types of different paths were detected, providing the crucial information for biosphere calculations.