1 Introduction

Excavation of a deep underground opening induces fractures around the opening due to stress redistribution. The development of such fractures may lead to spalling of the rock mass and instability of the opening. In particular, in a high-level radioactive waste (HLW) disposal project, the zone, where fracture development increases permeability, can provide a pathway for the migration of radionuclides around the facility (Japan Nuclear Cycle Development Institute 2000; Tsang et al. 2005). Such a zone is called excavation damaged zone (EDZ) (Tsang et al. 2005). To accurately understand the extent of the EDZ, it is important to create a method for predicting the hydro-mechanical characteristics in and around the EDZ. In particular, estimation of the spatial distribution of hydraulic conductivity in and around the EDZ is necessary for accurate assessment of the migration of radionuclides in an HLW disposal project.

Many studies have attempted to clarify the extent and hydro-mechanical characteristics of the EDZ around an excavation in underground research laboratories (URLs), based on in situ surveys such as borehole televiewer (BTV) surveys, rock core observations, fracture mapping, deformation measurements, and seismic surveys (e.g., Sato et al. 2000; Bossart et al. 2002, 2004; Shao et al. 2008; Baechler et al. 2011; Aoyagi et al. 2014, 2017; Armand et al. 2014). In addition to in situ surveys, some studies have performed hydro-mechanical coupling analysis to estimate the hydro-mechanical behaviour, such as the change in pore pressure (Guayacan-Carrilo et al. 2017; Souley et al. 2017) and hydraulic conductivity (Inoue et al. 2005; Rutqvist et al. 2009; Shen et al. 2011; Chen et al. 2014) around the EDZ.

Ishii (2015) proposed a method to estimate the highest potential transmissivity of fractures in fault zones using an index calculated as the ratio of the effective mean stress to the tensile strength of intact rock, known as the Ductility Index (DI) or the Mean Stress Index (MSI) (Fig. 1). In this paper, we use the term MSI. A correlation between MSI values and the highest potential transmissivity of fractures in fault zones has been reported, and an empirically derived power law relationship between transmissivity and MSI has been proposed (i.e., the MSI model) (Fig. 2). In a fault zone, permeable fractures can form as a result of shear-induced dilation due to stress concentrations generated at the asperities/tips of faults (Fig. 3a–c; Ishii 2015, 2016).

Fig. 1
figure 1

Concept of the Mean Stress Index (MSI), after Ishii (2015). \({\sigma _{\text{t}}}\) is the tensile strength, \({\sigma _{\text{c}}}\) is the uniaxial compressive strength, and B is equal to \({{{\sigma _{\text{t}}}} \mathord{\left/ {\vphantom {{{\sigma _{\text{t}}}} {{\sigma _{\text{c}}}}}} \right. \kern-0pt} {{\sigma _{\text{c}}}}}\)

Full size image
Fig. 2
figure 2

a Transmissivities of flow anomalies detected in fault zones in six regions worldwide (including discrete shear fracture systems in the diatomaceous mudstone of the Koetoi Formation) plotted against the corresponding value of the Mean Stress Index (MSI) (modified from Ishii 2017b). Horizontal bars show errors due to uncertainties in tensile strength. The MSI model line (solid line) and standard error (dashed lines) of Ishii (2015) are also shown. b Transmissivities of flow anomalies detected in fault zones in mudstone formations (i.e., the Wakkanai and Palfris formations, and shear fracture systems in the Koetoi Formation) plotted against the corresponding MSI

Full size image
Fig. 3
figure 3

Reproduced with permission from Diederichs (2003)

a Schematic diagram of the opening of secondary fractures caused by shear-induced tensile stresses in a damage zone. b Opening by mismatch of the shear fracture walls (modified from Ishii 2015). c Sketch of a sinistral strike-slip fault zone exposed in a horizontal outcrop of the Wakkanai Formation at the Horonobe site, annotated with the directions of the maximum and minimum principal stresses (Ishii 2015). d Schematic diagram of shear-induced tensile fractures developed around the wall of a gallery.

Full size image

Similarly, in the EDZ, permeable fractures have formed around gallery walls in various URLs as a result of shear-induced dilation due to stress concentrations generated at the asperities/tips of pre-existing flaws such as microcracks, grain boundaries, and low-aspect-ratio pores (Fig. 3d; Diederichs 2003, 2007; Yong et al. 2010, 2013; Armand et al. 2014; Aoyagi et al. 2017). However, factors such as the excavation method, excavation geometry, rock support, distribution of natural (pre-existing) fractures, and stress field can affect the failure behaviour of the EDZ, and tensile fracturing can occur if the differential stress is high. Therefore, the formation mechanism differs between permeable fractures in a fault zone and in the EDZ. However, openings in these fractures can be maintained by continuous shear-induced tensile stresses (Fig. 3), and the highest potential hydraulic conductivities of such fractures can be controlled by the MSI (Ishii 2015, 2017b). Therefore, estimation of the highest potential hydraulic conductivity using the MSI model can also be applied to the EDZ, utilising the similar degree of preservation of the opening in the EDZ and in a fault zone.

This paper focuses on the development of a method for estimating the highest potential hydraulic conductivity of the EDZ induced around a gallery in the Horonobe URL, Japan. BTV surveys, rock core observations, displacement measurements (i.e., convergence), and hydraulic tests were performed to understand the mechanical and hydraulic characteristics, such as the hydraulic conductivity, deformation behaviour, and extent of the EDZ. In addition, we performed a hydro-mechanical coupling analysis. The validity of the analysis was verified by comparing the results with those of the in situ surveys. Using the effective stress distribution derived from the analysis and distribution of fractures detected by the BTV survey, the highest potential hydraulic conductivity was estimated from the MSI model, and this was then compared with the in situ hydraulic test data. We also discuss the applicability of the MSI model in estimating the highest potential hydraulic conductivity in the EDZ.

2 Study Site

2.1 Horonobe Underground Research Laboratory

The Japan Atomic Energy Agency (JAEA) has been undertaking research and development at Horonobe (Hokkaido, Japan) to enhance the reliability of technology for the deep geological disposal of HLW in sedimentary rocks at this off-site URL. The layout of the Horonobe URL is shown in Fig. 4. The ventilation shaft and east access shaft were excavated to 380 m below ground level, and the west access shaft was excavated to 365 m below ground level. Excavations of horizontal galleries were also completed at depths of 140, 250, and 350 m.

Fig. 4
figure 4

Layout of the Horonobe URL

Full size image

We undertook in situ experiments in the gallery at 350 m depth, which was excavated in Neogene siliceous mudstone of the Wakkanai Formation that contains opal–cristobalite/tridymite. Table 1 lists the average values and standard deviations of the mechanical properties of the siliceous mudstone, as obtained from laboratory tests following ISRM methods (International Society for Rock Mechanics 1978a, b, 1979a, b) on rock cores collected from a borehole drilled in the 350 m gallery. The uniaxial compressive strength, Young’s modulus, and Poisson’s ratio were obtained from ten uniaxial compressive tests using samples measuring 30 mm in diameter and 60 mm in length. The tensile strength was obtained from 10 Brazil tests using samples measuring 30 mm in diameter and 30 mm in length. The density, effective porosity, water content, and seismic velocity were measured on ten samples.

Table 1 Mechanical and physical properties of siliceous mudstone in the 350 m gallery
Full size table

For the in situ experiments, we selected the gallery in Niche No. 3 for analysis, which is located at 350 m depth (Fig. 5). The niche had been excavated mechanically, without blasting, from July to August 2013, using an excavator equipped with a hydraulic impact hammer. Figure 6 shows a cross section of the horseshoe-shaped niche, which is 4.0 m wide, 3.2 m high, and 25 m long. The niche is supported mainly by 200 mm-thick shotcrete, 2 m-long rockbolts, and steel arch ribs. The support structures were installed every 1 m within the first 5 m section from the entrance along the excavation of the niche, and then every 1.5 m for the last 5 m section to the end of the niche. Rockbolts were not installed in the final 10 m of the 25 m-long excavation to ensure they did not penetrate the boreholes used for the in situ experiments.

Fig. 5
figure 5

Layout of the 350 m gallery and the borehole for surface-based investigation (SAB-1)

Full size image
Fig. 6
figure 6

Vertical cross section of Niche No. 3

Full size image

2.2 In Situ Experiments

2.2.1 Fracture Observations

We examined the distribution of fractures on the sidewall and below the floor of Niche No. 3 using BTV surveys of the boreholes, as shown in Figs. 7 and 8. The BTV has a borehole image processing (BIP) system with a high-resolution charge coupled device (CCD) camera that can capture a full-colour image of the borehole wall and detect fractures with an aperture of ~ 0.1 mm.

Fig. 7
figure 7

Locations of boreholes for BTV surveys and convergence measurements. a Plan view. b Cross-sectional view

Full size image
Fig. 8
figure 8

Locations of boreholes, sections subjected to hydraulic tests, and convergence measurements. a Vertical section perpendicular to the tunnel axis. b Cross-sectional view of the niches, boreholes, and test sections

Full size image

We also studied rock cores to analyse the fracture characteristics and their distribution. The fractures were classified into four types: shear, tensile, hybrid, and drilling/coring/handling-induced fractures. This classification follows the previous studies (Ishii 2016, 2017a), based on the structural geology of the siliceous mudstone. Shear fractures are characterized by slickenlines or slickensteps, whereas tensile fractures are characterized by natural plumose structures with arrest lines; hybrid fractures are defined by the coexistence of both these characteristics (Ishii 2016). Drilling/coring/handling-induced fractures are characterized by well-developed hackle plume components, which originate within the core, and by twist and inclusion hackle faces and steps that are geometrically related to the core boundary (Kulander et al. 1990).

Aoyagi et al. (2017) reported that all of the fractures in the niche are tensile or hybrid fractures with hackle plumes and arrest lines on their surfaces. These fractures are developed pervasively over the niche wall, subject to the wall geometry. As such, it was concluded that these fractures are EDZ fractures that were induced by gallery excavation. Furthermore, there are no natural fractures according to the fracture map on the sidewall. This observation is consistent with the results of the surface-based borehole investigation (SAB-1; Fig. 5), which did not identify any fractures in the core section near Niche No. 3 (Suko et al. 2013). Given this, we consider that the EDZ is marked by the zone, where tensile and hybrid fractures are densely developed in the vicinity of the sidewall and floor. Although the extent of the EDZ is difficult to precisely determine, a densely fractured EDZ would provide a pathway for the migration of radionuclides around an HLW disposal site. Thus, in this study, we assumed that the extent of the densely fractured zone is the same as that of the EDZ.

Figure 9 shows an example image of the EDZ taken from the BTV survey and rock core observations in borehole H2-1. The BTV image shows that fractures are densely developed at 0–0.2 m from the niche wall (Fig. 9a), but not beyond this zone. In this fractured zone, hackle plumes and an arrest line are present on the surface of the rock core (Fig. 9c). Given that the location of the origin of the hackle plumes is outside the core surface, the fracture around the niche wall must be a tensile fracture that was not induced by drilling. In contrast, core from the section farther than 0.2 m from the niche wall displays different types of fractures. For example, Fig. 9d shows a hackle plume component, whose origin was on the surface of the rock core and a twist hackle step, which is also geometrically related to the core boundary. These features represent evidence of drilling/coring/handling-induced fractures. Based on these observations, we consider that the extent of the EDZ in this borehole is 0–0.2 m from the niche wall. The other boreholes shown in Figs. 7 and 8 exhibit similar fracture types, distribution trends, and characteristics as borehole H2-1 (Aoyagi et al. 2017).

Fig. 9
figure 9

Fracture distribution observed by BTV surveys and rock core analysis (after Aoyagi et al. 2017). a BTV image of borehole H2-1, which is 0–0.6 m from the niche wall. An enlarged image of the region 0–0.25 m from the wall is shown above (a). b Photograph of core from borehole H2-1. c and d Photographs of the fracture surface at 0.15 and 0.50 m from the niche wall (see text for details)

Full size image

2.2.2 Convergence Measurements

To measure the deformation behaviour of the niche during excavation, the horizontal displacement (i.e., convergence) on the wall of the shotcrete after installation of support material was measured in four sections in Niche No. 3 (Figs. 7a, 8a), using a convergence measurement instrument (NH-15F). The resolution accuracy of the convergence measurement is 0.1 mm. The measurement point was attached to the surface of the shotcrete and the measured line was located along the spring line (Fig. 8a). The measurement in each section was performed for each face progressively after installation of the measurement points. The measurement was continued until the measured value was the same as the previous value.

2.2.3 Hydraulic Tests

We measured the change in hydraulic conductivity that was induced by excavation around Niche No. 3. Figure 8 shows the layout of boreholes H2-1, H3-2, and H4-2, and measured test sections. Boreholes H2-1 and H4-2 were drilled at angles of 15° and 16° below the horizontal. Borehole H3-2 was drilled at 15° below the horizontal from the floor of Niche No. 3, and parallel to the long axis of the niche. The lengths of boreholes H2-1, H3-2, and H4-2 were 10.3, 12.5, and 13.5 m, respectively, and their diameters were all 98 mm. Five test sections in borehole H2-1, one test section in borehole H4-2, and four test sections in H3-2 were isolated by 30 cm-long packers. To avoid relaxation of the packers, we monitored the pressure and held it constant. Note that the pressure measured in the hydraulic tests is a gauge pressure.

We injected water into section 5 of borehole H2-1, sections 3 and 4 of borehole H3-2, and section 5 of borehole H4-2 at a constant flow rate. When the pressure in the test section was stable for more than 1 h, we considered that it had achieved a steady state. The pressure in the section was increased using water injection until it reached a steady state of 0.05–0.10 MPa.

In sections 1–4 of borehole H2-1, located > 0.6 m from the niche wall, and sections 1 and 2 of borehole H3-2, located > 1.8 m from the floor of the gallery, we injected water at a constant pressure, because it took a long time to achieve steady state, due to the low hydraulic conductivity when injecting at a constant flow rate. The injection pressure in these sections was held at 0.1–0.5 MPa to avoid breaking the borehole wall.

Hydraulic conductivity was calculated using the Hvorslev equation (Hvorslev 1951). This assumes a steady-state laminar flow rate in a homogeneous and isotropic media around the borehole, and is formulated as follows:

$$k=\frac{Q}{{2\pi L}}\frac{{{\rho _w}g}}{p}\ln \left( {\frac{{2L}}{D}} \right),$$
(1)

where k is hydraulic conductivity (m/s), Q is the injection volume of water per second (m3/s), L is the length of the section (m), D is the diameter of the borehole (m), ρw is the density of water (kg/m3), g is acceleration due to gravity (9.8 m/s2), and p is the differential pressure between injection pressure and pore pressure (Pa). Relaxation of the packers did not occur during the tests, as the pressure of the packers was monitored and did not fall below 0.7 MPa.

In addition to the hydraulic tests, the pore pressures in all sections of the hydraulic tests were measured to monitor the change in pore pressure in and around the EDZ with time.

3 Hydro-mechanical Characteristics of the EDZ

3.1 Extent of the EDZ

Figure 10 shows the relationship between the cumulative number of fractures counted on the BTV images of each borehole (Fig. 7) and distance from the niche wall. A significant increase in the cumulative number of fractures was detected within 0.1–0.6 m of the wall. The fractures detected in the ranges by the BTV surveys correspond with the tensile or hybrid fractures observed in the rock cores. Thus, the extent of the EDZ is 0.1–0.6 m from the sidewall of the niche (e.g., Figs. 9, 12a, c).

Fig. 10
figure 10

Fracture distribution detected by the BTV survey. Labels are borehole numbers, as indicated in Fig. 7

Full size image

The distribution of below the niche floor was also investigated in borehole H3-2 by the BTV survey (Fig. 12e). Fractures are densely developed up to 1.6 m below the gallery floor, and the extent of the EDZ is 1.6 m below the gallery floor.

3.2 Convergence Measurements

Figure 11 shows the convergence curve during the excavation of Niche No. 3. In this figure, measured diameter is plotted as distance from the excavation face, which is normalized to the diameter of the Niche (D = 4.0 m). The convergence increased due to the excavation; once the distance from the excavation face was 1.5 to 2.0 × D, the rate of convergence increase was levelled off (Fig. 11). The maximum convergence in each measured section was 1.6–2.7 mm.

Fig. 11
figure 11

Relationship between measured convergence and distance from the excavation face

Full size image

3.3 Hydrological Characteristics of the EDZ

3.3.1 Hydraulic Conductivity

Figure 12a, b shows the locations of fractures and hydraulic conductivities measured in borehole H2-1 in October 2013, 1–2 months after completion of the excavation of Niche No. 3. In sections 2–4 of borehole H2-1, hydraulic conductivity varies between 2.88 × 10−11 and 4.06 × 10−11 m/s for the sections located farther than 0.6 m from the niche wall, which corresponds to the no-fracture zone (Fig. 12a). The presence of fractures at 7.6, 8.1, and 9.0 m from the niche wall increased the hydraulic conductivity in section 1 of borehole H2-1 to 6.17 × 10−10 m/s.

Fig. 12
figure 12

Results of hydraulic tests. a Fracture distribution in borehole H2-1. b Hydraulic conductivity in borehole H2-1. c Fracture distribution in borehole H4-2. d Hydraulic conductivity in borehole H4-2. e Fracture distribution in borehole H3-2. f Hydraulic conductivity in borehole H3-2

Full size image

In contrast, hydraulic conductivity is 3.53 × 10−6 m/s in section 5 of borehole H2-1, at distances of 0.06–0.26 m from the niche wall. This section corresponds to the fracture-developing zone (i.e., EDZ) detected by the BTV survey (Fig. 12a). The hydraulic conductivity is four-to-five orders of magnitude higher than that in sections 1–4 in borehole H2-1.

Figure 12c, d shows the locations of fractures and hydraulic conductivities measured in borehole H4-2 in April 2014, 8 months after completion of the excavation of Niche No. 3. The hydraulic conductivity is 8.43 × 10−7 m/s at distances of 0.14–0.74 m from the niche wall. This section corresponds to the EDZ detected by the BTV survey (Fig. 12c).

The hydraulic conductivities and fracture locations below the gallery floor (borehole H3-2) are shown in Fig. 12e, f. The hydraulic conductivities were measured in October 2013, as for borehole H2-1. Hydraulic conductivities in sections 1 and 2 of borehole H3-2 were 7.29 × 10−10 and 1.54 × 10−10 m/s, respectively, which are one order of magnitude higher than those of sections 1–4 of borehole H2-1, due to the presence of fractures at depths of 2.2 and 2.5 m below the niche floor. In sections 3 and 4, which correspond to the EDZ detected by the BTV surveys, the hydraulic conductivities were 6.41 × 10−7 and 1.24 × 10−8 m/s, respectively, which are two-to-three orders of magnitude higher than those in sections 1 and 2.

3.3.2 Pore Pressures

Figure 13 shows the pore pressures measured in the hydraulic test sections after the tests had been conducted. The pore pressure just before the excavation of Niche No. 3 was ca. 2.4 MPa (Mezawa et al. 2017), as shown by the blue dashed lines in Fig. 13. In section 5 of borehole H2-1, which corresponds to the fracture-developing and high hydraulic conductivity section (i.e., in the EDZ), the pore pressure is ca. 0 MPa (Fig. 13a). The pore pressure in section 5 of borehole H4-2 is also ca. 0 MPa. The pore pressures in sections 3 and 4 in borehole H3-2 are ca. 0 MPa, which is the same trend as observed in section 5 of boreholes H2-1 and H4-2. The pore pressures gradually increase with distance from the niche wall and outside of the EDZ (Fig. 13b).

Fig. 13
figure 13

Pore pressures in the hydraulic test sections. a and b Results of the measurements performed in boreholes H2-1 and H3-2, respectively

Full size image

4 Numerical Analysis

We performed hydro-mechanical coupling analysis, using the commercial code FLAC3D (v. 4.00.89) (Itasca Consulting Group 2009), to simulate the hydro-mechanical behaviour of a rock mass during excavation. The validity of the analysis results, such as the extent of fracturing, measured convergence, and pore pressure distribution, was verified by comparison with the results of in situ surveys.

4.1 Modelling Approach

4.1.1 Initial Stress and Pore Pressure Conditions

Initial stress conditions were estimated on the basis of the hydraulic fracturing tests performed in the boreholes around the Horonobe URL drilled prior to the excavation (surface-based investigation) and in the 350 m gallery (Aoyagi et al. 2013). The magnitude of stress around the facility decreased due to excavation, as a result of a pore pressure decrease (Aoyagi et al. 2013). It is also reported that the magnitude of the stress within an oil reservoir is reduced as a result of the decrease in pore pressure associated with depletion. This can be explained on the basis of linear poroelastic theory (Zoback 2007). Thus, in our analysis, the appropriate initial stress conditions took into account this change in pore pressure.

The pore pressure at 350 m depth prior to the exaction was ca. 3.5 MPa, based on pore pressure measurement in the boreholes. However, the pore pressure was ca. 0 MPa when measuring the stress in the 350 m gallery, because there was no water inflow from the measured borehole. The relationship between stress and pore pressure was estimated, as shown in Fig. 14, which shows the minimum principal stress \(\left( {{\sigma _{\text{x}}}} \right)\) and vertical stress \(\left( {{\sigma _{\text{z}}}} \right)\), which correspond to the horizontal and vertical stresses acting on the niche, respectively (Fig. 15). As described above, the pore pressure just before the excavation was ca. 2.4 MPa, and as such, the magnitude of initial pore pressure was set to 2.4 MPa in our analysis. The initial stress value corresponding to a pore pressure of 2.4 MPa was determined on the basis of the relationship between pore pressure and stress (Fig. 14). The magnitudes of the stresses are shown in the left panel of Fig. 15.

Fig. 14
figure 14

Initial stress and pore pressure conditions applied to the model analysis

Full size image
Fig. 15
figure 15

Analysed model and boundary conditions. The initial stress condition is shown at left, and an enlarged view of the model around Niche No. 3 is shown at right

Full size image

4.1.2 Analysis Model

Figure 15 shows the geometry of the model used for the numerical analysis. Only half of the domain was modelled for symmetry reasons. The niche was 4.0 m wide and 3.2 m high. The installation of shotcrete (width of 0.2 m) and the steel arch rib were also simulated. The rock mass and shotcrete were modelled as solid elements, whereas the steel arch rib was modelled as a beam element. The size of each element was 0.1 m within 1.0 m of the niche wall, and 0.2 m for distances of 1.0–2.0 m from the niche wall. Roller boundary conditions were applied to the left and bottom sides of the model. The pore pressure was fixed at the right and top sides of the model.

4.1.3 Hydro-mechanical Properties

Table 2 lists the properties of the rock, shotcrete, and steel arch rib used in the analysis. The rock properties were determined on the basis of the data in Table 1. The poroelastic coefficients (i.e., Biot’s and Skempton’s coefficients) were determined on the basis of the original laboratory tests described by Miyazawa et al. (2011). The dilation angle is assumed to be 0° on the basis of the post-failure behaviour of rock specimens measured in triaxial compression tests using siliceous mudstone (Yamamoto et al. 2005) and previous studies of mudstones (Wileveau and Bernier 2008; Xu et al. 2013). The hydraulic conductivity of intact rock was determined from laboratory permeability tests (Kurikami et al. 2008). The mechanical properties of shotcrete were determined from the standard specifications of concrete materials (Japan Society of Civil Engineers 2002), and those of steel arch ribs from the Japanese Chronological Scientific Tables (National Astronomical Observatory 2005).

Table 2 Properties of the rock, shotcrete, and steel arch rib used in the model analysis
Full size table

4.1.4 Failure Criterion and Post-failure Behaviour

The failure criterion for the rock mass was determined from the stress conditions under which a crack initiates along a developing failure within a specimen, as observed in laboratory tests. Crack initiation (CI) stress is defined as that when the volumetric or lateral stress–strain curve deviates from linearity, and this stress condition corresponds to the initiation of a wing crack (Bieniawski 1967; Lajtai 1974; Perras and Diederichs 2014). In these previous studies, the magnitude of CI stress was ~ 30–50% of the peak strength.

When the tangential stress exceeds the value of the CI stress in hard crystalline rock, visible tensile fractures develop and result in macroscopic spalling failure (Diederichs 2003, 2007). In galleries excavated in soft sedimentary rock formations, such as the Opalinus clay in the Mont Terri rock laboratory in Switzerland and in the siliceous mudstone formation in the Horonobe URL, tensile fractures develop when the tangential stress exceeds the CI stress (Yong et al. 2010, 2013; Aoyagi et al. 2017). Thus, in our study, the failure criterion based on the CI stress was determined and applied in estimating the extent of the EDZ.

To determine the failure criterion for the CI stress, we carried out ten uniaxial compression tests and 22 consolidated–drained triaxial compression tests using rock cores obtained from boreholes around the Horonobe URL and its galleries (Ota et al. 2010; Aoyagi et al. 2015). The CI stresses were determined by detecting the point at which the volumetric or lateral stress–strain curves deviated from linearity. Mohr circles corresponding to the stress conditions at CI obtained in the uniaxial and triaxial compression tests are shown in Fig. 16. Following the ISRM suggested method (International Society for Rock Mechanics 1983), the values of apparent cohesion c′ and apparent friction angle \(\phi ^{\prime}\) were calculated to be 2.37 MPa and 17.5°, respectively, as indicated in Table 2. Using these values, the Mohr–Coulomb criterion corresponding to CI stress employed with tension cutoff at the tensile strength \({\sigma _t}\) (1.83 MPa; Table 1) was determined (Fig. 16).

Fig. 16
figure 16

Mohr stress circles corresponding to CIs obtained from uniaxial and triaxial compression tests. The straight line indicates the determined Mohr–Coulomb criterion corresponding to CI stress with tension cutoff

Full size image

When the Mohr circle of each element intersects the failure criterion, post-failure behaviour is simulated, according to the plastic potential function described by Vermeer and de Borst (1984) and Itasca Consulting Group Inc (2009) (i.e., the non-associated flow rule). In addition, the pore pressure is reduced to 0 MPa, on the basis of the result that the measured pore pressure in the EDZ was ca. 0 MPa (Fig. 13).

4.1.5 Analysis Procedure

The workflow of the analysis procedure is shown in Fig. 17. In the analysis, initial stress and pore pressure conditions were simulated prior to the excavation analysis (Step 1). Excavation analysis was accomplished by reducing the forces that were required to maintain equilibrium with the initial stress state (i.e., equivalent excavation force on the excavation surface). Prior to the installation of the support, 78% of the equivalent excavation force was reduced to simulate the situation that corresponds to the advance of excavation faces 1.5 m from the representative section (Step 2). Subsequently, the shotcrete and steel arch rib were installed, and then, 100% of the equivalent excavation force was reduced (Step 3).

Fig. 17
figure 17

Workflow of the analysis procedure

Full size image

In a previous study, extensional EDZ fractures were observed in the Horonobe, Bure, and Mont Terri URLs, which were directly caused by short-term excavation-induced unloading (i.e., undrained elasto-plastic behaviour; Blümling et al. 2007). These phenomena can be observed in low-permeability formations (Martin and Lanyon 2003; Rutqvist et al. 2009). Therefore, in Steps 2 and 3 (i.e., the period of excavation), the excavated wall was set to be in an undrained condition.

After Step 3, the EDZ and excavated wall were set to a drained condition (i.e., the pore pressures of the nodes corresponding to the EDZ and excavated wall were decreased to 0 MPa). The analysis was then continued for 85 days, which corresponds to the period between the completion of the excavation and completion of the hydraulic tests (Step 4). The analysis results, such as extent of the EDZ, convergence, and pore pressure distribution, were then compared with the in situ measurements to confirm the model validity. Of note, when fluid flows through a porous medium, the solid weight, buoyancy, and seepage force should be considered to explain the force equilibrium (Wang 2000). These forces are taken into account in our analysis.

After the analysis, the highest potential hydraulic conductivity in the EDZ was estimated using the MSI model (Step 5), as described in section 5, and then compared with the results of the hydraulic tests.

4.2 Analysis Results

Figure 18 indicates the failure region or the distribution of the EDZ around Niche No. 3 after excavation (i.e., Step 3). The maximum extent of the EDZ is 0.9 m in the sidewall, 0.3 m at the crown, and 1.3 m below the floor. This result is largely consistent with the maximum extent of the EDZ detected by the in situ tests (sidewall: 0.6 m, floor: 1.6 m; Sect. 3.1).

Fig. 18
figure 18

Extent of the EDZ around Niche No. 3 (Step 3)

Full size image

The extent of the EDZ shows a marked increase at Step 3 (i.e., end of the excavation). In Step 4, the extent of the EDZ is unchanged from Step 3. This indicates that the EDZ developed due to the stress concentration and undrained behaviour of the niche wall during excavation.

Figure 19 shows the distributions of the maximum and minimum effective principal stresses in Steps 3 and 4. In the figure, the compressive stress is taken to be positive. In Step 3, the magnitude of the maximum effective principal stress increased around the niche wall due to the stress concentration produced by the excavation (Fig. 19a). The minimum effective principal stress (Fig. 19c) around the gallery wall exhibits tensile stress because of the undrained nature of the rock. Thus, the tensile stress and stress concentration around the niche led to the development of the EDZ around the niche. In Step 4, the concentration of the maximum effective principal stress was mitigated (Fig. 19b) and the tensile stress of the minimum effective principal stress decreased due to the drainage behaviour of the rock (Fig. 19d). These stress conditions mitigate rock mass failure, and the EDZ does not develop after excavation. In a previous study, it was reported that the measured hydraulic conductivity in and around the EDZ has not changed for 2 years since the excavation (Aoyagi et al. 2017). This result indicates that the EDZ does not further develop after excavation, and this trend is consistent with that of the stress analysis.

Fig. 19
figure 19

Distribution of the maximum and minimum effective principal stresses around Niche No. 3. a and b Maximum effective principal stress after the installation of the support material (Step 3) and after the drainage analysis (Step 4), respectively. c and d Minimum effective principal stress after the installation of the support material (Step 3) and after the drainage analysis (Step 4), respectively

Full size image

Figure 20 compares the measured and modelled convergence. In the analysis, the deformation of shotcrete at the site corresponding to the location of the measured point (i.e., along the spring line indicated in Fig. 8a) in Niche No. 3 was applied for comparison with the measured convergence. As the excavation analysis was conducted in 2D, the change in convergence with progressive excavation could not simulated. However, the model result (1.97 mm) is similar to the measured maximum convergence (1.6–2.7 mm; Sect. 3.2).

Fig. 20
figure 20

Comparison of the modelled and measured convergence results

Full size image

Figure 21 shows the pore pressure distribution around the gallery in Step 4. The pore pressure decreased over a wide area due to the drained nature of the rock around the EDZ. The pore pressure profiles in the sidewall and below the floor in boreholes H2-1 and H3-2 are shown in Fig. 22, along with measured results for comparison. The trends of pore pressure reduction obtained from numerical analysis after excavation of the niche are generally consistent with the measured results. The modelled results of the extent of the EDZ, convergences, and pore pressure distribution around the gallery are also consistent with the in situ surveys. As such, the numerical analysis reliably models the observed results.

Fig. 21
figure 21

Pore pressure distribution after the drainage analysis (Step 4)

Full size image
Fig. 22
figure 22

Pore pressure profiles estimated from model analysis and measured in the hydraulic test sections along the sidewall (a) and floor (b) of the niche (Step 4)

Full size image

5 Discussion: Estimation of the Highest Potential Hydraulic Conductivity

Ishii (2015) proposed that spatial–temporal predictions of the highest potential transmissivities of fractures in natural fault zones could be made on the basis of mapping or predicting the spatial distribution of the MSI. The MSI is the ratio of the effective mean stress \(\left( {{{\sigma ^{\prime}}_{\text{m}}}} \right)\) to the tensile strength of intact rock \(\left( {{\sigma _t}} \right)\), as follows:

$${\text{MSI}}\,=\,{\raise0.7ex\hbox{${\left( {\frac{{{\sigma _1}+{\sigma _3}}}{2} - p} \right)}$} \!\mathord{\left/ {\vphantom {{\left( {\frac{{{\sigma _1}+{\sigma _3}}}{2} - p} \right)} {{\sigma _t}}}}\right.\kern-0pt}\!\lower0.7ex\hbox{${{\sigma _t}}$}}\,=\,\frac{{{{\sigma ^{\prime}}_{\text{m}}}}}{{{\sigma _{\text{t}}}}}.$$
(2)

Ishii (2015) analysed the MSI values and transmissivities of flow anomalies in fault zones at six sites (Horonobe URL, Wellenberg in Switzerland, Forsmark in Sweden, Olkiluoto in Finland, northern Switzerland, and Sellafield in UK; Fig. 2a), and then formulated the relationship between MSI and transmissivity (T) considering the fault zone heterogeneity as a standard error, as follows:

$$\log T\,=\, - \,3.51\log {\text{MSI}}\, - \,6.54\,\,\,\left( {{\text{standard error }}={\text{ 1}}.{\text{25 in log}}T} \right).$$
(3)

This empirical power law relationship (MSI model) enables estimation of the highest potential transmissivity of fractures in a fault zone, with maximum errors of about two orders of magnitude, when taking into account fault zone heterogeneity (Fig. 2). This relationship is also applicable to estimation of the highest potential transmissivity of fractures in discrete shear fracture systems (Ishii 2017b). In particular, there is a strong correlation between transmissivity of flow anomalies in fault zones and MSI in mudstone formations, as indicated in Fig. 2b, which enhances the reliability of the model for estimating the highest potential hydraulic conductivity of the EDZ in the mudstone formation.

In this study, we assumed that the transmissivity derived from Eq. (3) is the transmissivity of one fracture detected by the BTV image in the EDZ. The estimation of the highest potential hydraulic conductivity in the EDZ was made on the basis of the numerical analysis, as described in Sect. 4, whereas the MSI value was calculated on the basis of Eq. (2). The effective mean stress can be calculated from the effective stress distribution (Fig. 19c, d). The tensile strength was obtained from Brazilian tension tests. The average and standard deviation of the tensile strength were 1.83 and 0.37 MPa, respectively (Table 1). Thus, we set the tensile strength in the MSI calculations to between 1.46 and 2.20 MPa, considering the range in tensile strength of the siliceous mudstone.

We investigated the MSI values in the analysis mesh that corresponds to the locations of the EDZ detected by the BTV survey. The estimated transmissivity was assumed to be the sum of the transmissivity of fractures in the EDZ within the hydraulic test section, as follows:

$${T_e}\,=\,\sum\limits_{i} {{T_i}} \left| {_{{x\,=\,{x_i}}}\,=\,\sum\limits_{i} {{{10}^{\left( { - \,3.51\,\log \,{\text{MS}}{{\text{I}}_i}\left| {_{{x\,=\,{x_i}}} - \,6.54} \right.} \right)}}} } \right.,$$
(4)

where xi is the fracture location, MSIi is the MSI value in the analysed mesh corresponding to the fracture location, Ti is the estimated transmissivity at the fracture site, and Te is the transmissivity in the hydraulic test section expressed as the sum of Ti. The estimated highest potential hydraulic conductivity in the hydraulic test section (ke) was determined by dividing the estimated transmissivity over the length of the test section (L), as follows:

$${k_e}=\frac{{{T_e}}}{L}.$$
(5)

The estimated hydraulic conductivities in each hydraulic test section in the EDZ, as obtained from Eq. (5), were then compared with the in situ test data. The estimated maximum and minimum hydraulic conductivities were calculated considering the standard error of the transmissivity, as described in Eq. (3).

The estimated ranges of the highest potential hydraulic conductivity corresponding to the hydraulic test sections in the EDZ are shown along with the hydraulic test results in Fig. 23. The estimated range is 1–3 orders of magnitude, as the MSI model considers the intrinsic heterogeneity of hydraulic conductivity among fractures (Fig. 2). The estimated ranges of hydraulic conductivities overlap with the measured results on both the sidewall and below the floor of the niche.

Fig. 23
figure 23

Comparison of the range of the highest potential hydraulic conductivities estimated from the MSI model with the hydraulic conductivities measured in the EDZ

Full size image

Multiple factors such as excavation method, excavation geometry, support pattern, distribution of natural (pre-existing) fractures, and stress field can affect the failure behaviour of the EDZ. In the Horonobe URL, mechanical excavation was performed, and the effect of blasting damage can be ignored. Furthermore, there were no or few natural fractures according to the fracture map on the sidewall of the niche and the distribution of fractures obtained from BTV images. This suggests that the effect of natural fractures can be ignored. Thus, for the conditions that the rock mass is relatively homogeneous and artificial damage such as that from blasting can be ignored, the MSI model is likely to be applicable in estimating the spatial distribution of the highest potential hydraulic conductivity around the facility.

In this analysis, the fracture locations based on the BTV surveys were used to estimate the highest potential hydraulic conductivity in the EDZ. For this methodology to be applied to the construction and safe excavation of a facility, the spatial distribution of fractures around the tunnel needs to be estimated in advance. This requires the future development of techniques to estimate the fracture density in the EDZ. In the present study, the number of hydraulic conductivity determinations was relatively small. Therefore, we will in future examine the relationship between transmissivities/hydraulic conductivities in the EDZ using MSI values under different MSI conditions, to improve the reliability of our methodology. These are important avenues for future research.

6 Conclusions

We have developed a method for estimating the hydro-mechanical characteristics of the EDZ induced around a gallery in the Horonobe URL, Japan. The hydraulic conductivity, pore pressure, and extent of the EDZ were investigated using in situ surveys and a hydro-mechanical coupling analysis. We demonstrated the applicability of the method to estimating the highest potential hydraulic conductivity in the EDZ based on the MSI model.

BTV surveys, rock core analysis, and hydraulic tests were performed on the sidewall and below the floor of Niche No. 3 in the Horonobe URL. The in situ tests and surveys showed that the EDZ extent in the niche sidewall and below the niche floor were ≤ 0.6 and ≤ 1.6 m, respectively. The hydraulic conductivity in the EDZ is 2–5 orders of magnitude higher than that around the EDZ. The numerical analysis yielded the EDZ extent of 0.9 m in the sidewall, 0.3 m above the crown, and 1.3 m below the floor. These model results are broadly consistent with the in situ survey data. Modelled pore pressure changes due to the excavation of the niche are also generally consistent with the measured results. Therefore, we conclude that the model results are consistent with the observed results.

The highest potential hydraulic conductivity in the EDZ was estimated on the basis of the empirical relationship between MSI and transmissivity (i.e., the MSI model). The estimated ranges of the highest potential hydraulic conductivities are comparable with the measured results. If the rock mass is relatively homogeneous and artificial damage such as blasting can be ignored, then the MSI model is likely to be applicable in estimating the spatial distribution of hydraulic conductivity around the facility. Our methodology potentially provides a novel means of assessing radionuclide migration around an HLW disposal site.