Numerical simulations of nuclide migration in highly fractured rock masses by the unified pipe-network method

Computers and Geotechnics 111 (2019) 261–276

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Numerical simulations of nuclide migration in highly fractured rock masses by the unified pipe-network method

T

⁎

Guowei Maa,b, Tuo Lib, , Yang Wangb, Yun Chenb a b

School of Civil and Transportation Engineering, Hebei University of Technology, Tianjin 300401, China School of Civil, Environmental and Mining Engineering, The University of Western Australia, Perth, WA 6009, Australia

A R T I C LE I N FO

A B S T R A C T

Keywords: Nuclide migration Fractured rock masses Unified pipe-network method Reactive advection-dispersion effects China’s R&D plan

On account of the numerous discontinuities and multiple physics involved, the numerical evaluation of pollution propagation is computationally sophisticated and time-consuming. This paper proposes a new approach based on the unified pipe-network method to simulate nuclide migration in highly fractured rock masses. This method explicitly represents fracture networks with interconnected nodes and pipes in a spatial domain and applies a set of universal constitutive models to both fractures and rock matrices. Whether fractures or rock matrix are the dominators of mass transport depends on the connectivity of fracture networks. Minor changes in connectivity may trigger significant percolation threshold effects.

1. Introduction Nuclear power is regarded as an alternative clean source of energy that can alleviate fuel shortages and climate change. In China, the installed capacity of nuclear power plants (NNP) will reach 48.6 GW by 2020 [1]. Such massive amounts of electricity generated by NNPs will result in 8,011.8 tons of spent nuclear fuel (SNF). The disposal of SNF and other radioactive waste has become a serious obstacle in the development of nuclear energy. This is evidenced by the fact that radioelements, such as neptunium, plutonium, americium, and technetium, are highly radioactive, toxic, heating-generating, and have long halflives. At present, deep geological disposal hundreds of meters underground is an internationally recognized method for the permanent repository of HLW [2], such as Onkalo in Finlan [3], the Gorleben laboratory in Germany [4], and Yucca Mountain in the U.S. [5]. Multiple layers of atificial and natural barriers have been engineered to contain waste, and host rock is considered the final and permanent contaminant vessel. However, discontinuities in typical geological formations provide preferential pathways where dissolved contaminants can migrate rapidly [6-12]. In terms of aquifer preservation, it is important to understand and quantify the mass transport of radioactive solutes in fractured rock masses prior to site selection. There are three challenges when evaluating the propagation of contamination: (1) the lack of available data from site investigations, (2) difficulties in the representation of numerous discontinuities and arbitrary boundary conditions, and (3) the coupled effects of multiple physics.

⁎

Uncertainty in the geometric and hydraulic properties compromises the accuracy of seepage simulation in rock masses. Fracture networks are commonly described by several statistical characteristics, including orientation, trace length, aperture, and location. These statistical characteristics are represented with distribution models, whose parameters are evaluated from fracture outcrops shown on excavation surfaces or the scanline and sampling windows. Previous researches suggest that the orientation of fracture set obeys fisher distributions [1315]. The cumulative size distribution of trace length follows either power law or lognormal distributions [12,16,17]. And the generation of fracture centers adopts either multiplicative cascade [18,19] or Poisson processes [20,21]. Based on an identical distribution model, infinite realizations of fracture networks can be reconstructed. To diminish randomness in individual simulations, a practical path is generating a series of stochastic realizations based on data collected and conducting statistical analysis with the results [22-24], which is commonly referred to as the Monte Carlo method. The representation of discontinuities can be categorized into two major types: the continuum and the discrete [25,26]. In the continuum framework, fracture networks are represented in less detail and implicitly as equivalent porous mediums. Parameters, such as hydraulic permeabilities and dispersion coefficients, are homogenized across the entire simulation domain. They are either regarded as a part of the single continuum [27-29] or as a separate continuum overlapping the rock matrix [30,31]. However, as researchers have emphasized [6,32], continuum approaches are theoretically valid only with the existence of

Corresponding author. E-mail address: [email protected] (T. Li).

https://doi.org/10.1016/j.compgeo.2019.03.024 Received 10 October 2018; Received in revised form 16 January 2019; Accepted 27 March 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

Nomenclature

piξ , pjξ

Abbreviations

t t1/2 uξ D∗ Dξ

DFN discrete fracture networks FEM finite element method HLW high-level radioactive waste LAD linear advection-diffusion NNP nuclear power plant REV representative elementary volume R&D Plan long-term research and development plan for geological disposal of high-level radioactive waste in China SNF spent nuclear fuel UPM unified pipe-network method

Dijξ

GL Jξ Kξ K ijξ Qijξ Rξ Viξ αLξ

List of Symbols

b cξ

ciξ , c jξ

c¯ij kξ ṁijξ n oc1 ec2 n c1 e

half of fracture aperture solute concentration in fluid, and the superscript ξ denotes the type of medium, ξ = f for fracture, and ξ = m for rock matrix nodal concentration at vertices i and j in the corresponding medium average concentration in pipeij permeability of corresponding medium total mass transport rate in pipe ij in the corresponding medium unit normal vector of plane Aoc1 ec2 in the tetrahedral element unit normal vector of plane l c1 e in the triangular element

μb

σb

λ τξ

ϕξ ϕiξ

nodal pressure at vertices i and j in the corresponding medium time variable half-life of radioactive matters velocity of fluid in the corresponding medium molecular diffusion coefficient in solvent hydrodynamic dispersion coefficient of corresponding medium resultant dispersion coefficient of pipe ij in the corresponding medium radiation intensity of solution dispersion mass flux in the corresponding medium distribution coefficient of corresponding medium resultant conductance of pipe ij in the corresponding medium fluid flow rate in the pipe ij in the corresponding medium retardation coefficient of corresponding medium volume controlled by i in the corresponding medium dispersivity of corresponding medium in the longitudinal direction mean value of hydraulic apertures in a fracture networks system standard deviation of hydraulic apertures in a fracture networks system decay constant of radioactive matters inverse of tortuosity of the flow channel in the corresponding medium the porosity of corresponding medium average porosity of corresponding medium in the volume controlled by node i

Eulerian scheme may accumulate errors to the next time step, it can produce “continuous” concentration profiles without significant statistical errors. With the concentration field, it is convenient to incorporate concentration-dependent chemical processes, like non-linear equilibrium sorption, reactions between different chemical species, or decay reactions [45]. The mass transport of radioactive effluent is a typical multi-physics process. To focus on the accuracy and efficiency of numerical models, the following assumptions are commonly adopted to idealize the scenarios considered: fluids are isothermal and incompressible; the solution is diluted so variations in concentration have minor effects on hydraulic properties; chemical reactions that could modify fracture openings, such as precipitation and dissolution, do not occur. Hence, the processes considered in fractured rock masses include: (1) advection and dispersion in the fractures and rock matrix, (2) sorption reactions on the surface of fracture walls and the voids in the rock matrix, (3) decay of radioactive matter, and (4) mass exchanges at the fracture intersections and on fracture-matrix interfaces. Efforts have been made to develop sophisticated constitutive models simulating the above processes. Fluid flow in fracture networks is likely non-Darcian with high in-situ stress [46,47] or with rough fracture walls [48]. Dispersion may not be Fickian when considering the effects of porosity, pore size distribution, and pore connectivity [49]. Numerous sorption mechanisms have been reported as first- or second-order kinetic processes [50]. In this study, the constitutive models were simplified to compare the model’s accuracy with analytical solutions for traditional linear advection-diffusion (LAD) problems [51,52]. It was assumed that fluid velocity obeys Darcy’s law and dispersion is Fickian. Dispersion in the transverse direction within each fracture ensures complete mixing across its width. Sorption reactions follow linear instantaneous equilibriums. This paper proposes a new path to solve LAD problems with the unified pipe-network method (UPM) [20]. Continuum models that use

representative elementary volume (REV), above which the properties of the medium are considered to be statistically homogeneous [33]. This assumption is less applicable to fractured rock masses since natural fracture networks often exhibit multi-scale or fractal organizations [10,12,21,34-36]. In these cases, it is difficult to identify the homogenization scale of the interested strata. Notwithstanding if REV existed, using homogenized parameters in the simulation of radioactive spreading would likely result in tragic consequences, since they tend to over-estimate the arrival time of pollutants and under-estimate solution concentrations at monitoring locations. In contrast, discrete methods are more logical and less abstract than continuum models [37]. The spacial configuration of fractures is explicitly represented in a heterogeneous and anisotropic geological environment for accurate predictions of contaminant transport [38,39]. Regardless of the model used to represent the fractured porous media, there are two schemes to numerically simulate the advectiondispersion transport: the Lagrangian and the Eulerian approaches. A good example of Lagrangian models is the random walk particle tracking method [40-43]. This method moves each mass parcel through the porous medium using the velocity field by advection and a random displacement by hydraulic dispersion. Although it avoids solving the transport equation directly, this approach still requires a large amount of parallel computing. Hassan and Mohamed [44] demonstrated that in a homogeneous aquifer with 20,000 cells approximately 2,500,000 particles were needed to achieve the smoothness of concentration contour. Not only a computational burden but also parameterizations, especially when concerning discontinuous properties and highly distorted grids, are the principal drawbacks that limit its application in field-scale problems. In comparison, the alternative is desirable for highly heterogeneous media and coupled multi-physics effects. Eulerian models void statistical errors and fulfill the local solute mass balance easily, especially in subsurface transport where discontinuities enhance the material heterogeneity. Although time-dependent modeling in the 262

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

2. Mathematical model

homogenized parameters to represent the effect of discontinuities sacrifice accuracy excessively, while conventional discrete approaches struggle with calculating complexity. Comparatively, with geometric equivalence, the proposed method improves computational efficiency and elevates the processing capability to thousands of fractures. As a method within the framework of discrete models, it well retains result precision. This method adopts the concept of node and pipe similar to the particle tracking approaches, and hence the entire fracture-matrix system is discretized into unified nodes and pipes with conforming meshes [53]. However, as a typical Eulerian model, UPM has advantages in the simulation of coupled multiphysics effects in a highly fractured porous media. The unity of discretized pipes is firstly reflected at the universal constitutive models that are applicable over the entire heterogeneous domain. It also avoids the constraint of mass exchange on the interface because of the consistency of pressure and concentration at the conforming interface nodes. And hence, the distribution of solute concentration can be described conveniently. Moreover, the virtual pipes are unified geometrically among various physical fields. Fluid flow, mass flux, and heat flux are simulated in a universal computer module, and the coupled effects are treated with greater flexibility and efficiency. The suggested method has been validated in simulations of oil production [54,55] and geothermal exploitation [56,57]. This paper aims to present an algorithm that enables the simulation and analysis of nuclide solute transport in complex fracture networks. Section 2 describes the mathematical model regarding UPM. Section 3 verifies the algorithm with analytical solutions to cases containing a single fracture and a set of parallel fractures. In Section 4, a segment of the HLW repository is simulated and analyzed to perform a safety assessment during a nuclear leak. Section 5 continues the discussion if this case. Finally, Section 6 delivers a brief conclusion and thoughts on future work.

2.1. Governing equation As reported in the literature, the governing LAD equations for an incompressible diluted solution are expressed as below, when the source and sink term are neglected.

∂ϕ ξ c ξ 1 + λϕ ξ c ξ = ξ [∇ ·(D ξ ∇c ξ − uξ c ξ )] ∂t R

(1)

where the superscript ξ represents the type of medium (i.e., ξ = f for fracture and ξ = m for rock matrix); ϕ [−] is the porosity of the volume; c [M·L−3] is the solute concentration; t [T ] is the time variable; and u [L·T −1] is the flow velocity. In contrast to conventional LAD problems, mass transport of radioactive solutions is also involved in reactions including radiation and adsorption. The radiation in spent fuel rods and other radioactive waste is defined by the decay constant λ [T −1], which yields

λ=

ln2 t1/2

(2)

where t1/2 [T ] is the half-life of the radioactive matter. Moreover, the effect of sorption on the fracture surface and the rock solid are represented by retardation factors, which yields

Rf = 1 + Rm = 1 +

Kf b Km

(3)

Kf

[L] is the distribution where b [L] is half of the fracture aperture; coefficient of the fracture, which is defined as the mass of solute adsorbed per unit area of the fracture surface divided by the concentration in solution; and K m [−] is the distribution coefficient in the rock matrix, defined as the mass of solute adsorbed in a unit of bulk volume of rock matrix divided by the solution concentration [58]. These coefficients are considered as constants when the adsorbed solute is assumed to be a linear proportion of the concentration. The hydrodynamic dispersion coefficient D [L2 ·T −1] combines the effects of mechanical dispersion and

Fig. 1. Typical discretized elements in rock matrix and fracture plane. 263

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

Qijm = ∫A oc

molecular diffusion, yielding

1 ec 2

D ξ = ϕ ξ (αLξ uξ + τ ξ D )

ṁ ijm = ∫A oc

(4)

1 ec 2

2.2. Unified pipe-network method

K ijm =

In UPM, the fracture-matrix domain is discretized into triangles and tetrahedrons, which represent fracture and rock matrix elements, respectively. Triangular and tetrahedral elements have conforming vertices on the fracture-matrix interface. All vertices are considered as representative nodes with controlled volume, and the interconnecting edges are regarded as equivalent pipes for mass transport between nodes. To determine the parameters of the discretized system, a FEMlike linear shape function was adopted [20]. Therefore, the distribution of fluid pressure and solute concentration within elements could be approximated with the linear interpolation of nodal pressures and concentrations. The fundamental concept of UPM is a geometric equivalence. Fig. 1 illustrates a typical discretized element in the rock matrix and its conforming fracture element, which are represented as a tetrahedron and a triangle, respectively. Point o is the circumcenter of the tetrahedron ijmn ; points c1, c2 , c3 , and c4 are the circumcenters of its triangular surfaces; and points e , f , g , h , r , and s are the midpoints of edges. Therefore, the tetrahedron is divided into four volumes by the crosssection boundaries oc1 ec2 , oc1 rc4 , oc1 sc3 , oc2 fc4 , oc4 hc3 , and oc3 gc2 . All scale variables, such as pressure and concentration, are assigned to these four controlled volumes represented by the vertices i , j , m , and n , respectively. Similarly, the conforming triangle ijm is a meshed element of the fracture plane. It is also divided into three quadrilaterals iec1 r , jec1 s , and mrc1 s . With individual apertures assigned to each quadrilateral as its thickness, the resulting volumes, represented by the vertices’ nodes i , j , and m , are derived in the fracture. Mass transport between two control volumes is equivalent to that in the pipe’s interconnected representative nodes. Take pipe ij as an example. The mass transport between the controlled volumes of nodes i and j consists of two parts: the mass flux in the rock matrix from volume iec1 rc2 oc4 f to volume jec1 sc2 oc3 g , and the mass flux in the fracture from volume iec1 r to volume jec1 s . Accordingly, the Darcian flow and mass flux in the fracture can be expressed as below

c1 e

f

f

+ cj 2

l c1 e

kf

2b μ

Dijf =

Aoc1 ec 2 k m lij μ

Dijm =

lij

∂ (ϕi f Vi f + ϕim Vim ) ci ∂t fi

=

∑

(6)

l c1 e lij

2bD f

Aoc1 ec 2 lij

Dm

(7)

[L2],

k m [L2],

+

+ λ (ϕi f Vi f + ϕim Vim ) ci

1 ⎡ f ⎛ ci + c wf ⎞ f f ⎤ Qiw ⎜ ⎟ − Diw (ci − c w ) ⎥ 2 Rf ⎢ ⎝ ⎠ ⎣ ⎦

w=1 mi

∑ v=1 f

1 ⎡ m ⎛ ci + c wm ⎞ Qiv − Divm (ci − c vm) ⎤ ⎥ Rm ⎢ ⎝ 2 ⎠ ⎦ ⎣ ⎜

[L3]

⎟

(8)

Vim [L3]

where Vi and are the controlled volumes represented by node i in the fracture (if exists) and in the rock matrix, respectively; fi and mi denote the total number of connected pipes in the fracture (if existing) and rock matrix. Consequently, a complex 3-D fracture-matrix domain is equivalent to a simple network of nodes and pipes, which can be solved readily with two-point flux approximation. 3. Verifications 3.1. Domain with a single fracture A fracture-matrix domain containing a single fracture was verified with the proposed method and compared with an analytical solution [51]. Table 1 lists the values of parameters. The original analytical solution was derived in two dimensions with infinite dimension in the longitudinal and transverse directions of the fracture, while the solute concentration converged to zero at 7 m along the fracture and 1.2 m into the rock matrix. Therefore, given the model thickness at 10 m, the size of the simulation domain was assumed to be 2.4 m × 10 m × 7.0 m, as shown in Fig. 2(a). A single fracture was placed vertically in the middle of a rock block, and the origin of the Table 1 Input parameters for a fracture-matrix domain containing a single fracture.

⎞

⎟

⎠

⎞ ⎠

The parameters for element characteristics D f [L2 ·T −1], and Dm [L2 ·T −1] account for saturation and porosity, especially when considering unsaturated or two-phase fluid flow. With the assumption of pressure and concentration consistency at the interface (i.e., pim = pi f and cim = ci f ), substituting Eqs. (5) and (6) into (1) would obtain a unified form of the discretized governing equation for all nodes.

(J f − u f c¯ijf )·n c1 e dl

ci = Dijf (ci f − c jf ) − Qijf ⎛⎜ ⎝

2

kf

u f ·n c1 e dl = K ijf (pi f − pjf )

ṁ ijf = 2b ∫l

cim + c m j

where Aoc1 ec2 [L2] and n oc1 ec2 [−] are the area and unit normal vector of cross-section boundary oc1 ec2 , respectively. Thus, the equivalent coefficients yield

K ijf =

c1 e

(J m − umc¯ijm)·n oc1 ec2 dA

= Dijm (cim − c jm) − Qijm ⎛ ⎝

where αL [L] is dispersivity in the longitudinal direction of fluid flow; u [L·T −1] is the absolute value of intrinsic velocity; τ [−] is the inverse of tortuosity of the flow channel; and D∗ [L2 ·T −1] is the molecular diffusion coefficient in the corresponding solvent, which was groundwater in this study.

Qijf = 2b ∫l

um·n oc1 ec2 dA = K ijm (pim − pjm )

(5)

where l c1 e [L] and n c1 e [−] are respectively the length and the unit normal vector of the cross-section boundary c1 e ; u [L·T −1] and J [M·L−2 ·T −1] are respectively the average fluid velocity and diffusion flux across the intersection boundary; pi [M·L·T −2] and pj [M·L·T −2] are the nodal pressures; ci [M·L−3] and cj [M·L−3] are the nodal concentrations; and Kij and Dij are respectively the equivalent conductance and dispersion coefficients of pipe ij . With the same approximation applied to the rock matrix, the mass transport in the rock matrix can be derived in a unified form:

Parameter

Symbols

Values

Units

Fracture aperture Fracture porosity

2b ϕf

100 1

μm −

Fracture tortuosity

τf

1

−

Dispersivity in fracture

aLf

0.5

m

Retardation coefficient in fracture

Rf

1

−

Matrix porosity Matrix tortuosity Dispersivity in matrix

ϕm τm aLm Rm D∗ t1/2

0.01 0.1 0

− −

Retardation coefficient in matrix Diffusion coefficient Half-life

264

1

m −

1.6 × 10−5 12.35

cm2/ s year

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

Fig. 2. 3-D fracture-matrix domain with single fracture: (a) geometric configuration of simulation domain; and contours of concentration at (b) 100 days; (c) 1,000 days; and (d) steady state.

is the breakthrough curve (i.e., concentration versus time) at a given location? Fig. 3(a) shows concentration profiles along the fracture at 100 and 1000 days, as well as in the steady state that occurs approximately at 10000 days. Excellent agreement between the numerical solution of UPM and the analytical solution was obtained for each recorded location at the middle of the fracture along the z axis. If fluid flows and first-order reactions were both overlooked, the concentration would reduce proportionally with distance from the pollution source, and the concentration profile would be constant in gradient. The applied advection would enhance solute motion and push the tracer front towards the fluid exit, while radioactive decay and sorption reaction would reduce the concentration downstream. Their coupled effects reshape the concentration curves as demonstrated. Based on the analytical solution, we derived that the location with 0.01 relative concentration is at 6.66 m, which is defined as 0.01 penetration depth (i.e., d 0.01) [51]. Fig. 3(b) shows the breakthrough curves for three monitoring locations at 1/2 , 1/4 , and 1/10 of d 0.01, respectively. The intensity of mass transport declines with time, since the concentration gap narrows as the pollutant spreads. The declination of mass transport rate

system O was located on the top of the fracture. A constant hydraulic head is applied at the ends of the fracture, and hence the flow velocity along the longitudinal direction of fracture remained 0.01 m/day, which is assumed saturated and incompressible. In contrast, the rock matrix remains impermeable. The initial and boundary conditions of concentration could be described as

c (t = 0) = 0

c (x = 0, z = 0) = c0

c (x = 0, z = 7 m) = 0 c (x = ± 1.2 m) = 0 where c0 [M·L−3] is a constant solute concentration at the pollution source. Fig. 2(b)–(d) demonstrates the contour of concentration along the fracture at 100 days, 1,000 days, and steady state. To describe the results in detail, the questions of immediate interest are therefore: how will the contaminant concentration distribute at a given time, and what 265

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

Fig. 4. Results for fracture-matrix domain with parallel fractures: (a) concentration profiles for the fracture at 100 days , 1000 days , 10000 days , and steady state; (b) breakthrough curves for the fracture at z = 1/2 , 1/4 , and 1/10 of.d 0.01.

Fig. 3. Results for fracture-matrix domain with single fracture: (a) concentration profiles for the fracture at 100 days, 1000 days , and steady state; (b) breakthrough curves for the fracture at z = 1/2 , 1/4 , and 1/10 of.d 0.01.

becomes increasingly faster at 1/10 of d 0.01, where the equilibrium of concentration is achieved sooner. Perfect agreements were again obtained for the process of concentration accumulation at all locations. Hence, UPM is not only useful for predicting the concentration distribution, but also for forecasting the reaching time when the radiation quantity at certain locations reaches the alert level.

concentration profiles along the fracture at 100, 1000, and 10000 days, as well as the steady state that occurs at approximately 20000 days. Owing to faster flow velocity, the downstream concentration is elevated, which pushes the penetration depth of tracer further. The 0.01 penetration depth d 0.01 in this case is 72 m. Although adsorption in rock matrix and radioactive decay are still effective at reducing concentration, the enhanced flowrate requires a thick layer of host rock for pollution isolation. Results demonstrate excellent agreement between UPM and analytical solutions at these four recorded times. Agreements were also achieved in the breakthrough curves for three monitoring locations at 1/2 , 1/4 , and 1/10 of d 0.01. When modelling fracture-matrix systems, rock matrix is not only the storage of the solute, but also the mass transport channels. The single fracture model can be considered a parallel fracture system with infinite spacing. As the spacing of two neighboring fracture decreases, the penetration depth of contaminant increases, since the effect of solute exchange through rock matrix is enhanced. The back-flux from the matrix is even more significant when considering advection flow in rock matrix in which complex fracture networks are embedded.

3.2. Domain with parallel fractures Another classic case verified was a fracture-matrix domain containing parallel fractures. Sudicky and Frind [52] provided analytical solutions and examined the influence of spacing between neighboring fractures. The parameters used were identical to those listed in Table 1, except that the longitudinal dispersity (i.e., αLf ) drops down to 0.1m and the flow velocity in fracture increases to 0.1 m/day. Rock matrix in this case is still assumed to be impermeable. When the fracture spacing is 0.5 m, the negligible concentration occurs at 80 m downstream. Similar to the previous case, the size of the simulation domain was set as 0.5 m × 10 m × 80 m. The initial and boundary conditions were listed as

c (t = 0) = 0

3.3. Robustness analysis

c (x = ± 2.5 m , z = 0) = c0

Two convergence tests with different element sizes and different time steps were carried out to analyze the robustness and feasibility of this method. The total grid number is an essential factor that may change the computational complexity and accuracy of results in

c (x = ± 2.5 m , z = 80 m) = 0 Fig. 4 plots the results of these cases. Fig. 4(a) indicates the 266

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

after pollution intrusion. As shown in Fig. 5(b), notwithstanding only one step takes place, the resulting concentrations in rock matrix are still very close to the analytical solutions. When the time steps decrease to 1, 000 or 300 days, the resulting deviation from analytical solutions becomes negligible. Since the suggested method provides implicit solutions, the results obtained are highly insensitive to the size of the time step. 4. Applications 4.1. Study area and statistical model Since 1985, the site selection of China’s first HLW repository has undergone nationwide screening, regional screening, and area screening [59]. At this stage, the Beishan area has been approved as “the priority area” for candidate disposal sites by the China Atomic Energy Authority and the Ministry of Environment Protection. Among eight granite intrusions in this area, Jiujing, Xinchang-Xiangyangshan, and Yemaquan are considered to have the most potential. Before final site confirmation, studies have been made on regional crust stability, tectonic evolution, lithological investigations, and hydrogeology. The preliminary geophysical survey was conducted through 11 deep boreholes and 8 shallow boreholes (BS01-BS19) in these three primary locations [59,60]. Site investigations show that the Beishan area possesses good integrity and favorable engineering conditions for an underground HWL repository in terms of geology, hydrogeology, and hydrology. This area is exceptionally inhabitable, has poor economic potential, and has convenient transportation conditions. The crust is stable with seismic intensity lower than Grade VI, and no moderate or significant earthquake has been recorded in history [59]. The average annual precipitation is only 70 mm, while the evaporation is over 3000 mm. Thus, groundwater resources are poor, and the water table is at 28-46 m subsurface [59]. Permeability tests suggest that the hydraulic conductivity of rock masses in the section of interest (i.e., 400 m-600 m below the ground surface) is extremely low, ranging between 3.9 × 10−14 m / s and 6.6 × 10−10 m / s . Researchers investigated statistical models of discontinuities in the nearly intact near-field granite and the fractured far-field granite. In this paper, our example is a homogeneous region in the Jiujing area, which is surrounded by the Shiyuejing and the Jiujing fault zones. Fractures in this region were categorized into three groups based on traces detected on the outcrop. Table 2 lists the distribution parameters for the configuration of fracture networks [61]. Accordingly, an infinite number of DFNs could be generated, and the results from different realizations vary to some extent. However, as the construction of an HLW repository proceeded, the precise fracture configuration was gradually revealed in the near-field. Therefore, a practical path to use the distribution pattern in the generation of DFN is the semi-deterministic approach, by which the fracture system consists of both deterministic fractures that are detectable on outcrops or with other noncontact methods and indeterminate fractures that are generated randomly with distribution models. Four seepage tests were also conducted on site, and the transmissibility of fractures is in the order of 10−7 m / s . Although the permeability of the fracture plane is commonly anisotropic, the fracture transmissibility was assumed to be 2.3 × 10−7 m / s in all directions in this study.

Fig. 5. Convergence tests with: (a) different mesh sizes; and (b) different time steps. Table 2 Input parameters for the generation of DFN. Parameter[unit ]

Set I

Set II

Set III

Distribution pattern of centers

Poisson 0.028

Poisson 0.024

Poisson 0.017

Fisher 48.6 61.0 6.79 Log-normal 1.42 1.55

Fisher 101.0 74.0 2.38 Log-normal 2.03 1.42

Fisher 316.0 68.0 11.8 Exponential 1.51 −

Number density [1/ m3] Distribution pattern of orientation Dip direction [deg ] Dip angle [deg ] Fisher constant [−] Distribution pattern of size Expectation [m] Standard deviation [m]

numerical simulations. Three tests, respectively with 11,896, 46,292, and 126,606 tetrahedral elements, were employed in the single fracture case to examine the sensitivity of mesh size. Concentration in rock matrix at 1/10 of d 0.01 (i.e., z = 0.66 m and y = 5 m) was measured at a steady state. As shown in Fig. 5(a), when the domain is discretized into 11, 896 elements, there exists some deviation between the analytical solution and that from UPM. However, the gap narrows rapidly when the number of tetrahedral elements rises to 46, 292 and 126, 606. We can conclude that the UPM can achieve stable results with a wide range of discretization size. The time step, alternatively, is another factor that contributes to errors in the transient state. We examined relative concentrations at the same monitoring locations (i.e., z = 0.66 m and y = 5 m) 3,000 days

4.2. Simulation domain and boundary conditions A section of a backfilled tunnel with a vertical burial hole in an HWL repository was investigated in this study, as shown in Fig. 6(a). Since the R&D plan is still in the stage of site selection, the geometric configuration of the artificial protection structure was assumed to be the design used in the Onkalo project [3,62]. The spacing of two consecutive tunnels is 25.0 m, and the separation of two holes is 11.0 m. 267

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

Fig. 6. Simulation domain: (a) 3-D geometric configuration of backfilled tunnel and disposal hole; (b) dimensions and boundary conditions on the 2-D cross section.

Fig. 7. Simulation domain with fracture networks: (a) geometric configuration of fracture networks; (b) discretization of simulation domain.

However, the access tunnel is backfilled with bentonite that can also be considered impermeable and non-diffusive. The proposed method specializes in a fracture-matrix domain with deterministic discontinuities, while the determinacy of fractures depends on their scale and the scale of the simulation domain. In this study, one realization of fracture networks was generated based on the distribution models from site statistics. Although statistical methods were required to eliminate the uncertainty induced by the lack of geometric and hydraulic data, we focused on the pollutant penetration with an explicitly represented fractured domain. In this realization, three sets of fractures were generated based on Table 2. The spatial

The tunnel is excavated to 3.5 m in width and 4.35 m in height, and the disposal hole is 1.75 m in diameter and 8.2 m in depth. Therefore, the host rock was assumed to be a 25 m × 11 m × 37.2 m fractured block. The tunnel center is located at the center of the simulation block. Fig. 6(b) illustrates the details of dimensions. A pressure gradient of 8 Pa/m was assumed along the x axis, while the rest of boundaries (i.e., y = ± 5.5 m and z = ± 18.6 m) were set to be impermeable. The impermeable boundaries were also assumed to be non-diffusive, while the concentration gradient on the pressure boundaries converges to zero. In the failure of the artificial barrier, the disposal hole was regarded as a source of pollution with a constant concentration c0 . 268

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

Fig. 8. Contour of solute concentration at: (a) 1000 days ; (b) 10, 000 days ; (c) 20, 000 days ; (d) steady state.

generated with 236, 347 nodes, 1, 427, 720 tetrahedron elements (i.e., rock matrix), and 190, 297 triangle elements (i.e., fractures), as shown in Fig. 7(b). Most of the hydraulic properties were set as shown in Table 1, while the retardation coefficients for both fracture and matrix were set at 1.2 . In addition, the permeabilities in the fracture and rock matrix were assumed to be 2.3 × 10−10 m2 and 1 × 10−15 m2 , respectively, according to the site tests.

Table 3 Half-lives of radioactive matters. Elements

Symbols

Values

Units

Tritium

3H

12.35

year

Plutonium − 241 Curium − 244

241Pu

year

244Cm

13.2 17.6

Californium − 249

249Cf

360

year

Americium − 241

241Am

432.2

year

Plutonium − 251

251Pu

898

year

year

4.3. Results of nuclide migration Fig. 8 demonstrates the resultant concentration. The concentration of pollutant decreases exponentially with distance from the source. Due to the advective flow along the x axis, mass transport is enhanced in

configuration of 706 fractures is illustrated in Fig. 7(a). With a selfprogrammed discretization algorithm, conforming meshes were 269

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

Fig. 9. Breakthrough curves of average concentration and radioactivity on the downstream boundary for: (a) nuclide matters with short half-lives; (b) nuclide matters with long half-lives.

this direction, and hence the estimation of the maximum hydraulic head in the designed working life is essential for scenario analysis. Additionally, discontinuities providing dominant conduits for pollution penetration deteriorate the pollution level downstream. One of the most crucial merits of discrete methods in numerical simulation is enabling the deterministic pattern of concentration distribution. At 1, 000 days after the failure of artificial protection, the most radioactive matter is still contained in the vicinity of the burial hole, and some pollutants may penetrate to the half-way point from the downstream boundary (i.e., x = 18.75 m) along fractures that intersect with the source. If the pollution source has a concentration of 1 mol/L (i.e., c0 = 1 mol/L), the average extent of pollution in the downstream fluid is as low as 1.88 × 10−22 mol/ L at this stage. At 10, 000 days after the nuclear leak, a significant amount of pollutant reaches the downstream boundary. The average pollutant concentration climbs exponentially to 1.65 × 10−6 mol/ L and continues expanding into rock matrix over the next 10, 000 days, when the concentration reaches 1.20 × 10−5 mol/ L downstream. Eventually, the system achieves a steady state at approximately60, 000 days. At this stage, the flowrate on the downstream boundary is 6.18 × 10−6 L/ s , and the total mass emitted on the downstream boundary is 2.25 × 10−10 mol/ s . Consequently, the average concentration of fluid flow is 3.64 × 10−5 mol/ L , and the intensity of radioactivity in this fluid is then measured by

GL = cλ where

GL [T −1·L−1]

Fig. 10. Effects of variance in fracture apertures: frequency of downstream concentrations with (a) μb = 0.1 mm ; (b) μb = 0.2 mm ; (c) μb = 0.5 mm .

volume of fluid. The radioactivity of fluid emitted downstream is 3.90 × 1010 Bq/ L when the pollutant is tritium. Therefore, the effectiveness of host rock highly depends on the decay constant of radioactivity. Relatively low permeability and apparent diffusion coefficient in the rock matrix slows down the speed of pollutant spreading, while radioactive decay continuously reduces the total amount of solute. Hence, there exists an upper limit of the half-life under which a certain thickness of host rock can contain radioactivity below the alert level.

(9) is the number of decayed atoms in unit time and unit 270

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

33.6%. The concentrations converge to a value that is defined by the non-radioactive solute. Contradicting to previous cases, solutes with longer half-lives tend to have lower radiation intensity. In addition, the time required to reach a steady state is more than 200 days for this group of radioactive matter, since the solution enables deeper penetration into rock matrix, thus taking longer to reach equilibrium. 5.2. Effect of fracture apertures Besides the half-lives of radioelements, statistical characteristics of fracture distribution models can also influence the safety and effectiveness of the host rock. This study overlooks variable openings in a single fracture but uses a universal hydraulic aperture to represent its overall conductivity. Apertures in various fractures are assumed to be normally distributed [63], which are characterized by the mean value, μb , and standard deviation, σb . With the identical geometric configuration of fracture networks shown in Fig. 7, we conduct 9 groups of numerical experiments to evaluate the effects of aperture variations in mean values and standard deviations. The mean values of hydraulic apertures are set at 0.10 mm, 0.20 mm, and 0.50 mm, while standard deviations are 10%, 50%, and 100% for each mean value. In each group of experiments, 100 realizations are simulated to determine downstream concentrations and diminish the randomness of results. Fig. 10 illustrates frequency histograms of relative concentrations at downstream boundary considering variance in mean values and standard deviations of hydraulic apertures. When the mean aperture is 0.10 mm, and the standard deviation is 0.01 mm, the expectation of downstream relative concentration after 100 realizations is converging to 3.55 × 10−5 . When the standard deviations increase up to 0.05 mm and 0.10 mm, the expectation values decrease down to 3.04 × 10−5 and 2.01 × 10−5 , respectively. Similar phenomena occur as the mean values reaching 0.20 mm and 0.50 mm. Expectations of relative concentration on the downstream boundary are as high as 1.12 × 10−4 and 3.28 × 10−4 , when the standard deviation is still low. However, with increase in standard deviations, the expectation values are sliding away. Apparently, increasement in standard deviation of hydraulic apertures will lead to more divergent resulting downstream concentration among different realizations. More importantly, the diversity of fracture apertures more likely results in lower downstream concentrations. The most probable cause of this phenomenon is that the overall conductivity in a complex fracture system is dominated by the minimum aperture along the strategic paths.

Fig. 11. Effects of variance in geometric configurations: frequency of downstream concentrations.

5. Discussion 5.1. Effects of half-lives The half-life of radioactive matter, which ranges from days (e.g., to millions of years (e.g., 235U ), has significant impacts on the effectiveness of natural barriers. Matter with a smaller half-life has a smaller penetration depth, since more active radiation reactions lead to faster lessening of concentration. Long-life radioelements may have higher remaining concentrations during spreading, while their relatively low intensity of radiation is likely less dangerous. Therefore, the type of radioelements that are taken into consideration depends on the scale of the simulation domain. In this study, radioelements with halflives of both tens of years and hundreds of years were evaluated and compared, as listed in Table 3. The concentrations and radiation for elements with half-lives on the order of days pose no threats to the effectiveness of repositories, while elements with half-lives of more than thousands of years will act as non-radioactive elements, which should be considered in larger-scale domains. Fig. 9(a) and (b) show the results of solute concentration and radiation intensity on the downstream boundary when the half-lives of radioelements are in the range of tens of years to hundreds of years, respectively. In the cases with half-lives of tens of years, concentration and radioactivity build up at an exponential rate until approximately 30, 000 days after the nuclear outbreak, and transit into a steady state at about 60, 000 days. Assuming c0 = 1 mol/L at the pollution source, the order of magnitude of average concentrations on the downstream fluid is between 10−5 mol/ L and 10−4 mol/ L , and the radioactivity is of the order of 1010 Bq/ L and 1011 Bq/ L . In solutions with a longer-half-life elements, the fluid concentration on the downstream boundary is higher, since the mass loss caused by the radioactive decay is less significant while the flowrate is identical. Although longer half-lives reduce the rate of radiation reactions, radioactivity on the downstream boundary is still higher when compared with radioactivity caused by elements with shorter half-lives. Among this group of tests, radioactive decay has a strong influence on the ultimate concentration. The half-life of 244Cm is 42.5% longer than that of tritium, while the resulting concentration is 3.65 times higher. However, when the half-lives of nuclide matters rise to hundreds of years, the resulting concentrations are elevated significantly to the order of 10−2 mol/ L , since the nuclide matters tend to perform as non-radioactive elements in which decay barely happens during the pollutant’s propagation from the source to the boundary. Note that elements with longer half-lives result in less significant changes in the ultimate concentrations. The half-life of 251Pu is 2.5 times that of 249Cf , while the resulting concentrations only differ by 131I )

5.3. Effects of geometric configurations Compared with variance in fracture apertures, the randomness of geometric configuration will cause even more significant divergence in results. Using the identical distribution model listed in Table 2, statistical analysis is conducted over 100 realizations of randomly generated fracture networks. Results show that the maximum and minimum downstream concentration can differentiate up to six orders of magnitude. In reality, the spacial configuration of fractures in the designed area becomes increasingly determinate as a project proceeds, especially for those intercepted with artificial structures (i.e., holes and tunnels), which have outcrops on the excavation surface. Therefore, in a semirandom model, a combination of determinate and indeterminate discontinuities was utilized when computationally reconstructing fracture networks [64,65]. Statistical analysis is also conducted and compared given that the intercepted fractures are known. Fig. 11 compares the histograms of resulting concentrations in simulations with completely random fracture networks and semi-random fracture networks. Apparently, the relative concentrations span in a large range from 5.93 × 10−8 to 0.01, and hence they are represented in the form of logarithm. When the fracture networks are generated randomly, the expectation value of relative concentration is 1.47 × 10−5 among 100 realizations. The expectation value in cases with semi271

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

Fig. 12. Typical fracture networks with length ratio at: (a) 0.5 ; (b) 0.8 ; (c) 1.2 ; and (d) 1.5.

random fracture networks is slightly smaller at 4.52 × 10−6 . What even more remarkable is that the standard deviation reduces from 1.55 × 10−4 to 3.18 × 10−5 when the fractures intersecting with holes and tunnels are determinate. The result variance among the semi-random models is still so significant that the difference between the minimum and maximum concentration is three orders of magnitude. This can be explained by the connectivity of fracture networks. It is emphasized that fractures provide the preferential path for pollutant propagation, and they substantially reform the hydraulic traveling distance between two locations. The uncertainty of fracture connectivity between the pollutant source and the downstream boundary is the major contributor to the variance in resulting concentration. When these two locations are connected by fractures, the downstream concentration is exponentially larger compared with less connected cases, since their hydraulic distance is significantly shortened.

5.4. Effect of fracture length Given constant fracture density, the connectivity of fracture networks is highly dependent on the fracture length. Since the service life of a disposal site is commonly over a century, and the designed lifecycle of a natural geological barrier is required to last an ice age (i.e., 250,000 years), a fracture system could hardly remain intact. The mechanism of crack growth is sophisticated. Previous researches made some efforts in two-dimensional fracturing processes [66,67], while three-dimensional simulations of cracking is still inconclusive. Hence, this study focuses on the effects of fracture length without considering causes that may alter fracture size. Fig. 12 shows typical geometric configurations when the mean value and standard deviation of fracture size in Table 2 are factored by a length ratio of 0.5, 0.8, 1.2 , and 1.5. At 0.5, fractures in the porous media are more likely isolated, and the connectivity of fracture networks is poor. As the length ratio increases 272

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

connectivity of the fracture system reaches the percolation threshold. A time-dependent analysis is also conducted in a typical realization with length ratio of 1.2 . Fig. 14 demonstrates the contour of concentration at different stages of solute propagation. Compared with the original domain, the solute travels a longer distance within the first 1, 000 days, since better connectivity of fracture networks provides a more preferred path for mass transport. The average extent of pollution in the downstream fluid rises to 2.76 × 10−10 mol/ L . At 10, 000 days after the nuclear outbreak, the polluted area on the downstream boundary is significantly larger than in the former case. The resulting concentration at 20, 000 days is also apparently elevated, not only along the fracture path but also in the rock matrix. As a result, the penetration depth into the rock matrix is much higher when the system achieves the steady state. The ultimate concentration and radioactivity enlarge by more than ten times reaching to 5.62 × 10−4 mol/ L and 6.02 × 1011 Bq/ L , respectively. More importantly, the arrival times of each pollution level are earlier. The concentration and radioactivity resulting from the nuclide matters 3H , 241Pu , and 244Cm were also compared in the fracture developed system. Fig. 15 shows the resulting concentrations with time. Owing to higher fluid velocity and larger mass flux, concentration and radioactivity on the downstream boundary increase by an order of magnitude. The increase of concentration in the tritium solution is more significant than that of the 244Cm solution, and the difference in radiation intensity is narrow. Consequently, the effectiveness of natural buffer fails. 6. Conclusion At present, computational constraints on the total number of fractures processed in a numerical simulation are still a significant obstacle for the applicability of discrete approaches to practical cases. Considering both computing complexity and the accuracy of results, UPM modeling demonstrates its capability in dealing with reactive mass transport problems in complicated fracture-matrix systems under sophisticated boundary conditions. This study involved the following work:

Fig. 13. Effects of variance in fracture lengths: (a) frequency of downstream concentrations with different length ratios; (b) mean values of downstream concentrations with different length ratios.

• We analyzed challenges in the numerical simulation and concluded

to 0.8, some fractures start to connect and cluster. When the length ratio is 1.2 , fractures are moderately connected, and the pollutant source and the downstream boundary are certainly connected by discontinuities. With the length ratio reaching 1.5, the porous media simulated are highly fractured, mass transport is completely dominated by fractures. Similar to previous experiments, statistical analysis with 100 realizations are conducted in each length ratio to compare the effects of fracture sizes. Fig. 13(a) illustrates the frequency histogram of downstream concentration. Apparently, increase in fracture size enhance the conductivity of the system exponentially. When fractures are isolated, mass transport is dominated by rock matrix, and relative concentration at downstream is 1.27 × 10−9 approximately. Despite the variance in fracture geometric configurations, the standard deviation in results stays minimal. When the length ratio rises to 0.8 and 1.0 , resulting concentrations become extremely diverse due to the uncertainty of connectivity between pollutant source and downstream boundary. As fracture networks are moderately connected (i.e., length ratio= 1.2 ), results of relative concentrations start to be increasingly certain again. The expectation value after 100 realizations is 6.07 × 10−4 . In the highly fractured rock masses (i.e., length ratio= 1.5), the expectation value rises exponentially to 3.97 × 10−3 . At this stage, the dominant contributor of mass transport transfers from rock matrix to fracture networks. The largest deviation occurs at the length ratio of 0.8 and 1.0 when the connectivity of fracture networks is at the transition stage. Therefore, a fractured porous media shows significant percolation threshold effects. As illustrated in Fig. 13(b), slightly enlarged fracture sizes may increase domain conductivity dramatically, once the

•

• •

• 273

that discrete models have advantages in solving problems with strong discontinuity accurately due to the explicit representation of fractures. UPM was proposed to solve reactive LAD problems. The fundamental idea in this model is using geometry equivalence and universal constitutive models to simplify complex three-dimensional problems into one-dimensional ones. Two metric cases with analytical solutions were studied and compared to examine the accuracy of the suggested method. Comparatively, result scatters with different fracture configurations and parameters achieved perfect agreement with the corresponding analytical solutions. We emphasize that rock matrix not only works as a safe buffer by absorbing contaminants but also by serving as transport channels. A scenario assessment was conducted, according to China’s R&D project, in which the fracture networks and artificial structures were generated based on the site statistics. Both transient and steady states were simulated to demonstrate the applicability of UPM to a complex DFN system with hundreds of highly connected fractures. The effects of radioelement half-lives were assessed. When the halflives are of the order of tens of years, the final concentration and radiation intensity of solute are both positively related to its halflife. However, when the half-lives increase to hundreds of years, radioactive solutions tend to perform non-radioactively. The ultimate concentration converges to an upper limit value, and the solute radioactivity is inversely correlated with its half-life. Monte Carlo analysis was conducted to assess the effects of variance in fracture hydraulic apertures. It concludes that the overall

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

Fig. 14. Contour of solute concentration with length ratio of 1.2 at: (a) 1000 days ; (b) 10, 000 days ; (c) 20, 000 days ; (d) steady state.

• •

conductivity of mass transport is dominated by the minimum aperture along the preferential paths. Statistical analysis also compares downstream concentrations resulting from random realizations and semi-random realizations. With some determinate fractures, the variance in results is narrowed, but still as large as three orders of magnitude, especially when the connectivity of fracture networks is uncertain. Lastly, the most critical factor of network connectivity is fracture sizes. Fracture networks show significant percolation threshold effects in the scenario assessment. When fracture sizes are enlarged by 20%, both fluid flowrate and mass flux are elevated by at least an order of magnitude, and the host rock may not be able to isolate the pollutant.

When selecting repository sites for HLW disposal, UPM enables forecasting of the ultimate depth of a contaminant propagation in explicitly represented fracture networks and predicting the time taken to reach this penetration. Eventually, the resultant concentration and intensity of radiation can be utilized in the assessment of protection effectiveness. Preferential paths consisting of fracture networks shorten the mass transport distance between two geometrically distant

Fig. 15. Breakthrough curves of average concentration and radioactivity on the downstream boundary with length ratio of.1.2 .

274

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

locations and hence deform the concentration field. The question of parameterization remains crucial but inconclusive, as numerical simulations are highly dependent on available data. The predicted dangerous zone over an alert level of solute concentration and radiation intensity might vary due to the uncertainty of fracture distribution. However, UPM can also be incorporated into the established regulations of statistical analysis for pollution zones. In subsequent papers, we will investigate the R&D project in more detail by considering thermal and multi-phase effects.

2001;23:247–58. https://doi.org/10.1016/S0191-8141(00)00094-8. [23] Belayneh MW, Matthai SK, Blunt MJ, Rogers SF. Comparison of deterministic with stochastic fracture models in water-flooding numerical simulations. Am Assoc Pet Geol Bull 2009;93:1633–48. https://doi.org/10.1306/07220909031. [24] Chesnaux R, Allen DM, Jenni S. Regional fracture network permeability using outcrop scale measurements. Eng Geol 2009;108:259–71. https://doi.org/10.1016/ J.ENGGEO.2009.06.024. [25] Maryška J, Severýn O, Vohralík M. Numerical simulation of fracture flow with a mixed-hybrid FEM stochastic discrete fracture network model. Comput Geosci 2005;8:217–34. https://doi.org/10.1007/s10596-005-0152-3. [26] Larsson M, Odén M, Niemi A, Neretnieks I, Tsang C-F. A new approach to account for fracture aperture variability when modeling solute transport in fracture networks. Water Resour Res 2013;49:2241–52. https://doi.org/10.1002/wrcr.20130. [27] Oda M. An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses. Water Resour Res 1986;22:1845–56. https://doi.org/10.1029/ WR022i013p01845. [28] Samardzioska T, Popov V. Numerical comparison of the equivalent continuum, nonhomogeneous and dual porosity models for flow and transport in fractured porous media. Adv Water Resour 2005;28:235–55. https://doi.org/10.1016/J. ADVWATRES.2004.11.002. [29] Pan J-B, Lee C-C, Lee C-H, Yeh H-F, Lin H-I. Application of fracture network model with crack permeability tensor on flow and transport in fractured rock. Eng Geol 2010;116:166–77. https://doi.org/10.1016/J.ENGGEO.2010.08.007. [30] Guo B, Tao Y, Bandilla K, Procedia MC-E, 2017 undefined. Vertically Integrated Dual-porosity and Dual-permeability Models for CO2 Sequestration in Fractured Geological Formation. UArizonaEdu n.d. [31] Jerbi C, Fourno A, Noetinger B, Delay F. A new estimation of equivalent matrix block sizes in fractured media with two-phase flow applications in dual porosity models. J Hydrol 2017;548:508–23. https://doi.org/10.1016/J.JHYDROL.2017.03. 028. [32] Long JCS, Remer JS, Wilson CR, Witherspoon PA. Porous media equivalents for networks of discontinuous fractures. Water Resour Res 1982;18:645–58. https:// doi.org/10.1029/WR018i003p00645. [33] Bear J. Dynamics of fluids in porous media. Dover; 1988. [34] Bodin J, Porel G, Delay F, Ubertosi F, Bernard S, de Dreuzy J-R. Simulation and analysis of solute transport in 2D fracture/pipe networks: The SOLFRAC program. J Contam Hydrol 2007;89:1–28. https://doi.org/10.1016/J.JCONHYD.2006.07.005. [35] Liu R, Yu L, Jiang Y. Fractal analysis of directional permeability of gas shale fracture networks: a numerical study. J Nat Gas Sci Eng 2016;33:1330–41. https://doi.org/ 10.1016/J.JNGSE.2016.05.043. [36] Li B, Liu R, Jiang Y. A multiple fractal model for estimating permeability of dualporosity media. J Hydrol 2016;540:659–69. https://doi.org/10.1016/J.JHYDROL. 2016.06.059. [37] Selroos J-O, Walker DD, Ström A, Gylling B, Follin S. Comparison of alternative modelling approaches for groundwater flow in fractured rock. J Hydrol 2002;257:174–88. https://doi.org/10.1016/S0022-1694(01)00551-0. [38] Poeter E, Gaylord DR. Influence of aquifer heterogeneity on contaminant transport at the Hanford Site. Ground Water 1990;28:900–9. https://doi.org/10.1111/j. 1745-6584.1990.tb01726.x. [39] Tsang C-F, Neretnieks I. Flow channeling in heterogeneous fractured rocks. Rev Geophys 1998;36:275–98. https://doi.org/10.1029/97RG03319. [40] Banton O, Delay F, Porel G. A new time domain random walk method for solute transport in 1-D heterogeneous media. Ground Water 1997;35:1008–13. https:// doi.org/10.1111/j.1745-6584.1997.tb00173.x. [41] Delay F, Bodin J. Time domain random walk method to simulate transport by advection-dispersion and matrix diffusion in fracture networks. Geophys Res Lett 2001;28:4051–4. https://doi.org/10.1029/2001GL013698. [42] Bodin J, Delay F, de Marsily G. Solute transport in a single fracture with negligible matrix permeability: 1. fundamental mechanisms. Hydrogeol J 2003;11:418–33. https://doi.org/10.1007/s10040-003-0268-2. [43] Bodin J, Delay F, de Marsily G. Solute transport in a single fracture with negligible matrix permeability: 2. mathematical formalism. Hydrogeol J 2003;11:434–54. https://doi.org/10.1007/s10040-003-0269-1. [44] Hassan AE, Mohamed MM. On using particle tracking methods to simulate transport in single-continuum and dual continua porous media. J Hydrol 2003;275:242–60. https://doi.org/10.1016/S0022-1694(03)00046-5. [45] Salamon P, Fernàndez-Garcia D, Gómez-Hernández JJ. A review and numerical assessment of the random walk particle tracking method. J Contam Hydrol 2006;87:277–305. https://doi.org/10.1016/J.JCONHYD.2006.05.005. [46] Derode B, Cappa F, Guglielmi Y, Rutqvist J. Coupled seismo-hydromechanical monitoring of inelastic effects on injection-induced fracture permeability. Int J Rock Mech Min Sci 2013;61:266–74. https://doi.org/10.1016/J.IJRMMS.2013.03.008. [47] Chen Y-F, Liu M-M, Hu S-H, Zhou C-B. Non-Darcy’s law-based analytical models for data interpretation of high-pressure packer tests in fractured rocks. Eng Geol 2015;199:91–106. https://doi.org/10.1016/J.ENGGEO.2015.10.011. [48] Zhang W, Dai B, Liu Z, Zhou C. A pore-scale numerical model for non-Darcy fluid flow through rough-walled fractures. Comput Geotech 2017;87:139–48. https:// doi.org/10.1016/J.COMPGEO.2017.02.011. [49] Mu D, Liu Z-S, Huang C, Djilali N. Determination of the effective diffusion coefficient in porous media including Knudsen effects. Microfluid Nanofluidics 2008;4:257–60. https://doi.org/10.1007/s10404-007-0182-3. [50] Ho Y, McKay G. Pseudo-second order model for sorption processes. Process Biochem 1999;34:451–65. https://doi.org/10.1016/S0032-9592(98)00112-5. [51] Tang DH, Frind EO, Sudicky EA. Contaminant transport in fractured porous media: Analytical solution for a single fracture. Water Resour Res 1981;17:555–64. https:// doi.org/10.1029/WR017i003p00555.

Acknowledgments This research was supported by the Research Training Program (RTP) scheme at UWA. The authors wish to thank the National Natural Science Foundation of China (NSFC) for their financial support (No.51778029 \& No.51627812). We would like to thank the reviewers for their constructive suggestions and discussion. References [1] Wang Q. China needing a cautious approach to nuclear power strategy. Energy Policy 2009;37:2487–91. https://doi.org/10.1016/J.ENPOL.2009.03.033. [2] Dlouhy Z. Disposal of radioactive wastes. 2009. [3] Aalto P, Aaltonen I, Kemppainen K. Programme for repository host rock characterisation in the ONKALO (ReRoC) 2009. [4] Kickmaier W, McKinley I. A review of research carried out in European rock laboratories. Nucl Eng Des 1997;176:75–81. https://doi.org/10.1016/S00295493(96)01344-1. [5] Walton JC. Effects of evaporation and solute concentration on presence and composition of water in and around the waste package at Yucca Mountain. Waste Manag 1993;13:293–301. https://doi.org/10.1016/0956-053X(93)90057-4. [6] Berkowitz B, Bear J, Braester C. Continuum models for contaminant transport in fractured porous formations. Water Resour Res 1988;24:1225–36. https://doi.org/ 10.1029/WR024i008p01225. [7] Willis-Richards J, Wallroth T. Approaches to the modelling of hdr reservoirs: a review. Geothermics 1995;24:307–32. https://doi.org/10.1016/0375-6505(95) 00011-E. [8] Li G. Tracer mixing at fracture intersections. Environ Geol 2002;42:137–44. https:// doi.org/10.1007/s00254-001-0483-x. [9] Matthäi SK, Belayneh M. Fluid flow partitioning between fractures and a permeable rock matrix. Geophys Res Lett 2004;31:n/a-n/a. https://doi.org/10.1029/ 2003GL019027. [10] Neuman SP. Trends, prospects and challenges in quantifying flow and transport through fractured rocks. Hydrogeol J 2005;13:124–47. https://doi.org/10.1007/ s10040-004-0397-2. [11] Faulkner DR, Jackson CAL, Lunn RJ, Schlische RW, Shipton ZK, Wibberley CAJ, et al. A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones. J Struct Geol 2010;32:1557–75. https://doi. org/10.1016/J.JSG.2010.06.009. [12] Liu R, Jiang Y, Li B, Wang X. A fractal model for characterizing fluid flow in fractured rock masses based on randomly distributed rock fracture networks. Comput Geotech 2015;65:45–55. https://doi.org/10.1016/J.COMPGEO.2014.11. 004. [13] Fisher R. Dispersion on a sphere. Proc R Soc A Math Phys Eng Sci 1953;217:295–305. https://doi.org/10.1098/rspa.1953.0064. [14] Priest SD. Discontinuity Analysis for Rock Engineering. Netherlands: Springer; 1993. [15] Kemeny J, Post R. Estimating three-dimensional rock discontinuity orientation from digital images of fracture traces. Comput Geosci 2003;29:65–77. https://doi.org/ 10.1016/S0098-3004(02)00106-1. [16] Baghbanan A, Jing L. Hydraulic properties of fractured rock masses with correlated fracture length and aperture. Int J Rock Mech Min Sci 2007;44:704–19. https://doi. org/10.1016/J.IJRMMS.2006.11.001. [17] Bertalan Z, Shekhawat A, Sethna JP, Zapperi S. Fracture strength: stress concentration, extreme value statistics, and the fate of the weibull distribution. Phys Rev Appl 2014;2:034008https://doi.org/10.1103/PhysRevApplied. 2.034008. [18] Darcel C, Bour O, Davy P. de Dreuzy JR. Connectivity properties of two-dimensional fracture networks with stochastic fractal correlation. Water Resour Res 2003;39.. https://doi.org/10.1029/2002WR001628. [19] Meakin P. Invasion percolation on substrates with correlated disorder. Phys A Stat Mech Its Appl 1991;173:305–24. https://doi.org/10.1016/0378-4371(91)90366-K. [20] Ren F, Ma G, Wang Y, Li T, Zhu H. Unified pipe network method for simulation of water flow in fractured porous rock. J Hydrol 2017;547. https://doi.org/10.1016/j. jhydrol.2017.01.044. [21] Liu R, Li B, Jiang Y. A fractal model based on a new governing equation of fluid flow in fractures for characterizing hydraulic properties of rock fracture networks. Comput Geotech 2016;75:57–68. https://doi.org/10.1016/J.COMPGEO.2016.01. 025. [22] Mauldon M, Dunne WM, Rohrbaugh MB. Circular scanlines and circular windows: new tools for characterizing the geometry of fracture traces. J Struct Geol

275

Computers and Geotechnics 111 (2019) 261–276

G. Ma, et al.

[52] Sudicky EA, Frind EO. Contaminant transport in fractured porous media: analytical solutions for a system of parallel fractures. Water Resour Res 1982;18:1634–42. https://doi.org/10.1029/WR018i006p01634. [53] Wang Y, Ma G, Ren F, Li T. A constrained Delaunay discretization method for adaptively meshing highly discontinuous geological media. Comput Geosci 2017;109:134–48. https://doi.org/10.1016/J.CAGEO.2017.07.010. [54] Chen Y, Ma G, Li T, Wang Y, Ren F. Simulation of wormhole propagation in fractured carbonate rocks with unified pipe-network method. Comput Geotech 2018;98:58–68. https://doi.org/10.1016/J.COMPGEO.2017.11.009. [55] Ma G, Chen Y, Jin Y, Wang H. Modelling temperature-influenced acidizing process in fractured carbonate rocks. Int J Rock Mech Min Sci 2018;105:73–84. https://doi. org/10.1016/J.IJRMMS.2018.03.019. [56] Chen Y, Ma G, Wang H, Li T. Evaluation of geothermal development in fractured hot dry rock based on three dimensional unified pipe-network method. Appl Therm Eng 2018;136:219–28. https://doi.org/10.1016/J.APPLTHERMALENG.2018.03. 008. [57] Chen Y, Ma G, Wang H. Heat extraction mechanism in a geothermal reservoir with rough-walled fracture networks. Int J Heat Mass Transf 2018;126:1083–93. https:// doi.org/10.1016/J.IJHEATMASSTRANSFER.2018.05.103. [58] Freeze RA, Cherry JA. Groundwater. Prentice-Hall; 1979. [59] Wang J. High-level radioactive waste disposal in China: update 2010. J Rock Mech Geotech Eng 2010;2:1–11. https://doi.org/10.3724/SP.J.1235.2010.00001. [60] Zhao XG, Wang J, Qin XH, Cai M, Su R, He JG, et al. In-situ stress measurements

[61]

[62] [63] [64]

[65] [66]

[67]

276

and regional stress field assessment in the Xinjiang candidate area for China’s HLW disposal. Eng Geol 2015;197:42–56. https://doi.org/10.1016/J.ENGGEO.2015.08. 015. Guo L, Li X. A dissertation submitted to the academic degree committee of nanjing university for applying the doctor degree in geological research on permeability of fracture systems and fault structure in high-level radioactive waste repository Beishan candidate site. Nanjing University; 2015. Toprak E. THM modelling of deposition tunnel in ONKALO project 2013. Hakami E, KTH T& kopiering). Aperture distribution of rock fractures. Tekniska högsk; 1995. Billaux D, Chiles JP, Hestir K, Long J. Three-dimensional statistical modelling of a fractured rock mass—an example from the Fanay-Augères mine. Int J Rock Mech Min Sci Geomech Abstr 1989;26:281–99. https://doi.org/10.1016/0148-9062(89) 91977-3. Ivanova V, Yu X, Veneziano D, Einstein H. Development of stochastic models for fracture systems 1995. Zhang Y, Zhuang X. Cracking elements: a self-propagating strong discontinuity embedded approach for quasi-brittle fracture. Finite Elem Anal Des 2018;144:84–100. https://doi.org/10.1016/J.FINEL.2017.10.007. Zhang Y, Lackner R, Zeiml M, Mang HA. Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations. Comput Methods Appl Mech Eng 2015;287:335–66. https://doi.org/10.1016/J.CMA.2015.02.001.