- Email: [email protected]
Inverse modeling of saturated-unsaturated flow in site-scale fractured rocks using the continuum approach: A case study at Baihetan dam site, Southwest China
Journal Pre-proofs Research papers Inverse modeling of saturated-unsaturated flow in site-scale fractured rocks using the continuum approach: A case study at Baihetan dam site, Southwest China Yi-Feng Chen, Hao Yu, Heng-Zhen Ma, Xing Li, Ran Hu, Zhibing Yang PII: DOI: Reference:
S0022-1694(20)30153-0 https://doi.org/10.1016/j.jhydrol.2020.124693 HYDROL 124693
To appear in:
Journal of Hydrology
Received Date: Revised Date: Accepted Date:
11 October 2019 19 January 2020 13 February 2020
Please cite this article as: Chen, Y-F., Yu, H., Ma, H-Z., Li, X., Hu, R., Yang, Z., Inverse modeling of saturatedunsaturated flow in site-scale fractured rocks using the continuum approach: A case study at Baihetan dam site, Southwest China, Journal of Hydrology (2020), doi: https://doi.org/10.1016/j.jhydrol.2020.124693
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2020 Published by Elsevier B.V.
Inverse modeling of saturated-unsaturated flow in site-scale fractured rocks using the continuum approach: A case study at Baihetan dam site, Southwest China
Yi-Feng Chen *, Hao Yu, Heng-Zhen Ma, Xing Li, Ran Hu *, Zhibing Yang State key Laboratory of Water Resources and Hydropower Engineering Science Wuhan University, Wuhan 430072, China Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering of the Ministry of Education, Wuhan 430072, China
Email: [email protected], [email protected]
Revised Manuscript January 20, 2020
A revised paper submitted for possible publication in Journal of Hydrology
*Corresponding
author. Tel.: +86-27-68774295; fax: +86-27-68774295. Email address: [email protected] (Y.-F. Chen), [email protected] (R. Hu). 1
Abstract
Modeling saturated-unsaturated flow in fractured rock formations remains a challenging issue due to the difficulties in properly calibrating the unsaturated flow properties for fractured rocks. On the basis of the continuum approach, this study uses inverse modeling to determine the unsaturated hydraulic parameters of fractured rocks, loose sediments, and geologic structures in a large-scale slope under difficult hydrogeologic conditions at the Baihetan dam site, Southwest China. The saturated hydraulic conductivities were determined a priori by packer tests. Field data, including time-series of groundwater level observations and available discharge measurements, is used. A combined procedure of orthogonal design, finite element forward analysis, artificial neural network, and genetic algorithm is adopted for the inverse analysis with high computational efficiency. The numerical results well reproduce the typical features of multiple water tables and wedge-shaped unsaturated zones controlled by gently-sloping tuff zones of low vertical permeability, showing 4665% of groundwater flowing subhorizontally along each tuff zone and the remaining 3554% infiltrating vertically and recharging deeper layers. The drastic drawdown of groundwater level after 10-year site characterization is also reproduced. The calculated discharge and water tables agree reasonably well with field observations, indicating that the estimated parameters are representative of the unsaturated hydraulic properties at the site. This work provides an efficient and effective methodology for seeking an optimal solution in large-scale inverse modeling of saturated-unsaturated flow, and further manifests the feasibility of the continuum 2
approach in predicting unsaturated flow in site-scale fractured rock formations. Keywords: Unsaturated flow; Hydraulic properties; Fractured rocks; Inverse modeling; Continuum approach
Nomenclature
Cw
Specific moisture capacity
D
Fractal dimension of flow paths
h*
Maximum allowable values of pressure head
h
Pressure head
h
Prescribed pressure head on the water head boundary h
hd
Minimum pressure head determined by temperature and air relative humidity on the evaporation boundary e
hp
Ponding depth of water on the infiltration boundary i
I(t)
Rainfall intensity as a function of time t
K
Saturated hydraulic conductivity tensor
kr
Relative permeability
M
Number of groundwater observation boreholes
N
Number of weirs
n
Outward unit normal vector to the boundary
Qj
Measured discharge data at weir j
Qj
Simulated discharge data at weir j
qn*
Maximum allowable flux on the third-type boundary
3
qn
Prescribed flux (positive for inflow) on the flux boundary q
Re
Maximum evaporation rate on the evaporation boundary e
Se
Effective saturation
Ss
Specific storage
s
Suction
t
Time
w
Weight coefficient
z
Positive upward vertical coordinate
, m, n
Empirical parameters in the van Genuchten model
r
Compressibility of rock
w
Compressibility of water
Volumetric water content
r
Residual volumetric water content
s
Saturated volumetric water content
e
Evaporation boundary
i
Rain infiltration boundary
s
Seepage face boundary
Sign indicator
Porosity
A step function of h
v
Flow velocity vector
w
Density of water
4
Positive parameter in the modified Mualem model
i
Measured groundwater level data in borehole i
i
Simulated groundwater level data in borehole i
1. Introduction
Modeling the saturated-unsaturated flow in fractured rocks has received considerable attention in the past decades in a wide range of applications such as nuclear waste disposal, contaminant transport, oil and gas extraction, underground oil storage, and landslide mitigation (Rulon et al., 1985; Liu et al., 1998; Li et al., 2017b). Accurately characterizing this flow process is extremely difficult and has been one of the most challenging issues in subsurface hydrology mainly because of the multi-scale geometries and complex interactions of the fracture-matrix system, the technical difficulties in measuring fracture flow in the field, and the multi-dimensional and heterogeneous flow behaviors (Evans and Rasmussen, 1991; Liu et al., 1998). In particular, the groundwater typically moves preferentially through fractures and faults of large scales (Peters and Klavetter, 1988; Illman and Hughson, 2005). Besides, multiple water tables and wedge-shaped unsaturated zones may develop if the stratigraphy consists of subhorizontal layers of contrasting permeabilities (Rulon et al., 1985). Various modeling approaches have been proposed to describe the multi-scale, multi-dimensional flow behaviors through the fracture-matrix system, including the equivalent continuum (Peters and Klavetter, 1988; Finsterle, 2000; Liu et al., 2003;
5
Lu and Kwicklis, 2012), dual-continuum (Illman and Hughson, 2005; Mathias et al., 2006), triple-continuum (Wu et al., 2004), discrete fracture network (Evans and Rasmussen, 1991), percolation network (Kueper and McWhorter, 1992) and stochastic representation approaches (Illman and Hughson, 2005; Lu and Kwicklis, 2012). The equivalent continuum model treats the fractured rock as a porous medium with an equivalent permeability obtained through homogenization of the flow process. This approach has been widely used for its simplicity and high computational efficiency (Bathe et al., 1979; Desai et al., 1983; Zheng et al., 2005; Chen et al., 2008, 2016; Li et al., 2017a). The discrete fracture network model can explicitly simulate flow in all fractures in a rock mass (Evans and Rasmussen, 1991). However, due to the large number of fractures involved, it is often computationally prohibitive to adopt this approach in site-scale applications. The dual-medium model simultaneously considers fast flow in the fractures and slow flow in the porous matrix (Barenblatt et al., 1960; Warren and Root, 1963; Zimmerman et al., 1993; Illman and Hughson, 2005; Mathias et al., 2006), but its application still has great limitations due to the uncertainty of the spatial distribution of the fracture system. The tripe-continuum model introduces an additional continuum (based on the dual continuum model) that requires one more parameter set for small fractures (Wu et al., 2004), which makes it even more difficult to use than the dual-continuum model. The stochastic approach assumes that the local hydraulic properties are realizations of spatially correlated random fields and derives the partial differential equations representing the large-scale flow conditions through averaging the local governing equations over the
6
ensemble of realizations (Yeh et al., 1985a, 1985b, 1985c; Mantoglou and Gelhar, 1987a, 1987b, 1987c). Among these approaches, the continuum model has dominated as the most commonly-used one for its efficiency in addressing large-scale problems because the presence of numerous fractures of different scales causes extreme difficulty in fracture network construction and flow modeling. The performance of a continuum model, however, largely depends on whether the unsaturated hydraulic properties (i.e., water retention curve and unsaturated hydraulic conductivity) of fractured rocks are properly determined and representative of the site conditions. Field measurements of permeability and porosity for unsaturated fractured rocks prove to be feasible by single-hole (e.g., Illman and Neuman, 2000) or cross-hole (Illman and Neuman, 2001, 2003) pneumatic injection tests without modifying the saturation as done by water infiltration tests. This is not the case, however, for direct measurement of unsaturated hydraulic properties representative of field scale behavior due to the complexity and variability of fracture geometries (National Research Council, 1996; Hughson and Yeh, 2000; Zhang and Fredlund, 2003; Wu et al., 2006), even though these measurements could be performed for core-scale samples of individual fracture (Wang and Narasimhan, 1985; Pruess and Tsang, 1990; Reitsma and Kueper, 1994; Persoff and Pruess, 1995; Bertels et al., 2001; Chen and Horne, 2006; Weerakone et al., 2012; Li et al., 2014) and rock matrix (Montazer and Wilson, 1984; Peters and Klavetter, 1988; Faybishenko and Finsterle, 2000; Flint, 2003). Therefore, the classic water retention curves (Brooks and Corey, 1964; van Genuchten, 1980; Fredlund and Xing, 1994) and relative permeability models
7
(Burdine, 1953; Mualem, 1976) developed for unsaturated soils and porous media are often borrowed for describing the unsaturated flow in fractured rocks. The development or extension of the water retention relations for unsaturated fractured rocks is based on the dual- or triple-porosity concept (Peters and Klavetter, 1988; Zhang and Fredlund, 2003; Wu et al., 2004), active fracture concept (Liu et al., 1998), fractal fracture geometry (Guarracino, 2006), and two-dimensional fracture network simulation and upscaling technique (Liu and Bodvarsson, 2001, 2003; Lu and Kwicklis, 2012). Inverse modeling provides a powerful tool for estimation of hydraulic parameters in water retention models from laboratory experiments, field groundwater observations and hydraulic tests (e.g., Kool and Parker, 1988; Liu et al., 1998, 2003; Bandurraga and Bodvarsson, 1999; Wu et al., 2004, 2006; Song et al., 2014; Younes et al., 2018). However, difficulties always arise in the inverse modeling due to the problems of ill-posedness, computational burden, and scales (Zhou et al., 2014; Hughson and Yeh, 2000). The objective of this study is to demonstrate the feasibility of modeling saturated-unsaturated flow in a large-scale rock slope of 5 km long located at the Baihetan dam site, Southwest China, based on the continuum approach and inverse modeling. This slope is characteristic of difficult hydrogeologic conditions, with multiple water levels in natural state due to the presence of gently-sloping, low-vertical permeability tuff zones and drastic depletion of water table after ten years of site characterization due to the excavation of exploratory adits and drilling of boreholes. The unsaturated hydraulic parameters in the water retention
8
relation and relative permeability model are estimated by inverse modeling, with the saturated hydraulic conductivities of rocks determined a priori by packer tests. The objective function is constructed by using the time-series observations of groundwater level and the available measurements of discharge. In order to achieve high computational efficiency, the combined procedure of orthogonal design (OD), finite element (FE) analysis, artificial neural network (ANN) and genetic algorithm (GA) presented by Zhou et al. (2015) is adopted to estimate the unsaturated hydraulic parameters of rocks at the site. The numerical results are compared with field observations to deepen our understanding of the saturated-unsaturated flow in the rock slope subjected to precipitation/evaporation cycles and exploratory adits excavation.
2. Site characterization
2.1. General description
The study site (Jiang et al., 2014) is situated on the left bank of the lower Jinsha River in Southwest China, where the world’s second largest hydropower project (i.e., Baihetan Project) consisting of a double curvature arch dam of 289 m high and two large-scale underground cavern systems is under construction (Fig. 1). The Jinsha River runs through the dam site from south to north (with a water level of about 591 m a.s.l. in dry seasons), forming a V-shaped valley by fast downcutting. The study area is characteristic of a mountainous landscape, with steep gradients of 5080 below the valley shoulder and gentler gradients towards the watershed between the Jinsha River and its tributary on the west side, the Yibu River (Fig. 1). The site lies in 9
the subtropical monsoon climate zone, where over 80% of the annual precipitation occurs during the wet season between May and October. Precipitation increases with elevation, with the measured mean annual rainfalls being 727.2 mm at 635 m a.s.l. near the dam site and 811.0 mm at 1,357 m, respectively. The mean monthly rainfall and evaporation at the study site are shown in Fig. 2. The average annual temperature is between 15 and 21. Surface runoff is infrequent and of short duration, with only two small perennial streams (Nos. 1 and 3 in Fig. 1) in the study area.
2.2. Geological settings
As shown in Fig. 3, the bedrocks outcropping in the study area are basaltic formations originated from multiple magmatic and volcanic eruption episodes. The formations belong to the Emei Mountain group of the upper Permian system (P2), which can be categorized into 4 basalt flow layers (P21P24) of 445470, 290340, 210260 and 85115 m thick, respectively, according to the eruption sequences. The rock layers mainly contain tholeiitic basalt, cryptocrystalline basalt, microcrystalline basalt and amygdaloidal basalt, typically with a layer of tuff developed on the top of each rock unit (Wang et al., 2009). The tuff layers are 0.1–1.75 m thick, and have regional continuity. The tuff is nonwelded, and is composed of vitric fragments and glass shards (40–50%), volcanic ash (30–40%), basalt clasts (5–10%) and plagioclase clasts (<5%) (Wang et al., 2009). Each tuff layer contains a shear zone (hereafter referred to as tuff zone) of 5–40 cm thick. The formations dip gently towards the right bank, with an orientation of N3050E/SE1525. The basalt rocks are overlain
10
by loose sediments (Q3 and Q4) of several to tens of meters thick on gentle hillslopes, and underlain by thick carbonate rocks (P1) deposited during the early Permian period (Fig. 3a). The main structures in the study area consist of subvertical faults (e.g., F17, F129, F133 and F134) and gently-sloping tuff zones (e.g., C2, C3-1 and C3) and shear zones (e.g., LS331, LS3318 and LS3319), as shown in Fig. 3b. Among them, F17 is a large-scale fault of 0.53 m thick, extending over 1400 m and mainly composed of well-cemented tectonic breccia and cataclasite of moderate permeability. Other faults (e.g., F129, F133 and F134) are of smaller sizes (about 0.11 m thick and extending up to hundreds of meters), and bounded by the tuff zones; these faults are also significant to groundwater infiltration as they outcrop and dip subvertically. The tuff zones are developed in the tuff layers, in which fault gouge layers of 1–5 cm thick are developed as a result of tectonic shear between basalt layers (Fig. S1 in Supplementary material), which makes the permeability of the tuff zones highly anisotropic, with a potential permeability contrast of 2 orders of magnitude between the directions parallel and perpendicular to the tuff layer. Except LS331, the shear zones within each basaltic formation are mostly of small scale, and mainly composed of basalt fragments. Besides the geologic structures mentioned above, four groups of critically-oriented fractures are developed in the basalt rocks, with their geometrical properties shown in Fig. 3b and Table S1 in Supplementary material.
11
2.3. Permeability of rocks
A total of 4413 conventional packer tests at 173 boreholes and 669 high pressure packer tests (Chen et al., 2015; Zhou et al., 2018) at 26 boreholes, together with a large number of pumping tests, recovery tests, slug tests, and in-situ seepage tests, were performed at the site to investigate the permeability of the fractured basalt rocks, structures (i.e., faults, tuff zones, and shear zones) and loose sediments. The boreholes were drilled vertically from the ground surface or from the floors of exploratory adits horizontally excavated at different elevations into the slope (Fig. 3b and Fig. 4). The diameter of the boreholes was 76 mm, and the length of the test sections varied between 3.7 and 6.4 m (Chen et al., 2018). The packer test data were interpreted with the Hvorslev (1951) equation. The test results indicate that the permeability of the fractured basalt rocks generally decreases with increasing depth (Fig. S2 in Supplementary material), and can be more correlated to the degree of weathering, the density of tensional joints and the level of in-situ stress, rather than the lithology. Downward from the ground surface, the hydraulic conductivity K can be classified into six categories or permeability zones (PZs, Fig. 3b): high (K 105), high-moderate (105 > K 3106), moderate (3106 > K 106), moderate-low (106 > K 3107), low (3107 > K 107), and very low (K < 107), all in units of m/s. The main structures are mostly of moderate-low to moderate permeability, but the permeability of the unweathered tuff zones and shear zones is highly anisotropic. Table 1 lists the statistics of the hydraulic conductivity of the rocks and sediments at the site. 12
2.4. Groundwater observations
Groundwater in the study area is principally recharged by precipitation, and discharged to streams and the Jinsha River. The groundwater flow and water table configuration in the basaltic formations are primarily controlled by the fracture system and the backbone structures consisting of subvertical faults and gently-sloping tuff zones and shear zones, as shown in Fig. 4. Due to the very low vertical permeability of the tuff zones and shear zones (Table 1), multiple water tables develop in the slope (Fig. 4), a phenomenon that has been frequently observed and well understood in layered hillsides (e.g., Rulon et al., 1985). Observations of groundwater in the slope started from 2001 through boreholes and exploratory adits, and lasted for over ten years. The exploratory adits were excavated horizontally into the slope, from which a few branches parallel to the river were also excavated (Fig. 4). The boreholes were drilled downward from the ground surface or from the floors of the adits. The mean groundwater levels in wet and dry seasons are shown in Fig. 4. The groundwater level in the slope, although influenced locally by fractures, generally increased with ground surface elevation. The groundwater level varied in a larger range (typically 1020 m) between wet and dry seasons in boreholes near the river bank (consistent with the change in river water level) than far away from the river bank. The presence of tuff zones and shear zones significantly influences the distribution and movement of groundwater in the slope. When an exploratory adit was excavated across the tuff zones or shear zones, a significant amount of discharge was 13
generally observed from the upper block of the tuff/shear zones, indicating that these structures are of low vertical permeability and hinder the vertical infiltration. The discharge into each adit decayed due to the advance of neighboring adits, but the total discharge into the whole adit system maintained a stable magnitude and was less influenced by precipitation events. For example, the total discharge during site characterization
from
the
exploratory
adit
PD61
and
its
five
branches
(PD61-1PD61-5), which cut through a series of faults and tuff/shear zones (e.g., F17 and C3), had rather stable measurements around 110 L/min by weirs. Furthermore, when the tuff zones (C3-1 and C3) were penetrated by the boreholes, a sharp drawdown of groundwater level generally occurred. Multiple groundwater levels were observed after the boreholes were isolated with concrete plugs at the intersection with the tuff zones. Table 2 lists the vertical distance between the groundwater levels observed above and below the tuff zones in four typical boreholes (Fig. 5). For example, in borehole ZK1119 a vertical distance of 24.4 m between the upper and lower groundwater levels was observed; similarly in borehole ZK412, a vertical distance of around 23.0 m between the groundwater levels was observed during the wet season. It should be noted, however, that as commonly done in dam engineering, the unsaturated zone was not monitored, even though techniques (e.g., tensiometers and neutron probes) are available for the measurements of flow quantities such as pressure head and water content (e.g., Ireson et al., 2006). This is a major limitation of this study, but motivates us to examine if the unsaturated zone parameters could be
14
reasonably estimated with the observations in the saturated zone.
3. Inverse modeling approach
3.1. Saturated-unsaturated flow model
In the framework of continuum approach (Peters and Klavetter, 1988) where the fractured rocks are treated as equivalent continua, the saturated-unsaturated water flow is governed by the following mass balance equation (e.g., Mao et al., 2011), also known as the Richards (1931) equation:
Cw + Ss
h v 0 t
(1)
where h is the pressure head, Cwh is the specific moisture capacity, is the volumetric water content, Ss is the specific storage, is a step function of h, = 0 for h < 0 (in unsaturated zone) and = 1 for h 0 (in saturated zone), t is time, and v is the flow velocity described by Darcy’s law:
v kr K (h z )
(2)
where K is the saturated hydraulic conductivity tensor, kr is the relative permeability depending on h or , z is the positive upward vertical coordinate, and h+z is the hydraulic head. The water mass balance equation Eq. (1) is subjected to the following initial condition in the domain of interest:
h(t = 0) h0 (in )
(3)
where h0 is the initial distribution of pressure head in , and the following boundary conditions: 15
(1) The water head boundary condition (Dirichlet boundary condition):
h(t ) h (t ) (on h )
(4a)
where h is the prescribed pressure head on the water head boundary h. (2) The flux boundary condition (Neumann boundary condition): qn (t ) v (t ) n qn (t ) (on q )
(4b)
where qn is the prescribed flux (positive for inflow) on the flux boundary q, and n is the outward unit normal vector to the boundary. (3) The unified lateral conditions on the third-type boundaries, including the seepage faces, rain infiltration surfaces and evaporation surfaces (Borsi et al., 2006; Hu et al., 2017): (h h* ) 0; (qn qn* ) 0 (on s i e ) * * (h h )(qn qn ) 0
(4c)
where s, i and e denote the seepage face boundary, rain infiltration surface, and the evaporation surface, respectively; is a sign indicator with = 1 on s and i, and = 1 on e; h* and qn* are the maximum allowable values of pressure head and flux, respectively, with h*0 and qn*0 on s, h*hp and qn*I(t) on i, and h*hd and qn*Re on e; hp is the ponding depth of water on the infiltration boundary i, and I(t) is the rainfall intensity; hd is the minimum pressure head determined by temperature and air relative humidity on the evaporation boundary e, and Re is a negative flux rate representing the maximum evaporation rate in the field. Eq. (4c) states that on the third-type boundaries (s, i and e), the boundary condition is either a prescribed pressure (h=h*) when the flux is no more than the allowable flux rate (qn<qn*), or a prescribed flux rate (qn=qn*) when the pressure is 16
no more than the allowable pressure (h<h*). The presence of the lateral boundary conditions in Eq. (4c) introduces strong nonlinearity of the partial differential equation (PDE) in Eq. (1) for saturated-unsaturated flow problems, which can be rigorously and effectively solved with a numerical procedure based on the parabolic variational inequality (PVI) method (Hu et al., 2017). The discretized PVI algorithm was detailed in Hu et al. (2017), and was implemented in the FE code THYME (Chen et al., 2009). As clarified by site characterization (Fig. 3 and Fig. 4), the groundwater at the site primarily flows along the direction perpendicular to the river (i.e., along the trace of cross-section I-I in Fig. 1), and the flow along the river direction is rather slow. A 2D FE mesh along section I-I was therefore created for saturated-unsaturated flow modeling, which contains 24,782 four-node isoparametric elements and 25,197 nodes, as plotted in Fig. 6. The size of the model domain is 5000 m long and 2190 m high. The topographic and geologic settings, such as the strata, faults, tuff/shear zones and PZs at cross-section I-I (Fig. 1), were well represented. It should be noted that the exploratory adits parallel to the profile (e.g., PD-61, Fig. 4) could not be accurately represented in the 2D mesh, and were modelled with equivalent spacing per unit length (Fig. 6b and Fig. S3 in Supplementary material). The boundary conditions were specified as follows: (1) the lateral boundaries on the mountain side and on the river side were both assumed to be groundwater divides, and thus were prescribed with a no-flow boundary; (2) the riverbed surface was imposed with a water head boundary according to the fluctuation of the river water
17
level; (3) the surfaces of the exploratory adits were taken as a seepage face boundary; (4) the slope surface was prescribed with a precipitation/evaporation boundary with hp=0 m and hd=200 m (corresponding to a saturation of 1618% for loose sediments and tuff zones). According to the site condition, precipitation was assumed to occur in the middle ten days of a month, and the mean monthly precipitation and evaporation rates were taken from Fig. 2; (5) the base of the model was assumed to be impermeable. The initial distribution of pressure head in the saturated zone of the slope (h00) was determined by estimating a water table that best fits the average groundwater level observations in boreholes with a steady-state flow model. In the unsaturated zone (h0<0), h0 was roughly estimated by assuming a vertically linear reduction of saturation to 60% at the ground surface. To eliminate the influence of h0 distribution on groundwater flow, a natural precipitation/evaporation process was first simulated until the groundwater flow in each season became stable. This process typically lasted for 360 months (from May 1971 to May 2001). The site exploration process was then modelled (see Table 5 for the excavation process of the exploratory adit system), lasting for 120 months (from May 2001 to May 2011). The initial time step and maximum time step in each simulation were taken as 7 hours and 7 days, respectively.
3.2. Constitutive relations and hydraulic parameters
The hydraulic properties that need to be determined for numerical modeling of water movement in saturated-unsaturated fractured rocks include the saturated
18
hydraulic conductivity K, the specific storage Ss, the water content and the relative permeability kr. The hydraulic conductivity tensor K (or its isotropic counterpart K) of fractured rocks at the field scale can be reasonably determined without technical difficulty by performing a sufficient number of field hydraulic tests (e.g., Hsieh and Neuman, 1985). Owing to the experiences that have been long accumulated, the determined hydraulic conductivity is generally of high confidence and representative of the field conditions (Table 1). The specific storage of rocks is estimated by Ss =wg(r+w), where r is the compressibility of rock, w is the compressibility of water, is the porosity, and w is the density of water (Table 1 and Table S2 in Supplementary material). The water content relation, also called the water retention curve, is defined as the relationship between the volumetric water content (or the degree of saturation Se) and suction s (i.e., the negative pressure head, s=h). Quite a number of water retention curves have been developed for unsaturated soils and porous media (Brooks and Corey, 1964; van Genuchten, 1980; Fredlund and Xing, 1994), and extended to incorporate the effects of hydraulic hysteresis and soil deformation (e.g., Assouline et al., 1998; Nuth and Laloui, 2008; Tarantino, 2009; Gallipoli, 2012; Zhou et al., 2012; Hu et al., 2013). These models are often borrowed for unsaturated fractured rocks, due to the technical difficulties in the field-scale measurements of water content relation in rocks of much higher heterogeneity. Numerical simulations and validations (Peters and Klavetter, 1988; Liu et al., 1998, 2003; Liu and Bodvarsson, 2001, 2003; Wu et al., 2004) have shown that the classic van Genuchten (1980) model, with a proper
19
parameterization representative of field conditions, is suitable for representing the water content relation in unsaturated fractured rocks, which reads: Se
r n m 1 | h | s r
(5)
where Se is the effective saturation, s and r are the saturated and residual volumetric water contents, respectively, and , m and n are empirical parameters with m=11/n. The relative permeability kr as a function of the degree of saturation Se is another constitutive relation of great importance in modeling the unsaturated flow through porous and fractured media. This relation can generally be obtained by using the water content relation in the classic conceptual models of Burdine (1953) and Mualem (1976). A major concern in the development of the relative permeability function is the tortuosity of the flow paths, which is often correlated to the fractal dimension of flow geometries (e.g., Guarracino, 2006) or the fraction of active fractures (e.g., Liu et al., 1998) to improve the prediction accuracy. Hu et al. (2015) proposed a modified Mualem relation for deformable soils by taking into account the change in pore size distribution and the tortuosity of flow paths, given by m kr 1 1 Se1/ m
2
(6)
where is a positive parameter related to the fractal dimension D of flow paths by
=2D2. A verification of the relation with experimental data (Hu et al., 2015) yields a global optimal estimate of =0.5 or D=1.25 with high confidence for five types of soils of quite different structures (including loess, silty loam, sandy loam and coarse sand). This suggests that similar to other modified van Genuchten or Brooks-Corey
20
permeability models (e.g., Peters and Klavetter, 1988; Liu et al., 1998; Liu and Bodvarsson, 2001; Guarracino, 2006), Eq. (6) may apply to fractured rocks if the parameter is properly determined. The constitutive relations in Eqs. (5) and (6) contain five parameters (i.e., m, ,
, s and r) for each rock unit, among which the saturated and residual volumetric water contents, s and r, can be empirically estimated with the porosity (or the volume fraction of connected fractures and pores) and the void structure of rocks. It has been understood that the seepage field in fractured media is generally less sensitive to the residual saturation r (Liu et al., 2003), and in the most simplified case, one may assume s = and r = 0. The other three parameters (m, and ) are to be representatively determined by inverse modeling, which will be stated in the next section. The values of s and r and the ranges of m, and for 6 rock zones are listed in Table 3.
3.3. The inverse modeling procedure As previously stated, the parameters K (or K), s and r in the unsaturated flow model can be reasonably measured or empirically estimated with engineering accuracy, but difficulty arises when the other parameters (m, and ) are to be measured for fractured rocks. This difficulty can be overcome with inverse modeling techniques by reproducing the dynamics of the seepage field and tuning the parameters to represent the site characteristics. 3.3.1. Objective function
21
Suppose r is the number of materials of contrasting unsaturated hydraulic properties in the domain . A vector P is defined to denote the hydraulic parameters of the media yet to be determined, with P = [m1, 1, 1, m2, 2, 2, ..., mr, r, r]T. In this study, r is chosen to be 6, representing the following six types of media: loose sediments, high-moderate/moderate permeability zones (PZs), moderate-low/low PZs, very low PZ, faults and tuff/shear zones, as listed in Table 3. It is noted here that some PZs and geologic structures are grouped to reduce the number of unknowns and hence the uncertainties and difficulties in the inverse modeling. To estimate the parameters P, an objective function is needed, which ensures a best fit between the field measurements and the numerical results (Carrera et al., 2005; Chen et al., 2016), given by
M (P) i i min f = 2 i 1 i 2
1
N Q (P) Q 2 j j 2 w 2 j 1 Qj 2 2
1
2 2 2
(7)
where M is the number of groundwater observation boreholes, i and i are the time series measurements of groundwater level and the corresponding numerical results in borehole i, respectively, N is the number of weirs, Q j and Qj are the time series measurements of discharge and the numerical results at weir j, respectively, w is a weight coefficient, and ||||2 denote the Euclidean norm of a vector. Eq. (7) utilizes all available time series measurements of both groundwater level and discharge to improve the inverse solution that is inevitably plagued by the problem of non-uniqueness. In this study, the available time series observations of groundwater level in 9 boreholes (Table 4) drilled in the near-bank region of the slope
22
(Fig. 4) and the available discharge measurements in May 2011 at 3 weirs (Table 5) installed in the exploratory adits (Fig. 4) were chosen for construction of the objective function. Note that except ZK929, the other boreholes were located in the neighbouring profiles, with their distances to cross-section I-I (Fig. 3 and 4) being listed in Table 4. With the increase of distance, the data may be less representative of the selected profile, which is the other weakness of this study, as discussed in Section 5. The weight coefficient w is introduced to ensure a balance between the relative errors of the groundwater level and discharge measurements, and a value of w=0.1 was adopted in this study by trials (Table S3 in Supplementary material). 3.3.2. Inverse modeling procedure To efficiently and reliably address the large-scale inverse problem in complex hydrogeologic conditions, this study adopted a combined procedure of orthogonal design (OD), finite element (FE) analysis, artificial neural network (ANN) and genetic algorithm (GA) recently presented in Zhou et al. (2015). The algorithms are given in Figs. S4-S6 in Supplementary material. OD is an effective statistics method to arrange multifactorial experiments (Gong et al., 2008). In this inverse modeling procedure, OD was employed to construct a small set of representative values for the hydraulic parameters P from a pre-defined parameter space. In this study, a total of 18 parameters (i.e., r=6) were to be estimated, with 7 levels specified for each parameter from the ranges empirically given in Table 3. Consequently, an orthogonal table L49(718) was constructed, which produces a set of 49 parameter combinations. These parameter combinations were then used as input parameters for forward modeling of
23
saturated-unsaturated flow (Hu et al., 2017; Chen et al., 2009), yielding the numerical time-series data of groundwater level and discharge at the corresponding observation points and times for each parameter combination. By doing so, only a limited number of time-consuming FE simulations were performed and hence the computational cost was greatly reduced. As a multivariate non-linear mapping and predictive tool that has been successfully applied in hydrology problems (e.g., Garcia and Shigidi, 2006), a back propagation neural network (BPNN) model was then constructed and trained (Fig. S5 in Supplementary material), using the parameter combinations and the prescribed time-dependent precipitation/evaporation boundary conditions (input data) as well as the numerical results (output data). Consequently, an implicit mapping from the parameter space to the unsaturated seepage fields at the observation points was obtained (Zhou et al., 2015). The BPNN model contains four layers, with 19 neurons in the input layer, 22 and 35 neurons in the two hidden layers and 13 neurons in the output layer, as shown in Fig. 7. As usual, the Sigmoid function was selected as the transfer function of the BPNN model. The Levenberg-Marquardt back-propagation algorithm combined with Bayesian regularization (Rafiai et al. 2013) was adopted for training of the model, in order to obtain quick training time and high generalization accuracy. Finally, as an appropriate optimization tool in hydrology problems (e.g., Karpouzos et al., 2001), GA was used to obtain a globally optimal estimate of P such that the objective function, Eq. (7), is minimized and a best fit between the field measurements and the numerical results calculated by the trained BPNN model is
24
ensured. This is achieved by randomly generating an initial population of hydraulic parameter combinations and then iteratively invoking the genetic operations of selection (reproduction), crossover and mutation (Fig. S6 in Supplementary material). In this study, the size of the initial population was set to 100, and the probabilities for crossover and mutation operations were specified to 0.9 and 0.1, respectively. The GA was run at least 10 times, and the best solution was selected as the optimal result.
4. Results
4.1. The results of inverse modeling Table 3 lists the estimated unsaturated hydraulic parameters for rocks in different zones, and the corresponding water retention curves are plotted in Fig. 8. The inverse modeling showed that the unsaturated flow behavior in the slope is most sensitive to the value of α, which controls the minimum threshold pressure at which air starts to penetrate fractured rocks. A larger α results in a higher infiltration rate into the slope and a stronger fluctuation of groundwater level during the precipitation-evaporation cycles. The estimated value of α for loose sediments (several to tens of meters thick) is 0.60 m1, which is comparable to the findings of Perkins et al. (2014), who showed that α is dependent on depth, varying from 0.63 m1 to 65.2 m1 for depth from 39.2 m to 84.3 m. The representative values of α for fractured rocks in different PZs ranges from 0.010 m1 to 0.034 m1, with lower value for rocks in lower PZ (Table 3). This result is comparable to the calibrated values (α=0.0011.63 m1) for rock matrix continuum at Yucca Mountain (Bandurraga and Bodvarsson, 1999; Flint, 2003). For faults and tuff/shear zones, the estimated values of α are 3.03 m1 and 4.65 m1,
25
respectively, again comparable to the calibrated results (α=5.2658.3 m1) for fractures at Yucca Mountain (Bandurraga and Bodvarsson, 1999; Liu and Bodvarsson, 2001), where the α value was found to be highly sensitive to the average aperture of fractures. The parameter m determines the slope of the water retention curve. Here, it is found that increasing m tends to suppress the fluctuation of groundwater level. The estimated values of m for fractured rocks are in a narrow range between 0.31 and 0.43 (corresponding to n=1.451.75), comparable to the calibrated results (m=0.140.43) for rock samples collected from various depths in different hydrogeologic units at Yucca Mountain (Flint, 2003). The representative value of m is moderate for loose sediments (m=0.44), and higher for faults and tuff/shear zones (m=0.630.69). The unsaturated flow process is less sensitive to the parameter β, but a larger β may result in smaller infiltration rate and groundwater level fluctuation. The estimated values of β for different rock zones are in a narrow range between 0.53 and 0.79 (corresponding to D=1.271.40), slightly larger but rather close to the optimal estimate for soils of quite different structures (Hu et al., 2015). This indicates that, as a first approximation, a value of β=0.5 similarly applies to fractured rocks.
4.2. Comparison of field measurements With the estimated parameters given in Table 3 and Fig. 8 as well as the numerical setup given in Section 3.1, simulations of groundwater flow in the slope were performed for 30 years of natural precipitation/evaporation cycles and 10 years of site exploration. Fig. 9 shows the fluctuations of the observed and calculated groundwater levels over tuff zone C3 (Fig. 4) in nine boreholes listed in Table 4. A comparison between the curves shows that the numerical simulation roughly 26
reproduces the groundwater level observations (Fig. 9a and Fig. 9b), except at boreholes ZK37 and ZK112 (Fig. 9c and Fig. 9d) located over 300 m away from the cross-section (Table 4). This difference is attributed to local terrain changes (Fig. 1). If the groundwater level is transformed to groundwater depth (i.e., with the topography changes compensated), then the difference becomes smaller. The predicted groundwater level at ZK1123 (Fig. 9a) is consistent with the rainfall time-series, but seems to have an opposite trend with the measurements. The groundwater level over tuff zone C3 either remained rather stable (Fig. 9a), or exhibited seasonal fluctuations (Fig. 9bd) as influenced by the rainfall infiltration through faults (e.g., F17 and F129, Fig. 4). Even for the latter, the numerical result well captures the trend of groundwater level variation. Simulated water pressure head distributions (in long-term stable state) in wet and dry seasons are shown in Fig. 10. It is interesting to observe that the presence of tuff zones C3 and C3-1 of low vertical permeability caused the distribution of groundwater in multiple layers. This phenomenon was well reproduced in the numerical model. The development of faults (e.g., F17 and F129) facilitated the infiltration of rainfall into the slope, and resulted in the rise of groundwater level over the tuff zones in wet seasons. The infiltrated groundwater discharged to the river either by subhorizontal movement along the tuff zones (for their high in-plane permeability) or by vertical infiltration across the tuff zones. The numerical simulation shows that for the groundwater across the ceiling of C3, about 46% of the water flowed along the tuff layer and 54% further infiltrated downward to recharge the layers below; similarly for C3-1, 65% of the groundwater flowed along it and the rest 35% infiltrated vertically through (see inset in Fig. 4). Consequently, the water table below the tuff zones remained rather stable, except the near-bank region influenced by river water level
27
fluctuations. Table 2 lists the vertical distance between the groundwater levels above and below the tuff zone C3 or C3-1 in four typical boreholes in the wet season, showing that the numerical result is comparable to the field observations. The weathering of tuff zones (C3 and C3-1) makes their permeability isotropic (K=K//) in subzone B (Table 1 and Fig. 6b), which facilitated the subvertical discharge of groundwater in the shallow region of the slope, causing the tuff zones to mostly remain unsaturated and hence the unsaturated zone to be wedge-shaped. As shown in Fig. 11, the water head contours on May 15, 2011, two years after the exploratory adit PD61 and its five branches PD61-1PD61-5 were excavated (Table 5), indicated a marked alteration of the seepage field (compare Fig. 10 and Fig. 11). The evolution of the calculated discharges out of PD61 and its three branches PD61-3PD61-5 is shown in Fig. 12, where the discharges are compared with the available measurements on May 15, 2011 (Table 5). The excavation of the adits led to increased vertical discharge of groundwater across the tuff zones and hence significant drawdown of water table. The discharges fluctuated significantly during the excavation of the adits, but overall tended to decrease to some stable magnitudes (Fig. 12). In May 2011, the simulated discharge from the adit system was about 114 L/min, which is in excellent agreement with the field measurements (Table 5). At that time, field observation showed that the 340 m near-face section of adit PD61 and the branches PD61-1 and PD61-2 were in the unsaturated zone, while the surrounding rocks of the rest section of PD61 and the branches PD61-3PD61-5 remained saturated. The numerical result (see Fig. 11) well reproduces this observation.
5. Discussion This study has demonstrated the feasibility of modeling saturated-unsaturated 28
flow in fractured rock formations with the continuum approach. This relies on proper estimation of the unsaturated hydraulic parameters of rocks in the van Genuchten model and the modified Mualem relation. The estimates were obtained by inverse modeling of saturated-unsaturated flow using the measurements of groundwater levels and discharges collected in the saturated zone. A major concern is whether the saturated flow is sensitive enough to the hydraulic parameters in the unsaturated zone. An additional sensitivity analysis (Table S4 in Supplementary material) confirms the sensitivities of , m and in the unsaturated zone to the input data obtained from the saturated zone. The parameter is most sensitively influenced, followed by m and , respectively, in sensitivity, which is consistent with the findings of Gribb et al. (1996) and Liu et al. (2003); besides, the parameters are more sensitive in the rock zones where the water tables are located. This partially justifies the inverse modeling strategy adopted in this study. It should be noted, however, that the proposed approach contains quite a number of simplifications and assumptions. First, the flow was modeled with a continuum model, where the details of the fracture system, the flow in individual fractures, and the fracture-matrix interaction were not explicitly accounted for. However, given the difficulty in obtaining (and thus the lack of) sufficient fracture network information from site characterization, the continuum model is arguably acceptable for this large-scale problem. This is partially justified by the capability of the continuum model in capturing the changes in groundwater level induced by changes in rainfall intensity during a typical period (Table S5 and Fig. S7 in Supplementary material). In fact, the equivalent continuum approach has been successfully applied in a large number of studies (e.g., Peters and Klavetter, 1988; Doughty, 1999; Wu et al., 1999; Finsterle, 2000; Liu et al., 2003).
29
Second, only the unsaturated hydraulic parameters were estimated by inverse modeling, where the hydraulic conductivities of rocks were directly determined from the field hydraulic tests. We made this choice because saturated hydraulic conductivities estimated from hydraulic tests have been shown to be representative of field scale flow properties of the rock (Gribb et al., 1998; Zhou et al., 2015; Hong et al., 2017). It should be pointed out that the heterogeneity and scale effect of the parameters within each rock zone were not considered. Since the field test data do not support a more detailed delineation of the heterogeneity at the site, it was technically infeasible to account for heterogeneity within each zone in the simulations (Eaton and McCord, 1995). Third, the flow was assumed to be two-dimensional in the slope, where the flow component perpendicular to the cross-section, albeit slow, was neglected. Due to insufficient data in the selected profile, the use of measurements from the neighboring profiles further gave rise to uncertainties in the 2D modeling. The 2D assumption is justifiable, given the fact that the groundwater at the site (the left bank) primarily flows along the direction perpendicular to the river (i.e., along the trace of section I-I in Fig. 1). In addition, the chosen cross-section (I-I in Fig. 1) roughly represents the geological conditions (including the topography, lithology, and geological structures) at the site. Therefore, the 2D model could be considered representative of the site characteristics. The simplification of the model from 3D to 2D is common in seepage analysis due to the constraint of computational resource in inverse calculation of 3D model (Yan and Jiao, 2018; Guo et al., 2015). A 3D model will undoubtedly provide a more realistic and effective evaluation, which is a focus of our future work. Finally, the use of ANN in the inverse modeling procedure significantly reduced the computational cost (Rafiai et al., 2013; Zhou et al., 2015; Hong et al., 2017). As a
30
surrogate model, ANN has become very popular and efficient in engineering problems (Zhang et al., 2017). It is difficult to estimate the unsaturated parameters in large-scale groundwater flow problems without using surrogate models (Asher et al., 2015; Salem and Tomaso, 2018). One should, however, be cautious when using surrogate models in inverse analysis because the reduced computational time possibly comes at the cost of oversmoothing the relation between the input and output data, depending on the choice and size of the training samples. In this study, we used a combination of orthogonal design and forward modeling to create proper training samples for the ANN model. It is nevertheless acknowledged that further improvement is possible in the surrogate modeling. Despite of these limitations, the proposed approach effectively reproduced the multiple water tables and wedge-shaped unsaturated zones in a large-scale rock slope. The change of groundwater levels was predicted within a relatively small error by making full use of site characterization data. The inversed unsaturated hydraulic parameters are comparable to existing results for rocks, fractures, and soils (Bandurraga and Bodvarsson, 1999; Flint, 2003; Hu et al., 2015). Without introducing sophisticated flow models and constitutive relations, the numerical results are of engineering accuracy and useful for, e.g., the design and performance assessment of impervious barriers (Li et al., 2017a; Zhou et al., 2018).
6. Conclusions This study investigated the groundwater flow behavior in a large-scale rock slope by means of inverse modeling using the continuum approach based on field observations spanning ten years at the Baihetan dam site, Southwest China. The inverse modeling showed that the unsaturated hydraulic parameters are in the range of
31
=0.014.65 m1, m=0.310.69 (n=1.453.23) and =0.530.79 for fractured rocks, loose sediments and geologic structures at the site, comparable to existing results for rocks and fractures at Yucca Mountain (Bandurraga and Bodvarsson, 1999; Flint, 2003) and/or soils of different types (Hu et al., 2015). The multiple water levels above and below the tuff zones and the drawdown of water table in the slope during site exploration were reasonably reproduced, and the estimated parameters are overall representative of the unsaturated hydraulic properties of rock formations at the site. It has been challenging to obtain the field scale unsaturated hydraulic parameters in fractured rocks. Here, it is demonstrated that without using the difficult-to-acquire flow and saturation data in the unsaturated zone, the unsaturated hydraulic parameters can be effectively obtained through the groundwater level and discharge data in the saturated zone. It is concluded that the inverse modeling procedure adopted is efficient and effective in seeking an optimal solution in a complex, large-scale inverse problem. Furthermore, this study boosts our confidence in the continuum approach for approximating the site-scale saturated-unsaturated flow behavior in fractured rock formations, an important but tough issue that is frequently encountered in various subsurface applications.
Acknowledgments
The financial supports from the National Natural Science Foundation of China (No.
51925906)
and
the
National
Key
R&D
Program
of
China
(No.
2018YFC0407001) are gratefully acknowledged. The authors gratefully thank POWERCHINA HUADONG Engineering Corporation Limited for the field data supporting this study.
32
Appendix A. Supplementary data
Supplementary data associated with this article can be found, in the online version. References Asher, M.J., Croke, B.F., Jakeman, A.J., Peeters, L.J., 2015. A review of surrogate models and their application to groundwater modeling. Water Resour. Res. 51(8), 5957-5973. https://doi.org/10.1002/2015WR016967. Assouline, S., Tessier, D., Bruand, A., 1998. A conceptual model of the soil water retention
curve.
Water
Resour.
Res.
34
(2),
223–231.
https://doi.org/10.1029/97WR03039. Bandurraga, T., Bodvarsson, G., 1999. Calibrating hydrogeologic parameters for the 3-D site-scale unsaturated zone model of Yucca Mountain, Nevada. J. Contam. Hydrol. 38, 25–46. https://doi.org/10.1016/S0169-7722(99)00010-8. Barenblatt, G.I., Zheltov, I.P., Kochina, I.N., 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]. J APPL MATH MEC. 24(5), 1286-1303. https://doi.org/10.1016/0021-8928(60)90107-6. Bathe, K.J., Khoshgoftaar, M.R., 1979. Finite element free surface seepage analysis without
mesh
iteration.
Int
J
Numer
Anal
Met.
3(1),
13–22.
https://doi.org/10.1002/nag.1610030103. Bertels, S.P., Dicarlo, D.A., Blunt, M.J., 2001. Measurement of aperture distribution, capillary pressure, relative permeability, and in situ saturation in a rock
33
fracture using computed tomography scanning. Water Resour. Res. 37 (3), 649–662. https://doi.org/10.1029/2000WR900316. Borsi, I., Farina, A., Primicerio, M., 2006. A rain water infiltration model with unilateral boundary condition: qualitative analysis and numerical simulations. Math. Methods Appl. Sci. 29, 2047–2077. https://doi.org/10.1002/mma.763. Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media. Hydrology Paper No. 3. Colorado State University, Fort Collins, CO. Burdine, N.T., 1953. Relative permeability calculations from pore-size distribution data.
Trans.
Am.
Inst.
Min.
Metall.
Pet.
Eng.
198,
71–77.
https://doi.org/10.2118/225-G. Carrera, J., Alcolea, A., Medina, A., Hidalgo, J., Slooten, L.J., 2005. Inverse problem in
hydrogeology.
Hydrogeol
J.
13
(1),
206–222.
https://doi.org/10.1007/s10040-004-0404-7. Chen, C., Horne, R.N., 2006. Two-phase flow in rough-walled fractures: Experiments and
a
flow
structure
model.
Water
Resour.
Res.
42,
W03430.
https://doi.org/10.1029/2004WR003837. Chen, Y.F., Hong, J.M., Zheng, H.K., Li, Y., Hu, R., Zhou, C.B., 2016. Evaluation of groundwater leakage into a drainage tunnel in Jinping-I arch dam foundation in Southwestern China: A case study. Rock Mech. Rock Eng. 49 (3), 961–979. https://doi.org/10.1007/s00603-015-0786-y. Chen, Y.F., Hu, S.H., Hu, R., Zhou, C.B., 2015. Estimating hydraulic conductivity of fractured rocks from high pressure packer tests with an Izbash’s law-based
34
empirical
model.
Water
Resour.
Res.
51,
2096–2118.
https://doi.org/10.1002/2014WR016458. Chen, Y.F., Ling, X.M., Liu, M.M., Hu, R., Yang, Z.B., 2018. Statistical distribution of hydraulic conductivity of rocks in deep-incised valleys, southwest China. J. Hydrol. 566, 216-226. https://doi.org/10.1016/j.jhydrol.2018.09.016. Chen, Y.F., Zhou, C.B., Jing, L., 2009. Modeling coupled THM processes of geological porous media with multiphase flow: theory and validation against laboratory and field scale experiments. Comput. Geotech. 36 (8), 1308–1329. https://doi.org/10.1016/j.compgeo.2009.06.001. Chen, Y.F., Zhou, C.B., Zheng, H., 2008. A numerical solution to seepage problems with complex drainage systems. Comput Geotech. 35(3), 383–393. https://doi.org/10.1016/j.compgeo.2007.08.005. Desai, C.S., Li, G.C., 1983. A residual flow procedure and application for free surface flow
in
porous
media.
Adv
Water
Resour.
6(1),
27–35.
https://doi.org/10.1016/0309-1708(83)90076-3. Doughty, C., 1999. Investigation of conceptual and numerical approaches for evaluating moisture, gas, chemical, and heat transport in fractured unsaturated rock.
J
Contam
Hydrol.
38(1-3),
69-106.
https://doi.org/10.1016/S0169-7722(99)00012-1. Eaton, R.R., McCord, J.T., 1995. Effective hydraulic conductivities in unsaturated heterogeneous media by Monte Carlo simulation. Transport Porous Med. 18(3), 203-216. https://xs.scihub.ltd/https://doi.org/10.1007/BF00616931.
35
Evans, D.D., Rasmussen, T.C., 1991. Unsaturated flow and transport through fractured rock related to high-level waste repositories, Final Report: Phase III. NUREG/CR-5581, 67 pp. https://www.osti.gov/servlets/purl/137966. Faybishenko, B., and Finsterle, S., 2000. Tensiometry in fractured rocks, in Zhang, D., Winter, C.L., eds., Theory, Modeling, and Field Investigation in Hydeogeology: A Special Volume in Honor of Shlomo P. Neuman’s 60th Birthday: Boulder, Colorado, Geological Society of America Special Paper 348, p. 161-174. Finsterle, S., 2000. Using the continuum approach to model unsaturated flow in fractured
rock.
Water
Resour.
Res.
36
(8),
2055–2066.
https://doi.org/10.1029/2000WR900122. Flint, L.E., 2003. Physical and hydraulic properties of volcanic rocks from Yucca Mountain,
Nevada.
Water
Resour.
Res.
39
(5).
https://doi.org/10.1029/2001WR001010. Fredlund, D.G., Xing, A.Q., 1994. Equations for the soil-water characteristic curve. Can. Geotech. J. 31, 521–532. https://doi.org/10.1139/t94-061. Gallipoli, D., 2012. A hysteretic soil-water retention model accounting for cyclic variations of suction and void ratio. Géotechnique 62 (7), 605–616. http://dx.doi.org/10.1680/geot.11.P.007. Garcia, L.A., Shigidi, A., 2006. Using neural networks for parameter estimation in ground
water.
J.
Hydrol.
https://doi.org/10.1016/j.jhydrol.2005.05.028.
36
318,
215–231.
Gong, W.Y., Cai, Z.H., Jiang, L.X., 2008. Enhancing the performance of differential evolution using orthogonal design method. Appl. Math. Comput. 206, 56–69. https://doi.org/10.1016/j.amc.2008.08.053. Gribb, M.M., 1996. Parameter estimation for determining hydraulic properties of a fine sand from transient flow measurements. Water Resour. Res. 32(7), 1965– 1974. https://doi.org/10.1029/96WR00894. Gribb, M.M., Šimůnek, J., Leonard, M.F., 1998. Development of cone penetrometer method to determine soil hydraulic properties. J Geotech Geoenviron. 124(9), 820–829. https://doi.org/10.1061/(ASCE)1090-0241(1998)124:9(820). Guarracino, L., 2006. A fractal constitutive model for unsaturated flow in fractured hard
rocks.
J.
Hydrol.
324,
154–162.
https://doi.org/10.1016/j.jhydrol.2005.10.004. Guo, T., Zhang, S., Zou, Y., Xiao, B., 2015. Numerical simulation of hydraulic fracture propagation in shale gas reservoir. J Nat Gas Sci Eng. 26, 847-856. https://doi.org/10.1016/j.jngse.2015.07.024. Hvorslev, M.L., 1951. Time lag and soil permeability in groundwater observations. Bull 36. Waterways Experiment Station Corps of Engineers. U. S. Army, Vicksburg, Mississippi. Hughson, D.L., Yeh, T.-C.J., 2000. An inverse model for three-dimensional flow in variably saturated porous media. Water Resour. Res. 36 (4), 829–839. https://doi.org/10.1029/2000WR900001. Hong, J.M., Chen, Y.F., Liu, M.M., Zhou, C.B., 2017. Inverse modelling of
37
groundwater flow around a large-scale underground cavern system considering the excavation-induced hydraulic conductivity variation. Comput Geotech. 81, 346-359. https://doi.org/10.1016/j.compgeo.2016.09.008. Hsieh, P.A., Neuman, S.P., 1985. Field determination of the three-dimensional hydraulic conductivity tensor of anisotropic media. Water Resour. Res. 21 (11), 1655–1665. https://doi.org/10.1029/WR021i011p01667. Hu, R., Chen, Y.F., Liu, H.H., Zhou, C.B., 2013. A water retention curve and unsaturated hydraulic conductivity model for deformable soils: consideration of the change in pore-size distribution. Géotechnique 63 (16), 13891405. https://doi.org/10.1680/geot.12.P.182. Hu, R., Chen, Y.F., Liu, H.H., Zhou, C.B., 2015. A relative permeability model for deformable soils and its impact on coupled unsaturated flow and elastoplastic deformation processes. Sci. China Tech. Sci. 58 (11), 1971–1982. https://xs.scihub.ltd/https://doi.org/10.1007/s11431-015-5948-3. Hu, R., Chen, Y.F., Zhou, C.B., Liu, H.H., 2017. A numerical formulation with unified unilateral boundary condition for unsaturated flow problems in porous media.
Acta
Geotech.
12
(2),
277291.
https://xs.scihub.ltd/https://doi.org/10.1007/s11440-016-0475-3. Illman, W.A., Hughson, D.L., 2005. Stochastic simulations of steady state unsaturated flow in a three-layer, heterogeneous, dual continuum model of fractured rock. J. Hydrol. 307, 17–37. https://doi.org/10.1016/j.jhydrol.2004.09.015. Illman, W.A., Neuman, S.P., 2000. Type-curve interpretation of multi-rate single-hole
38
pneumatic injection tests in unsaturated fractured rock. Groundwater 38(6), 899–911. https://doi.org/10.1111/j.1745-6584.2000.tb00690.x. Illman, W.A., Neuman, S.P., 2001. Type-curve interpretation of a cross-hole pneumatic test in unsaturated fractured tuff. Water Resour. Res. 37(3), 583– 604. https://doi.org/10.1029/2000WR900273. Illman, W.A., Neuman, S.P., 2003. Steady-state analyses of cross-hole pneumatic injection tests in unsaturated fractured tuff. J. Hydrol. 281(1-2), 36–54. https://doi.org/10.1016/S0022-1694(03)00199-9. Ireson, A.M., Wheater, H.S., Butler, A.P., Mathias, S.A., Finch, J., Cooper, J.D., 2006. Hydrological processes in the Chalk unsaturated zone–insights from an intensive field monitoring programme. J. Hydrol. 330(1-2), 29–43. https://doi.org/10.1016/j.jhydrol.2006.04.021. Jiang, Q., Feng, X.T., Hatzor, Y.H., Hao, X.J., Li, S.J., 2014. Mechanical anisotropy of columnar jointed basalts: An example from the Baihetan hydropower station,
China.
Eng.
Geol.
175,
35–45.
https://doi.org/10.1016/j.enggeo.2014.03.019. Karpouzos, D.K., Delay, F., Katsifarakis, K.L., de Marsily, G., 2001. A multipopulation genetic algorithm to solve the inverse problem in hydrogeology.
Water
Resour.
Res.
37,
2291–2302.
https://doi.org/10.1029/2000WR900411. Kool, J.B., Parker, J.C., 1988. Analysis of the inverse problem for transient unsaturated
flow.
Water
Resour.
39
Res.
24
(6),
817–830.
https://doi.org/10.1029/WR024i006p00817. Kueper, B.H., McWhorter, D.B., 1992. The use of macroscopic percolation theory to construct large-scale capillary pressure curves. Water Resour. Res. 28 (9), 2425–2436. https://doi.org/10.1029/92WR01176. Li, X., Chen, Y.F., Hu, R., Yang, Z.B., 2017a. Towards an optimization design of seepage control: A case study in dam engineering. Sci China Technol Sc. 60(12), 1903–1916. https://doi.org/10.1007/s11431-016-9160-y. Li, Y., Chen, Y.F., Zhang, G.J., Liu, Y., Zhou, C.B., 2017b. A numerical procedure for modeling the seepage field of water-sealed underground oil and gas storage caverns.
Tunn.
Undergr.
Space
Tech.
66,
5663.
https://doi.org/10.1016/j.tust.2017.04.002. Li, Y., Chen, Y.F., Zhou, C.B., 2014. Hydraulic properties of partially saturated rock fractures subjected to mechanical loading. Eng. Geol. 179, 2431. https://doi.org/10.1016/j.enggeo.2014.06.019. Liu, H.H., Bodvarsson, G.S., 2001. Constitutive relations for unsaturated flow in fracture
network.
J.
Hydrol.
252,
116–125.
https://doi.org/10.1016/S0022-1694(01)00449-8. Liu, H.H., Bodvarsson, G.S., 2003. Upscaling of constitutive relations in unsaturated heterogeneous
tuff
matrix.
J.
Hydrol.
276,
198–209.
https://doi.org/10.1016/S0022-1694(03)00071-4. Liu, H.H., Doughty, C., Bodvarsson, G.S., 1998. An active fracture model for unsaturated flow and transport in fractured rocks. Water Resour. Res. 34 (10),
40
2633-2646. https://doi.org/10.1029/98WR02040. Liu, H.H., Haukwa, C.B., Ahlers, C.F., Bodvarsson, G.S., Flint, A.L., Guertal, W.B., 2003. Modeling flow and transport in unsaturated fractured rock: an evaluation of
the
continuum
approach.
J.
Contam.
Hydrol.
62,
173–188.
https://doi.org/10.1016/S0169-7722(02)00170-5. Lu, Z., Kwicklis, E.M., 2012. Numerical evaluation of effective unsaturated hydraulic properties of fractured rocks using a stochastic continuum approach. Vadose Zone J. 11 (4). https://doi.org/10.2136/vzj2011.0164. Mantoglou, A., Gelhar, L.W., 1987a. Stochastic modeling of large-scale transient unsaturated
flow
systems.
Water
Resour.
Res.
23(1),
37-46.
https://doi.org/10.1029/WR023i001p00037. Mantoglou, A., Gelhar, L., 1987b. Capillary tension head variance, mean soil moisture content, and effective specific soil moisture capacity of transient unsaturated flow in stratified soils. Water Resour. Res. 23(1), 47-56. https://doi.org/10.1029/WR023i001p00047. Mantoglou, A., Gelhar, L.W., 1987c. Effective hydraulic conductivities of transient unsaturated flow in stratified soils. Water Resour. Res. 23(1), 57-67. https://doi.org/10.1029/WR023i001p00057. Mao, D., Wan, L., Yeh, T.-C.J., Lee, C.-H., Hsu, K.-C., Wen, J.-C., Lu, W., 2011. A revisit of drawdown behavior during pumping in unconfined aquifers. Water Resour. Res. 47(5). https://doi.org/10.1029/2010WR009326. Mathias, S.A., Butler, A.P., Jackson, B.M., Wheater, H.S., 2006. Transient
41
simulations of flow and transport in the Chalk unsaturated zone. J. Hydrol. 330, 10–28. https://doi.org/10.1016/j.jhydrol.2006.04.010. Montazer, P., Wilson, W.E., 1984. Conceptual hydrological model of flow in the unsaturated zone, Yucca Mountain, Nevada. USGS Water-Resources Investigation Report, 84-4345, 55 pp. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous
media.
Water
Resour.
Res.
12
(3),
513–522.
https://doi.org/10.1029/WR012i003p00513. National Research Council, 1996. Rock Fractures and Fluid Flow-Contemporary Understanding and Applications. National Academy Press, Washington D.C. Nuth, M., Laloui, L., 2008. Advances in modelling hysteretic water retention curve in deformable
soils,
Comput.
Geotech.,
35
(6),
835–844.
https://doi.org/10.1016/j.compgeo.2008.08.001. Perkins, K., Johnson, B.D., Mirus, B.B., 2014. Measurement of unsaturated hydraulic properties and evaluation of property-transfer models for deep sedimentary interbeds, Idaho National Laboratory, Idaho. Report 2328–0328. Persoff, P., Pruess, K., 1995. Two-phase flow visualization and relative permeability measurement in natural rough-walled rock fractures. Water Resour. Res. 31 (5), 1175–1186. https://doi.org/10.1029/95WR00171. Peters, R.P., Klavetter, E.A., 1988. A continuum model for water movement in an unsaturated fractured rock mass. Water Resour. Res. 24 (3), 416–430. https://doi.org/10.1029/WR024i003p00416.
42
Pruess, K., Tsang, Y.W., 1990. On two-phase relative permeability and capillary pressure of rough-walled rock fractures. Water Resour. Res. 26, 1915–1926. https://doi.org/10.1029/WR026i009p01915. Rafiai, H., Jafari, A., Mahmoudi, A., 2013. Application of ANN-based failure criteria to rocks under polyaxial stress conditions. Int. J. Rock Mech. Min. Sci. 59, 42–49. https://doi.org/10.1016/j.ijrmms.2012.12.003. Reitsma, S., Kueper, B. H., 1994. Laboratory measurement of capillary pressure-saturation relationships in a rock fracture. Water Resour. Res. 30, 865–878. https://doi.org/10.1029/93WR03451. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics 1, 318–333. https://doi.org/10.1063/1.1745010. Rulon, J.J., Rodway, R., Freeze, R.A., 1985. The development of multiple seepage faces
on
layered
slopes.
Water
Resour.
Res.
21,
1625–1636.
https://doi.org/10.1029/WR021i011p01625. Salem, M.B., Tomaso, L., 2018. Automatic selection for general surrogate models. Struct
Multidiscip
O.
58(2),
719-734.
https://xs.scihub.ltd/https://doi.org/10.1007/s00158-018-1925-3. Song, X., Shi, L., Ye, M., Yang, J., Navon, I.M., 2014. Numerical comparison of iterative ensemble Kalman filters for unsaturated flow inverse modeling. Vadose Zone J. 13(2). https://doi.org/10.2136/vzj2013.05.0083. Tarantino, A., 2009. A water retention model for deformable soils. Géotechnique 59 (9), 751–762. https://doi.org/10.1680/geot.7.00118.
43
van Genuchten, M., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soil. Soil Sci. Soc. Am. J. 44 (5), 892–898. https://doi.org/10.2136/sssaj1980.03615995004400050002x. Wang, J.S.Y., Narasimhan, T.N., 1985. Hydrologic mechanisms governing fluid flow in a partially saturated, fractured porous medium. Water Resour. Res. 20 (12), 1861–1874. https://doi.org/10.1029/WR021i012p01861. Wang, X.S., Jiang, X.W., Wan, L., Song, G., Xia, Q., 2009. Evaluation of depth-dependent
porosity
and
bulk
modulus
of
a
shear
using
permeability-depth trends. Int. J. Rock Mech. Min. Sci. 46, 1175–1181. https://doi.org/10.1016/j.ijrmms.2009.02.002. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE. J. 3(03), 245-255. https://doi.org/10.2118/426-PA. Weerakone, W., Wong, R.C.K., Mehrotra, A.K., 2012. Measurement of capillary pressure curve of DNAPL in a water-saturated sandstone fracture. J. Geotech. Geoenviron.
Eng.
138,
614–624.
https://doi.org/10.1061/(ASCE)GT.1943-5606.0000590. Wu, Y.S., Liu, H.H., Bodvarsson, G.S., 2004. A triple-continuum approach for modeling flow and transport processes in fractured rock. J. Contam. Hydrol. 73, 145–179. https://doi.org/10.1016/j.jconhyd.2004.01.002. Wu, Y.S., Ritcey, A.C., Bodvarsson, G.S., 1999. A modeling study of perched water phenomena in the unsaturated zone at Yucca Mountain. J. Contam. Hydrol. 38, 157–184. https://doi.org/10.1016/S0169-7722(99)00015-7.
44
Wu, Y.S., Zhang, K., Liu, H.H., 2006. Estimating large-scale fracture permeability of unsaturated rock using barometric pressure data. Vadose Zone J. 5(4), 1129– 1142. https://doi.org/10.2136/vzj2006.0015. Yan, C., Jiao, Y.Y., 2018. A 2D fully coupled hydro-mechanical finite-discrete element model with real pore seepage for simulating the deformation and fracture of porous medium driven by fluid. Comput Struct. 196, 311-326. https://doi.org/10.1016/j.compstruc.2017.10.005. Yeh, T.C.J., Gelhar, L.W., Gutjahr, A.L., 1985a. Stochastic analysis of unsaturated flow in heterogeneous soils: 1. Statistically isotropic media. Water Resour. Res. 21(4), 447-456. https://doi.org/10.1029/WR021i004p00447 Yeh, T.C.J., Gelhar, L.W., Gutjahr, A.L., 1985b. Stochastic analysis of unsaturated flow in heterogeneous soils: 2. Statistically anisotropic media with variable α, Water
Resour.
Res.
21(4),
457-464.
https://doi.org/10.1029/WR021i004p00457. Yeh, T.C.J., Gelhar, L.W., Gutjahr, A.L., 1985c. Stochastic analysis of unsaturated flow in heterogeneous soils: 3. Observations and applications. Water Resour. Res. 21(4), 465-471. https://doi.org/10.1029/WR021i004p00465. Younes, A., Mara, T.A., Voltz, M., Guellouz, L., Baalousha, H.M., Fahs, M., 2018. A new efficient Bayesian parameter inference strategy: Application to flow and pesticide transport through unsaturated porous media. J. Hydro. 563, 887–899. https://doi.org/10.1016/j.jhydrol.2018.06.043. Zhang, M., Hu, L., Yao, L., Yin, W., 2017. Surrogate models for sub-region
45
groundwater management in the beijing plain, China. Water-Sui. 9(10), 766. https://doi.org/10.3390/w9100766. Zhang, L.M., Fredlund, D.G., 2003. Characteristics of water retention curves for unsaturated fractured rocks. Second Asian Conference on Unsaturated Soils (UNSAT-ASIA 2003), Osaka, Japan, p.425–428. Zheng, H., Liu, D.F., Lee, C.F., Tham, L.G., 2005. A new formulation of Signorini's type for seepage problems with free surfaces. Int J Numer Meth Eng. 64(1), 1– 16. https://doi.org/10.1002/nme.1345. Zhou, A., Sheng, D., Carter, J., 2012. Modelling the effect of initial density on soil-water
characteristic
curves.
Géotechnique
62
(8),
669–680.
https://doi.org/10.1680/geot.10.P.120. Zhou, C.B., Liu, W., Chen, Y.F., Hu, R., Wei, K., 2015. Inverse modeling of leakage through a rockfill dam foundation during its construction stage using transient flow model, neural network and genetic algorithm. Eng. Geol. 187, 183–195. https://doi.org/10.1016/j.enggeo.2015.01.008. Zhou, C.B., Zhao, X.J., Chen, Y.F., Liao, L., Liu, M.M., 2018. Interpretation of high pressure pack tests for design of impervious barriers under high-head conditions.
Eng.
Geol.
234,
112–121.
https://doi.org/10.1016/j.enggeo.2018.01.006. Zhou, H., Gómez-Hernández, J.J., Li, L., 2014. Inverse methods in hydrogeology: Evolution
and
recent
trends.
Adv.
Water
https://doi.org/10.1016/j.advwatres.2013.10.014.
46
Resour.
63,
22–37.
Zimmerman, R.W., Chen, G., Hadgu, T., Bodvarsson, G.S., 1993. A numerical dual-porosity model with semianalytical treatment of fracture/matrix flow. Water Resour. Res. 29(7), 2127-2137. https://doi.org/10.1029/93WR00749.
47
48
List of Figures (12 figures in total)
Fig.1. Location and landscape of the Baihetan dam site.
Precipitation and evaporation rate (mm)
200 Precipitation Evaporation
150 100 50 0 -50 -100 1
2
3
4
5
6
7 Month
8
9
10
11
12
Fig.2. Mean monthly precipitation measured at Baihetan weather station (635 m a.s.l.) and mean monthly evaporation from the land surface in the study area.
49
Fig.3. Geologic cross-section I-I at the dam site (with its trace shown in Fig. 1): (a) regional view and (b) enlarged view for the near-bank region on the left bank.
50
Fig.4. Groundwater flow system and locations of boreholes and exploratory adits.
Fig.5. Schematic diagram of groundwater level observation at boreholes penetrating a tuff zone. The drilling of boreholes across the tuff zone resulted in drawdown of groundwater level above it. The groundwater levels were measured after a concrete plug was constructed at the intersection between the borehole and the tuff zone, and the vertical distance between the upper and lower groundwater levels were then calculated. 51
Fig.6. 2D finite element mesh for the study area: (a) overall view, and (b) enlarged view for the near-bank region.
52
Fig.7. Sketch of the two hidden-layer BPNN model. The input layer contains 19 neurons for input of the unsaturated hydraulic parameters P and the time-dependent precipitation/evaporation boundary conditions. The output layer contains 12 neurons for output of the time series data of groundwater level in 9 boreholes and discharge at 3 weirs. The numbers of neurons at the two hidden layers (22 and 35) were determined by minimizing an error function on a test data set with a trial-and-error method.
Fig.8. The estimated water retention curves for rocks in different rock zones. 53
(a) 770
Measured water level Predicted water level Precipitation
760
Water level (m)
750
760 750
ZK1123
740
740 730
730 720
720
ZK1119
710
710
700
700
ZK9319
690 680 640
Precipitation (mm)
770
690 680 640
ZK929
630
630
620
620
610
610
150
150
100
100
50
50
0 2009/12/15 2010/02/15
2010/05/15
2010/08/15
Time
2010/11/15
0 2011/02/15
(b)
750
ZK512
Water level (m)
740
Measured water level Predicted water level Precipitation
730
ZK412
740 730
720
720
710
710
700
700
ZK312
690
Precipitation (mm)
750
690
680
680
670
670
150
150
100
100
50
50
0
0
2005/12/15 2006/04/15 2006/08/15 2006/12/15 2007/04/15 2007/08/15
Time
54
(c)
Precipitation (mm)
Water level (m)
748
Water level (m) Precipitation (mm)
748
744
744
ZK37 726
726
724
724
722
722
720
720
150
150
100
100
50
50
0
2001/12/15
(d)
Measured water level Predicted water level Precipitation
2002/03/15
2002/06/15
2002/09/15
Time
0
2002/12/15 2003/02/15
Measured water level Predicted water level Precipitation
780 774
780 774
768
768
762
762
ZK112
756
756
710 700 690 680
710 700 690 680
150
150
100
100
50
50
0
2005/09/15 2006/03/15
2006/11/15
2007/07/15
Time
2010/05/15
0
2010/11/15
Fig.9. Comparison of the observed (dot points) and calculated (solid lines) groundwater levels at observation boreholes in the slope. The rainfall intensity is also plotted in each figure.
55
Fig. 10. Water pressure head contours (a) in wet season (on July 15, 2000) and (b) in dry season (on January 15, 2001).
56
Fig.11. Water head contours on May 15, 2011 after the exploratory adits were excavated.
300
200
2.0
May, 2011
1.5 150 109
100
0.64
50 0
Excavation began May, 2002
2002
2004
1.0 0.5
Flow Rate (L/min)
250
Flow Rate (L/min)
2.5
Prediction Measurement PD61 PD61-3 PD61-4 PD61-5
0.036
2006 2008 Time (year)
2010
0.0 2012
Fig.12. Variation of discharge into the exploratory adits during the exploration period. Also plotted are the available measurements in May 2011.
57
58
List of Tables (5 tables in total)
Table 1 Specific storage and hydraulic conductivity for rocks at the site. Statistics of hydraulic conductivity (K) Standard Rock No. of Range (m/s) Mean (m/s) deviation tests (m/s) Loose sediments 4.31×106 1.41×1048.20×104 7.00104 2.40104 63 6 6 5 5 High-moderate PZ 2.94×10 5.69×10 1.94×10 1.00×10 3.48×106 8 Moderate PZ 1.91×106 1.15×1071.29×105 3.00106 1.18106 128 Moderate-low PZ 1.53×106 1.55×1089.32×106 9.00107 4.96107 561 6 9 6 7 7 Low PZ 1.53×10 9.06×10 4.14×10 3.0010 3.2210 1603 Very low PZ 9.42×107 6.47×1091.26×106 8.00108 9.39108 694 Fault F17 5.78×106 5.18×1071.61×105 2.75106 1.95107 11 Faults F129, F133, F134* 4.22×106 / 2.00106 / / * 6 8 6 7 7 Tuff zone C2 5.55×10 1.8110 1.1010 2.2010 4.4310 31 Tuff zone C3(A) * 5.55×106 1.811083.40106 3.01×107 4.07107 34 Tuff zone C3(B) * 5.55×106 3.101071.30×105 2.98106 / / Tuff zone C3-1(A) * 5.55×106 1.681082.70106 4.50×107 9.07107 33 Tuff zone C3-1(B) * 5.55×106 1.601067.20106 1.68106 / / 6 8 5 Shear zone LS331 5.55×10 2.3310 5.5010 8.00107 2.35106 23 Shear zone LS3318 5.55×106 3.621086.00106 8.15107 2.26107 8 Shear zone LS3319 5.55×106 7.121081.20106 3.73107 3.14107 18 * The tuff zones C and C 3 3-1 are divided into an unweathered subzone A and a weathered subzone B. The hydraulic conductivity of C2, C3(A) and C3-1(A) is highly anisotropic. The listed K values of C2, C3(A) and C3-1(A) represent the in-plane hydraulic conductivity (K//), and the values of K Specific storage Ss (m1)
perpendicular to the plane (K) are about two orders of magnitude smaller than the in-plane values (i.e., K=K///100). The hydraulic conductivity of all other rock zones is assumed to be isotropic. The hydraulic conductivity of the tuff zones C3(B) and C3-1(B) and the small-scale faults F129, F133 and F134 were estimated by hydraulic tests other than the packer tests.
Table 2 Comparison of the observed and calculated vertical distance between the upper and lower groundwater levels in wet season (on June 18, 2000). Borehole ZK925 ZK1119 ZK1121 ZK412
Distance to section I-I (m)
Tuff zone penetrated
0 111.6 111.6 259.6
C3 C3-1 C3-1 C3 59
Vertical distance of groundwater levels (m) Observed
Calculated
13.20 24.40 38.90 22.99
9.68 16.26 23.74 19.68
Table 3 Unsaturated hydraulic parameters for rocks at the site. Rock
θs
θr
Loose sediments High-moderate/moderate PZ Moderate-low/low PZ Very low PZ Tuff and shear zones Faults
0.32 0.20 0.12 0.08 0.25 0.30
0.05 0.03 0.02 0.01 0.04 0.05
α (m1) Range Estimated 0.110 0.600 0.0010.1 0.034 0.0010.1 0.012 0.0010.1 0.010 0.110 3.033 0.110 4.646
m Range Estimated 01 0.437 01 0.430 01 0.425 01 0.305 01 0.689 01 0.634
β Range Estimated 03 0.533 03 0.580 03 0.611 03 0.756 03 0.782 03 0.792
Table 4 Groundwater level observations in nine boreholes used for inverse modelling. Borehole
Distance to section I-I (m)
Collar elevation (m)
Depth (m)
Groundwater level (m)
Water level fluctuation (m)
Observation period
ZK929 ZK9319 ZK1123 ZK1119 ZK312 ZK412 ZK37 ZK112 ZK512
0 49.9 111.6 111.6 205.8 259.6 309.2 309.2 361.9
695.45 712.63 853.04 840.96 726.81 877.65 873.71 859.50 883.23
130.72 181.10 200.89 200.00 140.28 377.81 200.90 120.20 165.40
624.3620.6 687.0685.3 750.1745.9 717.3712.8 686.5677.5 715.0707.1 745.4741.9 767.7752.0 740.7733.2
3.7 1.7 4.2 4.5 9.0 7.9 3.5 15.7 7.5
2010/022011/01 2010/052011/01 2010/052011/01 2010/052011/01 2006/062007/07 2006/082007/10 2004/012005/01 2010/022011/01 2005/122007/10
Table 5 Comparison of the measured and calculated discharges in May 2011 at three weirs installed in the exploratory adit system. No.
Location of weir
Excavation time of adit
1 2 3 4 5 6
PD61 PD61-1 PD61-2 PD61-3 PD61-4 PD61-5
2002/052006/05 2006/032006/05 2006/052006/06 2006/072006/08 2006/092006/10 2009/012009/05
60
Discharge in May 2011 (L/min) Measured
Calculated
109 N/A N/A 0.036 N/A 0.64
113.08 0 0 0.041 0.259 0.841
61
evaporation
C3
46% 54%
C3-1
faults
65% 35%
C3
tuff zones
C3-1
multiple wa ter tables
river 50 m
Graphical Abstract: We investigate the saturated-unsaturated flow in a large-scale rock slope at the Baihetan dam site, Southwest China, using the continuum approach and inverse modeling. The unsaturated hydraulic parameters in water retention models are determined based on the time-series observations of groundwater level and the available discharge measurements. The typical features of the groundwater flow mainly governed by gently sloping tuff zones of low vertical permeability, such as multiple water tables, wedge-shaped unsaturated zones and significant groundwater level drawdown during ten-year site characterization, are well simulated.
62
63
Highlights
A methodology was presented for field-scale inverse modeling of unsaturated flow (82/85)
The determined unsaturated parameters are well representative of rocks at the site (84/85)
Numerical results reproduce the main features of multiple water tables at the site (84/85)
64
S0022-1694(20)30153-0 https://doi.org/10.1016/j.jhydrol.2020.124693 HYDROL 124693
To appear in:
Journal of Hydrology
Received Date: Revised Date: Accepted Date:
11 October 2019 19 January 2020 13 February 2020
Please cite this article as: Chen, Y-F., Yu, H., Ma, H-Z., Li, X., Hu, R., Yang, Z., Inverse modeling of saturatedunsaturated flow in site-scale fractured rocks using the continuum approach: A case study at Baihetan dam site, Southwest China, Journal of Hydrology (2020), doi: https://doi.org/10.1016/j.jhydrol.2020.124693
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2020 Published by Elsevier B.V.
Inverse modeling of saturated-unsaturated flow in site-scale fractured rocks using the continuum approach: A case study at Baihetan dam site, Southwest China
Yi-Feng Chen *, Hao Yu, Heng-Zhen Ma, Xing Li, Ran Hu *, Zhibing Yang State key Laboratory of Water Resources and Hydropower Engineering Science Wuhan University, Wuhan 430072, China Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering of the Ministry of Education, Wuhan 430072, China
Email: [email protected], [email protected]
Revised Manuscript January 20, 2020
A revised paper submitted for possible publication in Journal of Hydrology
*Corresponding
author. Tel.: +86-27-68774295; fax: +86-27-68774295. Email address: [email protected] (Y.-F. Chen), [email protected] (R. Hu). 1
Abstract
Modeling saturated-unsaturated flow in fractured rock formations remains a challenging issue due to the difficulties in properly calibrating the unsaturated flow properties for fractured rocks. On the basis of the continuum approach, this study uses inverse modeling to determine the unsaturated hydraulic parameters of fractured rocks, loose sediments, and geologic structures in a large-scale slope under difficult hydrogeologic conditions at the Baihetan dam site, Southwest China. The saturated hydraulic conductivities were determined a priori by packer tests. Field data, including time-series of groundwater level observations and available discharge measurements, is used. A combined procedure of orthogonal design, finite element forward analysis, artificial neural network, and genetic algorithm is adopted for the inverse analysis with high computational efficiency. The numerical results well reproduce the typical features of multiple water tables and wedge-shaped unsaturated zones controlled by gently-sloping tuff zones of low vertical permeability, showing 4665% of groundwater flowing subhorizontally along each tuff zone and the remaining 3554% infiltrating vertically and recharging deeper layers. The drastic drawdown of groundwater level after 10-year site characterization is also reproduced. The calculated discharge and water tables agree reasonably well with field observations, indicating that the estimated parameters are representative of the unsaturated hydraulic properties at the site. This work provides an efficient and effective methodology for seeking an optimal solution in large-scale inverse modeling of saturated-unsaturated flow, and further manifests the feasibility of the continuum 2
approach in predicting unsaturated flow in site-scale fractured rock formations. Keywords: Unsaturated flow; Hydraulic properties; Fractured rocks; Inverse modeling; Continuum approach
Nomenclature
Cw
Specific moisture capacity
D
Fractal dimension of flow paths
h*
Maximum allowable values of pressure head
h
Pressure head
h
Prescribed pressure head on the water head boundary h
hd
Minimum pressure head determined by temperature and air relative humidity on the evaporation boundary e
hp
Ponding depth of water on the infiltration boundary i
I(t)
Rainfall intensity as a function of time t
K
Saturated hydraulic conductivity tensor
kr
Relative permeability
M
Number of groundwater observation boreholes
N
Number of weirs
n
Outward unit normal vector to the boundary
Qj
Measured discharge data at weir j
Qj
Simulated discharge data at weir j
qn*
Maximum allowable flux on the third-type boundary
3
qn
Prescribed flux (positive for inflow) on the flux boundary q
Re
Maximum evaporation rate on the evaporation boundary e
Se
Effective saturation
Ss
Specific storage
s
Suction
t
Time
w
Weight coefficient
z
Positive upward vertical coordinate
, m, n
Empirical parameters in the van Genuchten model
r
Compressibility of rock
w
Compressibility of water
Volumetric water content
r
Residual volumetric water content
s
Saturated volumetric water content
e
Evaporation boundary
i
Rain infiltration boundary
s
Seepage face boundary
Sign indicator
Porosity
A step function of h
v
Flow velocity vector
w
Density of water
4
Positive parameter in the modified Mualem model
i
Measured groundwater level data in borehole i
i
Simulated groundwater level data in borehole i
1. Introduction
Modeling the saturated-unsaturated flow in fractured rocks has received considerable attention in the past decades in a wide range of applications such as nuclear waste disposal, contaminant transport, oil and gas extraction, underground oil storage, and landslide mitigation (Rulon et al., 1985; Liu et al., 1998; Li et al., 2017b). Accurately characterizing this flow process is extremely difficult and has been one of the most challenging issues in subsurface hydrology mainly because of the multi-scale geometries and complex interactions of the fracture-matrix system, the technical difficulties in measuring fracture flow in the field, and the multi-dimensional and heterogeneous flow behaviors (Evans and Rasmussen, 1991; Liu et al., 1998). In particular, the groundwater typically moves preferentially through fractures and faults of large scales (Peters and Klavetter, 1988; Illman and Hughson, 2005). Besides, multiple water tables and wedge-shaped unsaturated zones may develop if the stratigraphy consists of subhorizontal layers of contrasting permeabilities (Rulon et al., 1985). Various modeling approaches have been proposed to describe the multi-scale, multi-dimensional flow behaviors through the fracture-matrix system, including the equivalent continuum (Peters and Klavetter, 1988; Finsterle, 2000; Liu et al., 2003;
5
Lu and Kwicklis, 2012), dual-continuum (Illman and Hughson, 2005; Mathias et al., 2006), triple-continuum (Wu et al., 2004), discrete fracture network (Evans and Rasmussen, 1991), percolation network (Kueper and McWhorter, 1992) and stochastic representation approaches (Illman and Hughson, 2005; Lu and Kwicklis, 2012). The equivalent continuum model treats the fractured rock as a porous medium with an equivalent permeability obtained through homogenization of the flow process. This approach has been widely used for its simplicity and high computational efficiency (Bathe et al., 1979; Desai et al., 1983; Zheng et al., 2005; Chen et al., 2008, 2016; Li et al., 2017a). The discrete fracture network model can explicitly simulate flow in all fractures in a rock mass (Evans and Rasmussen, 1991). However, due to the large number of fractures involved, it is often computationally prohibitive to adopt this approach in site-scale applications. The dual-medium model simultaneously considers fast flow in the fractures and slow flow in the porous matrix (Barenblatt et al., 1960; Warren and Root, 1963; Zimmerman et al., 1993; Illman and Hughson, 2005; Mathias et al., 2006), but its application still has great limitations due to the uncertainty of the spatial distribution of the fracture system. The tripe-continuum model introduces an additional continuum (based on the dual continuum model) that requires one more parameter set for small fractures (Wu et al., 2004), which makes it even more difficult to use than the dual-continuum model. The stochastic approach assumes that the local hydraulic properties are realizations of spatially correlated random fields and derives the partial differential equations representing the large-scale flow conditions through averaging the local governing equations over the
6
ensemble of realizations (Yeh et al., 1985a, 1985b, 1985c; Mantoglou and Gelhar, 1987a, 1987b, 1987c). Among these approaches, the continuum model has dominated as the most commonly-used one for its efficiency in addressing large-scale problems because the presence of numerous fractures of different scales causes extreme difficulty in fracture network construction and flow modeling. The performance of a continuum model, however, largely depends on whether the unsaturated hydraulic properties (i.e., water retention curve and unsaturated hydraulic conductivity) of fractured rocks are properly determined and representative of the site conditions. Field measurements of permeability and porosity for unsaturated fractured rocks prove to be feasible by single-hole (e.g., Illman and Neuman, 2000) or cross-hole (Illman and Neuman, 2001, 2003) pneumatic injection tests without modifying the saturation as done by water infiltration tests. This is not the case, however, for direct measurement of unsaturated hydraulic properties representative of field scale behavior due to the complexity and variability of fracture geometries (National Research Council, 1996; Hughson and Yeh, 2000; Zhang and Fredlund, 2003; Wu et al., 2006), even though these measurements could be performed for core-scale samples of individual fracture (Wang and Narasimhan, 1985; Pruess and Tsang, 1990; Reitsma and Kueper, 1994; Persoff and Pruess, 1995; Bertels et al., 2001; Chen and Horne, 2006; Weerakone et al., 2012; Li et al., 2014) and rock matrix (Montazer and Wilson, 1984; Peters and Klavetter, 1988; Faybishenko and Finsterle, 2000; Flint, 2003). Therefore, the classic water retention curves (Brooks and Corey, 1964; van Genuchten, 1980; Fredlund and Xing, 1994) and relative permeability models
7
(Burdine, 1953; Mualem, 1976) developed for unsaturated soils and porous media are often borrowed for describing the unsaturated flow in fractured rocks. The development or extension of the water retention relations for unsaturated fractured rocks is based on the dual- or triple-porosity concept (Peters and Klavetter, 1988; Zhang and Fredlund, 2003; Wu et al., 2004), active fracture concept (Liu et al., 1998), fractal fracture geometry (Guarracino, 2006), and two-dimensional fracture network simulation and upscaling technique (Liu and Bodvarsson, 2001, 2003; Lu and Kwicklis, 2012). Inverse modeling provides a powerful tool for estimation of hydraulic parameters in water retention models from laboratory experiments, field groundwater observations and hydraulic tests (e.g., Kool and Parker, 1988; Liu et al., 1998, 2003; Bandurraga and Bodvarsson, 1999; Wu et al., 2004, 2006; Song et al., 2014; Younes et al., 2018). However, difficulties always arise in the inverse modeling due to the problems of ill-posedness, computational burden, and scales (Zhou et al., 2014; Hughson and Yeh, 2000). The objective of this study is to demonstrate the feasibility of modeling saturated-unsaturated flow in a large-scale rock slope of 5 km long located at the Baihetan dam site, Southwest China, based on the continuum approach and inverse modeling. This slope is characteristic of difficult hydrogeologic conditions, with multiple water levels in natural state due to the presence of gently-sloping, low-vertical permeability tuff zones and drastic depletion of water table after ten years of site characterization due to the excavation of exploratory adits and drilling of boreholes. The unsaturated hydraulic parameters in the water retention
8
relation and relative permeability model are estimated by inverse modeling, with the saturated hydraulic conductivities of rocks determined a priori by packer tests. The objective function is constructed by using the time-series observations of groundwater level and the available measurements of discharge. In order to achieve high computational efficiency, the combined procedure of orthogonal design (OD), finite element (FE) analysis, artificial neural network (ANN) and genetic algorithm (GA) presented by Zhou et al. (2015) is adopted to estimate the unsaturated hydraulic parameters of rocks at the site. The numerical results are compared with field observations to deepen our understanding of the saturated-unsaturated flow in the rock slope subjected to precipitation/evaporation cycles and exploratory adits excavation.
2. Site characterization
2.1. General description
The study site (Jiang et al., 2014) is situated on the left bank of the lower Jinsha River in Southwest China, where the world’s second largest hydropower project (i.e., Baihetan Project) consisting of a double curvature arch dam of 289 m high and two large-scale underground cavern systems is under construction (Fig. 1). The Jinsha River runs through the dam site from south to north (with a water level of about 591 m a.s.l. in dry seasons), forming a V-shaped valley by fast downcutting. The study area is characteristic of a mountainous landscape, with steep gradients of 5080 below the valley shoulder and gentler gradients towards the watershed between the Jinsha River and its tributary on the west side, the Yibu River (Fig. 1). The site lies in 9
the subtropical monsoon climate zone, where over 80% of the annual precipitation occurs during the wet season between May and October. Precipitation increases with elevation, with the measured mean annual rainfalls being 727.2 mm at 635 m a.s.l. near the dam site and 811.0 mm at 1,357 m, respectively. The mean monthly rainfall and evaporation at the study site are shown in Fig. 2. The average annual temperature is between 15 and 21. Surface runoff is infrequent and of short duration, with only two small perennial streams (Nos. 1 and 3 in Fig. 1) in the study area.
2.2. Geological settings
As shown in Fig. 3, the bedrocks outcropping in the study area are basaltic formations originated from multiple magmatic and volcanic eruption episodes. The formations belong to the Emei Mountain group of the upper Permian system (P2), which can be categorized into 4 basalt flow layers (P21P24) of 445470, 290340, 210260 and 85115 m thick, respectively, according to the eruption sequences. The rock layers mainly contain tholeiitic basalt, cryptocrystalline basalt, microcrystalline basalt and amygdaloidal basalt, typically with a layer of tuff developed on the top of each rock unit (Wang et al., 2009). The tuff layers are 0.1–1.75 m thick, and have regional continuity. The tuff is nonwelded, and is composed of vitric fragments and glass shards (40–50%), volcanic ash (30–40%), basalt clasts (5–10%) and plagioclase clasts (<5%) (Wang et al., 2009). Each tuff layer contains a shear zone (hereafter referred to as tuff zone) of 5–40 cm thick. The formations dip gently towards the right bank, with an orientation of N3050E/SE1525. The basalt rocks are overlain
10
by loose sediments (Q3 and Q4) of several to tens of meters thick on gentle hillslopes, and underlain by thick carbonate rocks (P1) deposited during the early Permian period (Fig. 3a). The main structures in the study area consist of subvertical faults (e.g., F17, F129, F133 and F134) and gently-sloping tuff zones (e.g., C2, C3-1 and C3) and shear zones (e.g., LS331, LS3318 and LS3319), as shown in Fig. 3b. Among them, F17 is a large-scale fault of 0.53 m thick, extending over 1400 m and mainly composed of well-cemented tectonic breccia and cataclasite of moderate permeability. Other faults (e.g., F129, F133 and F134) are of smaller sizes (about 0.11 m thick and extending up to hundreds of meters), and bounded by the tuff zones; these faults are also significant to groundwater infiltration as they outcrop and dip subvertically. The tuff zones are developed in the tuff layers, in which fault gouge layers of 1–5 cm thick are developed as a result of tectonic shear between basalt layers (Fig. S1 in Supplementary material), which makes the permeability of the tuff zones highly anisotropic, with a potential permeability contrast of 2 orders of magnitude between the directions parallel and perpendicular to the tuff layer. Except LS331, the shear zones within each basaltic formation are mostly of small scale, and mainly composed of basalt fragments. Besides the geologic structures mentioned above, four groups of critically-oriented fractures are developed in the basalt rocks, with their geometrical properties shown in Fig. 3b and Table S1 in Supplementary material.
11
2.3. Permeability of rocks
A total of 4413 conventional packer tests at 173 boreholes and 669 high pressure packer tests (Chen et al., 2015; Zhou et al., 2018) at 26 boreholes, together with a large number of pumping tests, recovery tests, slug tests, and in-situ seepage tests, were performed at the site to investigate the permeability of the fractured basalt rocks, structures (i.e., faults, tuff zones, and shear zones) and loose sediments. The boreholes were drilled vertically from the ground surface or from the floors of exploratory adits horizontally excavated at different elevations into the slope (Fig. 3b and Fig. 4). The diameter of the boreholes was 76 mm, and the length of the test sections varied between 3.7 and 6.4 m (Chen et al., 2018). The packer test data were interpreted with the Hvorslev (1951) equation. The test results indicate that the permeability of the fractured basalt rocks generally decreases with increasing depth (Fig. S2 in Supplementary material), and can be more correlated to the degree of weathering, the density of tensional joints and the level of in-situ stress, rather than the lithology. Downward from the ground surface, the hydraulic conductivity K can be classified into six categories or permeability zones (PZs, Fig. 3b): high (K 105), high-moderate (105 > K 3106), moderate (3106 > K 106), moderate-low (106 > K 3107), low (3107 > K 107), and very low (K < 107), all in units of m/s. The main structures are mostly of moderate-low to moderate permeability, but the permeability of the unweathered tuff zones and shear zones is highly anisotropic. Table 1 lists the statistics of the hydraulic conductivity of the rocks and sediments at the site. 12
2.4. Groundwater observations
Groundwater in the study area is principally recharged by precipitation, and discharged to streams and the Jinsha River. The groundwater flow and water table configuration in the basaltic formations are primarily controlled by the fracture system and the backbone structures consisting of subvertical faults and gently-sloping tuff zones and shear zones, as shown in Fig. 4. Due to the very low vertical permeability of the tuff zones and shear zones (Table 1), multiple water tables develop in the slope (Fig. 4), a phenomenon that has been frequently observed and well understood in layered hillsides (e.g., Rulon et al., 1985). Observations of groundwater in the slope started from 2001 through boreholes and exploratory adits, and lasted for over ten years. The exploratory adits were excavated horizontally into the slope, from which a few branches parallel to the river were also excavated (Fig. 4). The boreholes were drilled downward from the ground surface or from the floors of the adits. The mean groundwater levels in wet and dry seasons are shown in Fig. 4. The groundwater level in the slope, although influenced locally by fractures, generally increased with ground surface elevation. The groundwater level varied in a larger range (typically 1020 m) between wet and dry seasons in boreholes near the river bank (consistent with the change in river water level) than far away from the river bank. The presence of tuff zones and shear zones significantly influences the distribution and movement of groundwater in the slope. When an exploratory adit was excavated across the tuff zones or shear zones, a significant amount of discharge was 13
generally observed from the upper block of the tuff/shear zones, indicating that these structures are of low vertical permeability and hinder the vertical infiltration. The discharge into each adit decayed due to the advance of neighboring adits, but the total discharge into the whole adit system maintained a stable magnitude and was less influenced by precipitation events. For example, the total discharge during site characterization
from
the
exploratory
adit
PD61
and
its
five
branches
(PD61-1PD61-5), which cut through a series of faults and tuff/shear zones (e.g., F17 and C3), had rather stable measurements around 110 L/min by weirs. Furthermore, when the tuff zones (C3-1 and C3) were penetrated by the boreholes, a sharp drawdown of groundwater level generally occurred. Multiple groundwater levels were observed after the boreholes were isolated with concrete plugs at the intersection with the tuff zones. Table 2 lists the vertical distance between the groundwater levels observed above and below the tuff zones in four typical boreholes (Fig. 5). For example, in borehole ZK1119 a vertical distance of 24.4 m between the upper and lower groundwater levels was observed; similarly in borehole ZK412, a vertical distance of around 23.0 m between the groundwater levels was observed during the wet season. It should be noted, however, that as commonly done in dam engineering, the unsaturated zone was not monitored, even though techniques (e.g., tensiometers and neutron probes) are available for the measurements of flow quantities such as pressure head and water content (e.g., Ireson et al., 2006). This is a major limitation of this study, but motivates us to examine if the unsaturated zone parameters could be
14
reasonably estimated with the observations in the saturated zone.
3. Inverse modeling approach
3.1. Saturated-unsaturated flow model
In the framework of continuum approach (Peters and Klavetter, 1988) where the fractured rocks are treated as equivalent continua, the saturated-unsaturated water flow is governed by the following mass balance equation (e.g., Mao et al., 2011), also known as the Richards (1931) equation:
Cw + Ss
h v 0 t
(1)
where h is the pressure head, Cwh is the specific moisture capacity, is the volumetric water content, Ss is the specific storage, is a step function of h, = 0 for h < 0 (in unsaturated zone) and = 1 for h 0 (in saturated zone), t is time, and v is the flow velocity described by Darcy’s law:
v kr K (h z )
(2)
where K is the saturated hydraulic conductivity tensor, kr is the relative permeability depending on h or , z is the positive upward vertical coordinate, and h+z is the hydraulic head. The water mass balance equation Eq. (1) is subjected to the following initial condition in the domain of interest:
h(t = 0) h0 (in )
(3)
where h0 is the initial distribution of pressure head in , and the following boundary conditions: 15
(1) The water head boundary condition (Dirichlet boundary condition):
h(t ) h (t ) (on h )
(4a)
where h is the prescribed pressure head on the water head boundary h. (2) The flux boundary condition (Neumann boundary condition): qn (t ) v (t ) n qn (t ) (on q )
(4b)
where qn is the prescribed flux (positive for inflow) on the flux boundary q, and n is the outward unit normal vector to the boundary. (3) The unified lateral conditions on the third-type boundaries, including the seepage faces, rain infiltration surfaces and evaporation surfaces (Borsi et al., 2006; Hu et al., 2017): (h h* ) 0; (qn qn* ) 0 (on s i e ) * * (h h )(qn qn ) 0
(4c)
where s, i and e denote the seepage face boundary, rain infiltration surface, and the evaporation surface, respectively; is a sign indicator with = 1 on s and i, and = 1 on e; h* and qn* are the maximum allowable values of pressure head and flux, respectively, with h*0 and qn*0 on s, h*hp and qn*I(t) on i, and h*hd and qn*Re on e; hp is the ponding depth of water on the infiltration boundary i, and I(t) is the rainfall intensity; hd is the minimum pressure head determined by temperature and air relative humidity on the evaporation boundary e, and Re is a negative flux rate representing the maximum evaporation rate in the field. Eq. (4c) states that on the third-type boundaries (s, i and e), the boundary condition is either a prescribed pressure (h=h*) when the flux is no more than the allowable flux rate (qn<qn*), or a prescribed flux rate (qn=qn*) when the pressure is 16
no more than the allowable pressure (h<h*). The presence of the lateral boundary conditions in Eq. (4c) introduces strong nonlinearity of the partial differential equation (PDE) in Eq. (1) for saturated-unsaturated flow problems, which can be rigorously and effectively solved with a numerical procedure based on the parabolic variational inequality (PVI) method (Hu et al., 2017). The discretized PVI algorithm was detailed in Hu et al. (2017), and was implemented in the FE code THYME (Chen et al., 2009). As clarified by site characterization (Fig. 3 and Fig. 4), the groundwater at the site primarily flows along the direction perpendicular to the river (i.e., along the trace of cross-section I-I in Fig. 1), and the flow along the river direction is rather slow. A 2D FE mesh along section I-I was therefore created for saturated-unsaturated flow modeling, which contains 24,782 four-node isoparametric elements and 25,197 nodes, as plotted in Fig. 6. The size of the model domain is 5000 m long and 2190 m high. The topographic and geologic settings, such as the strata, faults, tuff/shear zones and PZs at cross-section I-I (Fig. 1), were well represented. It should be noted that the exploratory adits parallel to the profile (e.g., PD-61, Fig. 4) could not be accurately represented in the 2D mesh, and were modelled with equivalent spacing per unit length (Fig. 6b and Fig. S3 in Supplementary material). The boundary conditions were specified as follows: (1) the lateral boundaries on the mountain side and on the river side were both assumed to be groundwater divides, and thus were prescribed with a no-flow boundary; (2) the riverbed surface was imposed with a water head boundary according to the fluctuation of the river water
17
level; (3) the surfaces of the exploratory adits were taken as a seepage face boundary; (4) the slope surface was prescribed with a precipitation/evaporation boundary with hp=0 m and hd=200 m (corresponding to a saturation of 1618% for loose sediments and tuff zones). According to the site condition, precipitation was assumed to occur in the middle ten days of a month, and the mean monthly precipitation and evaporation rates were taken from Fig. 2; (5) the base of the model was assumed to be impermeable. The initial distribution of pressure head in the saturated zone of the slope (h00) was determined by estimating a water table that best fits the average groundwater level observations in boreholes with a steady-state flow model. In the unsaturated zone (h0<0), h0 was roughly estimated by assuming a vertically linear reduction of saturation to 60% at the ground surface. To eliminate the influence of h0 distribution on groundwater flow, a natural precipitation/evaporation process was first simulated until the groundwater flow in each season became stable. This process typically lasted for 360 months (from May 1971 to May 2001). The site exploration process was then modelled (see Table 5 for the excavation process of the exploratory adit system), lasting for 120 months (from May 2001 to May 2011). The initial time step and maximum time step in each simulation were taken as 7 hours and 7 days, respectively.
3.2. Constitutive relations and hydraulic parameters
The hydraulic properties that need to be determined for numerical modeling of water movement in saturated-unsaturated fractured rocks include the saturated
18
hydraulic conductivity K, the specific storage Ss, the water content and the relative permeability kr. The hydraulic conductivity tensor K (or its isotropic counterpart K) of fractured rocks at the field scale can be reasonably determined without technical difficulty by performing a sufficient number of field hydraulic tests (e.g., Hsieh and Neuman, 1985). Owing to the experiences that have been long accumulated, the determined hydraulic conductivity is generally of high confidence and representative of the field conditions (Table 1). The specific storage of rocks is estimated by Ss =wg(r+w), where r is the compressibility of rock, w is the compressibility of water, is the porosity, and w is the density of water (Table 1 and Table S2 in Supplementary material). The water content relation, also called the water retention curve, is defined as the relationship between the volumetric water content (or the degree of saturation Se) and suction s (i.e., the negative pressure head, s=h). Quite a number of water retention curves have been developed for unsaturated soils and porous media (Brooks and Corey, 1964; van Genuchten, 1980; Fredlund and Xing, 1994), and extended to incorporate the effects of hydraulic hysteresis and soil deformation (e.g., Assouline et al., 1998; Nuth and Laloui, 2008; Tarantino, 2009; Gallipoli, 2012; Zhou et al., 2012; Hu et al., 2013). These models are often borrowed for unsaturated fractured rocks, due to the technical difficulties in the field-scale measurements of water content relation in rocks of much higher heterogeneity. Numerical simulations and validations (Peters and Klavetter, 1988; Liu et al., 1998, 2003; Liu and Bodvarsson, 2001, 2003; Wu et al., 2004) have shown that the classic van Genuchten (1980) model, with a proper
19
parameterization representative of field conditions, is suitable for representing the water content relation in unsaturated fractured rocks, which reads: Se
r n m 1 | h | s r
(5)
where Se is the effective saturation, s and r are the saturated and residual volumetric water contents, respectively, and , m and n are empirical parameters with m=11/n. The relative permeability kr as a function of the degree of saturation Se is another constitutive relation of great importance in modeling the unsaturated flow through porous and fractured media. This relation can generally be obtained by using the water content relation in the classic conceptual models of Burdine (1953) and Mualem (1976). A major concern in the development of the relative permeability function is the tortuosity of the flow paths, which is often correlated to the fractal dimension of flow geometries (e.g., Guarracino, 2006) or the fraction of active fractures (e.g., Liu et al., 1998) to improve the prediction accuracy. Hu et al. (2015) proposed a modified Mualem relation for deformable soils by taking into account the change in pore size distribution and the tortuosity of flow paths, given by m kr 1 1 Se1/ m
2
(6)
where is a positive parameter related to the fractal dimension D of flow paths by
=2D2. A verification of the relation with experimental data (Hu et al., 2015) yields a global optimal estimate of =0.5 or D=1.25 with high confidence for five types of soils of quite different structures (including loess, silty loam, sandy loam and coarse sand). This suggests that similar to other modified van Genuchten or Brooks-Corey
20
permeability models (e.g., Peters and Klavetter, 1988; Liu et al., 1998; Liu and Bodvarsson, 2001; Guarracino, 2006), Eq. (6) may apply to fractured rocks if the parameter is properly determined. The constitutive relations in Eqs. (5) and (6) contain five parameters (i.e., m, ,
, s and r) for each rock unit, among which the saturated and residual volumetric water contents, s and r, can be empirically estimated with the porosity (or the volume fraction of connected fractures and pores) and the void structure of rocks. It has been understood that the seepage field in fractured media is generally less sensitive to the residual saturation r (Liu et al., 2003), and in the most simplified case, one may assume s = and r = 0. The other three parameters (m, and ) are to be representatively determined by inverse modeling, which will be stated in the next section. The values of s and r and the ranges of m, and for 6 rock zones are listed in Table 3.
3.3. The inverse modeling procedure As previously stated, the parameters K (or K), s and r in the unsaturated flow model can be reasonably measured or empirically estimated with engineering accuracy, but difficulty arises when the other parameters (m, and ) are to be measured for fractured rocks. This difficulty can be overcome with inverse modeling techniques by reproducing the dynamics of the seepage field and tuning the parameters to represent the site characteristics. 3.3.1. Objective function
21
Suppose r is the number of materials of contrasting unsaturated hydraulic properties in the domain . A vector P is defined to denote the hydraulic parameters of the media yet to be determined, with P = [m1, 1, 1, m2, 2, 2, ..., mr, r, r]T. In this study, r is chosen to be 6, representing the following six types of media: loose sediments, high-moderate/moderate permeability zones (PZs), moderate-low/low PZs, very low PZ, faults and tuff/shear zones, as listed in Table 3. It is noted here that some PZs and geologic structures are grouped to reduce the number of unknowns and hence the uncertainties and difficulties in the inverse modeling. To estimate the parameters P, an objective function is needed, which ensures a best fit between the field measurements and the numerical results (Carrera et al., 2005; Chen et al., 2016), given by
M (P) i i min f = 2 i 1 i 2
1
N Q (P) Q 2 j j 2 w 2 j 1 Qj 2 2
1
2 2 2
(7)
where M is the number of groundwater observation boreholes, i and i are the time series measurements of groundwater level and the corresponding numerical results in borehole i, respectively, N is the number of weirs, Q j and Qj are the time series measurements of discharge and the numerical results at weir j, respectively, w is a weight coefficient, and ||||2 denote the Euclidean norm of a vector. Eq. (7) utilizes all available time series measurements of both groundwater level and discharge to improve the inverse solution that is inevitably plagued by the problem of non-uniqueness. In this study, the available time series observations of groundwater level in 9 boreholes (Table 4) drilled in the near-bank region of the slope
22
(Fig. 4) and the available discharge measurements in May 2011 at 3 weirs (Table 5) installed in the exploratory adits (Fig. 4) were chosen for construction of the objective function. Note that except ZK929, the other boreholes were located in the neighbouring profiles, with their distances to cross-section I-I (Fig. 3 and 4) being listed in Table 4. With the increase of distance, the data may be less representative of the selected profile, which is the other weakness of this study, as discussed in Section 5. The weight coefficient w is introduced to ensure a balance between the relative errors of the groundwater level and discharge measurements, and a value of w=0.1 was adopted in this study by trials (Table S3 in Supplementary material). 3.3.2. Inverse modeling procedure To efficiently and reliably address the large-scale inverse problem in complex hydrogeologic conditions, this study adopted a combined procedure of orthogonal design (OD), finite element (FE) analysis, artificial neural network (ANN) and genetic algorithm (GA) recently presented in Zhou et al. (2015). The algorithms are given in Figs. S4-S6 in Supplementary material. OD is an effective statistics method to arrange multifactorial experiments (Gong et al., 2008). In this inverse modeling procedure, OD was employed to construct a small set of representative values for the hydraulic parameters P from a pre-defined parameter space. In this study, a total of 18 parameters (i.e., r=6) were to be estimated, with 7 levels specified for each parameter from the ranges empirically given in Table 3. Consequently, an orthogonal table L49(718) was constructed, which produces a set of 49 parameter combinations. These parameter combinations were then used as input parameters for forward modeling of
23
saturated-unsaturated flow (Hu et al., 2017; Chen et al., 2009), yielding the numerical time-series data of groundwater level and discharge at the corresponding observation points and times for each parameter combination. By doing so, only a limited number of time-consuming FE simulations were performed and hence the computational cost was greatly reduced. As a multivariate non-linear mapping and predictive tool that has been successfully applied in hydrology problems (e.g., Garcia and Shigidi, 2006), a back propagation neural network (BPNN) model was then constructed and trained (Fig. S5 in Supplementary material), using the parameter combinations and the prescribed time-dependent precipitation/evaporation boundary conditions (input data) as well as the numerical results (output data). Consequently, an implicit mapping from the parameter space to the unsaturated seepage fields at the observation points was obtained (Zhou et al., 2015). The BPNN model contains four layers, with 19 neurons in the input layer, 22 and 35 neurons in the two hidden layers and 13 neurons in the output layer, as shown in Fig. 7. As usual, the Sigmoid function was selected as the transfer function of the BPNN model. The Levenberg-Marquardt back-propagation algorithm combined with Bayesian regularization (Rafiai et al. 2013) was adopted for training of the model, in order to obtain quick training time and high generalization accuracy. Finally, as an appropriate optimization tool in hydrology problems (e.g., Karpouzos et al., 2001), GA was used to obtain a globally optimal estimate of P such that the objective function, Eq. (7), is minimized and a best fit between the field measurements and the numerical results calculated by the trained BPNN model is
24
ensured. This is achieved by randomly generating an initial population of hydraulic parameter combinations and then iteratively invoking the genetic operations of selection (reproduction), crossover and mutation (Fig. S6 in Supplementary material). In this study, the size of the initial population was set to 100, and the probabilities for crossover and mutation operations were specified to 0.9 and 0.1, respectively. The GA was run at least 10 times, and the best solution was selected as the optimal result.
4. Results
4.1. The results of inverse modeling Table 3 lists the estimated unsaturated hydraulic parameters for rocks in different zones, and the corresponding water retention curves are plotted in Fig. 8. The inverse modeling showed that the unsaturated flow behavior in the slope is most sensitive to the value of α, which controls the minimum threshold pressure at which air starts to penetrate fractured rocks. A larger α results in a higher infiltration rate into the slope and a stronger fluctuation of groundwater level during the precipitation-evaporation cycles. The estimated value of α for loose sediments (several to tens of meters thick) is 0.60 m1, which is comparable to the findings of Perkins et al. (2014), who showed that α is dependent on depth, varying from 0.63 m1 to 65.2 m1 for depth from 39.2 m to 84.3 m. The representative values of α for fractured rocks in different PZs ranges from 0.010 m1 to 0.034 m1, with lower value for rocks in lower PZ (Table 3). This result is comparable to the calibrated values (α=0.0011.63 m1) for rock matrix continuum at Yucca Mountain (Bandurraga and Bodvarsson, 1999; Flint, 2003). For faults and tuff/shear zones, the estimated values of α are 3.03 m1 and 4.65 m1,
25
respectively, again comparable to the calibrated results (α=5.2658.3 m1) for fractures at Yucca Mountain (Bandurraga and Bodvarsson, 1999; Liu and Bodvarsson, 2001), where the α value was found to be highly sensitive to the average aperture of fractures. The parameter m determines the slope of the water retention curve. Here, it is found that increasing m tends to suppress the fluctuation of groundwater level. The estimated values of m for fractured rocks are in a narrow range between 0.31 and 0.43 (corresponding to n=1.451.75), comparable to the calibrated results (m=0.140.43) for rock samples collected from various depths in different hydrogeologic units at Yucca Mountain (Flint, 2003). The representative value of m is moderate for loose sediments (m=0.44), and higher for faults and tuff/shear zones (m=0.630.69). The unsaturated flow process is less sensitive to the parameter β, but a larger β may result in smaller infiltration rate and groundwater level fluctuation. The estimated values of β for different rock zones are in a narrow range between 0.53 and 0.79 (corresponding to D=1.271.40), slightly larger but rather close to the optimal estimate for soils of quite different structures (Hu et al., 2015). This indicates that, as a first approximation, a value of β=0.5 similarly applies to fractured rocks.
4.2. Comparison of field measurements With the estimated parameters given in Table 3 and Fig. 8 as well as the numerical setup given in Section 3.1, simulations of groundwater flow in the slope were performed for 30 years of natural precipitation/evaporation cycles and 10 years of site exploration. Fig. 9 shows the fluctuations of the observed and calculated groundwater levels over tuff zone C3 (Fig. 4) in nine boreholes listed in Table 4. A comparison between the curves shows that the numerical simulation roughly 26
reproduces the groundwater level observations (Fig. 9a and Fig. 9b), except at boreholes ZK37 and ZK112 (Fig. 9c and Fig. 9d) located over 300 m away from the cross-section (Table 4). This difference is attributed to local terrain changes (Fig. 1). If the groundwater level is transformed to groundwater depth (i.e., with the topography changes compensated), then the difference becomes smaller. The predicted groundwater level at ZK1123 (Fig. 9a) is consistent with the rainfall time-series, but seems to have an opposite trend with the measurements. The groundwater level over tuff zone C3 either remained rather stable (Fig. 9a), or exhibited seasonal fluctuations (Fig. 9bd) as influenced by the rainfall infiltration through faults (e.g., F17 and F129, Fig. 4). Even for the latter, the numerical result well captures the trend of groundwater level variation. Simulated water pressure head distributions (in long-term stable state) in wet and dry seasons are shown in Fig. 10. It is interesting to observe that the presence of tuff zones C3 and C3-1 of low vertical permeability caused the distribution of groundwater in multiple layers. This phenomenon was well reproduced in the numerical model. The development of faults (e.g., F17 and F129) facilitated the infiltration of rainfall into the slope, and resulted in the rise of groundwater level over the tuff zones in wet seasons. The infiltrated groundwater discharged to the river either by subhorizontal movement along the tuff zones (for their high in-plane permeability) or by vertical infiltration across the tuff zones. The numerical simulation shows that for the groundwater across the ceiling of C3, about 46% of the water flowed along the tuff layer and 54% further infiltrated downward to recharge the layers below; similarly for C3-1, 65% of the groundwater flowed along it and the rest 35% infiltrated vertically through (see inset in Fig. 4). Consequently, the water table below the tuff zones remained rather stable, except the near-bank region influenced by river water level
27
fluctuations. Table 2 lists the vertical distance between the groundwater levels above and below the tuff zone C3 or C3-1 in four typical boreholes in the wet season, showing that the numerical result is comparable to the field observations. The weathering of tuff zones (C3 and C3-1) makes their permeability isotropic (K=K//) in subzone B (Table 1 and Fig. 6b), which facilitated the subvertical discharge of groundwater in the shallow region of the slope, causing the tuff zones to mostly remain unsaturated and hence the unsaturated zone to be wedge-shaped. As shown in Fig. 11, the water head contours on May 15, 2011, two years after the exploratory adit PD61 and its five branches PD61-1PD61-5 were excavated (Table 5), indicated a marked alteration of the seepage field (compare Fig. 10 and Fig. 11). The evolution of the calculated discharges out of PD61 and its three branches PD61-3PD61-5 is shown in Fig. 12, where the discharges are compared with the available measurements on May 15, 2011 (Table 5). The excavation of the adits led to increased vertical discharge of groundwater across the tuff zones and hence significant drawdown of water table. The discharges fluctuated significantly during the excavation of the adits, but overall tended to decrease to some stable magnitudes (Fig. 12). In May 2011, the simulated discharge from the adit system was about 114 L/min, which is in excellent agreement with the field measurements (Table 5). At that time, field observation showed that the 340 m near-face section of adit PD61 and the branches PD61-1 and PD61-2 were in the unsaturated zone, while the surrounding rocks of the rest section of PD61 and the branches PD61-3PD61-5 remained saturated. The numerical result (see Fig. 11) well reproduces this observation.
5. Discussion This study has demonstrated the feasibility of modeling saturated-unsaturated 28
flow in fractured rock formations with the continuum approach. This relies on proper estimation of the unsaturated hydraulic parameters of rocks in the van Genuchten model and the modified Mualem relation. The estimates were obtained by inverse modeling of saturated-unsaturated flow using the measurements of groundwater levels and discharges collected in the saturated zone. A major concern is whether the saturated flow is sensitive enough to the hydraulic parameters in the unsaturated zone. An additional sensitivity analysis (Table S4 in Supplementary material) confirms the sensitivities of , m and in the unsaturated zone to the input data obtained from the saturated zone. The parameter is most sensitively influenced, followed by m and , respectively, in sensitivity, which is consistent with the findings of Gribb et al. (1996) and Liu et al. (2003); besides, the parameters are more sensitive in the rock zones where the water tables are located. This partially justifies the inverse modeling strategy adopted in this study. It should be noted, however, that the proposed approach contains quite a number of simplifications and assumptions. First, the flow was modeled with a continuum model, where the details of the fracture system, the flow in individual fractures, and the fracture-matrix interaction were not explicitly accounted for. However, given the difficulty in obtaining (and thus the lack of) sufficient fracture network information from site characterization, the continuum model is arguably acceptable for this large-scale problem. This is partially justified by the capability of the continuum model in capturing the changes in groundwater level induced by changes in rainfall intensity during a typical period (Table S5 and Fig. S7 in Supplementary material). In fact, the equivalent continuum approach has been successfully applied in a large number of studies (e.g., Peters and Klavetter, 1988; Doughty, 1999; Wu et al., 1999; Finsterle, 2000; Liu et al., 2003).
29
Second, only the unsaturated hydraulic parameters were estimated by inverse modeling, where the hydraulic conductivities of rocks were directly determined from the field hydraulic tests. We made this choice because saturated hydraulic conductivities estimated from hydraulic tests have been shown to be representative of field scale flow properties of the rock (Gribb et al., 1998; Zhou et al., 2015; Hong et al., 2017). It should be pointed out that the heterogeneity and scale effect of the parameters within each rock zone were not considered. Since the field test data do not support a more detailed delineation of the heterogeneity at the site, it was technically infeasible to account for heterogeneity within each zone in the simulations (Eaton and McCord, 1995). Third, the flow was assumed to be two-dimensional in the slope, where the flow component perpendicular to the cross-section, albeit slow, was neglected. Due to insufficient data in the selected profile, the use of measurements from the neighboring profiles further gave rise to uncertainties in the 2D modeling. The 2D assumption is justifiable, given the fact that the groundwater at the site (the left bank) primarily flows along the direction perpendicular to the river (i.e., along the trace of section I-I in Fig. 1). In addition, the chosen cross-section (I-I in Fig. 1) roughly represents the geological conditions (including the topography, lithology, and geological structures) at the site. Therefore, the 2D model could be considered representative of the site characteristics. The simplification of the model from 3D to 2D is common in seepage analysis due to the constraint of computational resource in inverse calculation of 3D model (Yan and Jiao, 2018; Guo et al., 2015). A 3D model will undoubtedly provide a more realistic and effective evaluation, which is a focus of our future work. Finally, the use of ANN in the inverse modeling procedure significantly reduced the computational cost (Rafiai et al., 2013; Zhou et al., 2015; Hong et al., 2017). As a
30
surrogate model, ANN has become very popular and efficient in engineering problems (Zhang et al., 2017). It is difficult to estimate the unsaturated parameters in large-scale groundwater flow problems without using surrogate models (Asher et al., 2015; Salem and Tomaso, 2018). One should, however, be cautious when using surrogate models in inverse analysis because the reduced computational time possibly comes at the cost of oversmoothing the relation between the input and output data, depending on the choice and size of the training samples. In this study, we used a combination of orthogonal design and forward modeling to create proper training samples for the ANN model. It is nevertheless acknowledged that further improvement is possible in the surrogate modeling. Despite of these limitations, the proposed approach effectively reproduced the multiple water tables and wedge-shaped unsaturated zones in a large-scale rock slope. The change of groundwater levels was predicted within a relatively small error by making full use of site characterization data. The inversed unsaturated hydraulic parameters are comparable to existing results for rocks, fractures, and soils (Bandurraga and Bodvarsson, 1999; Flint, 2003; Hu et al., 2015). Without introducing sophisticated flow models and constitutive relations, the numerical results are of engineering accuracy and useful for, e.g., the design and performance assessment of impervious barriers (Li et al., 2017a; Zhou et al., 2018).
6. Conclusions This study investigated the groundwater flow behavior in a large-scale rock slope by means of inverse modeling using the continuum approach based on field observations spanning ten years at the Baihetan dam site, Southwest China. The inverse modeling showed that the unsaturated hydraulic parameters are in the range of
31
=0.014.65 m1, m=0.310.69 (n=1.453.23) and =0.530.79 for fractured rocks, loose sediments and geologic structures at the site, comparable to existing results for rocks and fractures at Yucca Mountain (Bandurraga and Bodvarsson, 1999; Flint, 2003) and/or soils of different types (Hu et al., 2015). The multiple water levels above and below the tuff zones and the drawdown of water table in the slope during site exploration were reasonably reproduced, and the estimated parameters are overall representative of the unsaturated hydraulic properties of rock formations at the site. It has been challenging to obtain the field scale unsaturated hydraulic parameters in fractured rocks. Here, it is demonstrated that without using the difficult-to-acquire flow and saturation data in the unsaturated zone, the unsaturated hydraulic parameters can be effectively obtained through the groundwater level and discharge data in the saturated zone. It is concluded that the inverse modeling procedure adopted is efficient and effective in seeking an optimal solution in a complex, large-scale inverse problem. Furthermore, this study boosts our confidence in the continuum approach for approximating the site-scale saturated-unsaturated flow behavior in fractured rock formations, an important but tough issue that is frequently encountered in various subsurface applications.
Acknowledgments
The financial supports from the National Natural Science Foundation of China (No.
51925906)
and
the
National
Key
R&D
Program
of
China
(No.
2018YFC0407001) are gratefully acknowledged. The authors gratefully thank POWERCHINA HUADONG Engineering Corporation Limited for the field data supporting this study.
32
Appendix A. Supplementary data
Supplementary data associated with this article can be found, in the online version. References Asher, M.J., Croke, B.F., Jakeman, A.J., Peeters, L.J., 2015. A review of surrogate models and their application to groundwater modeling. Water Resour. Res. 51(8), 5957-5973. https://doi.org/10.1002/2015WR016967. Assouline, S., Tessier, D., Bruand, A., 1998. A conceptual model of the soil water retention
curve.
Water
Resour.
Res.
34
(2),
223–231.
https://doi.org/10.1029/97WR03039. Bandurraga, T., Bodvarsson, G., 1999. Calibrating hydrogeologic parameters for the 3-D site-scale unsaturated zone model of Yucca Mountain, Nevada. J. Contam. Hydrol. 38, 25–46. https://doi.org/10.1016/S0169-7722(99)00010-8. Barenblatt, G.I., Zheltov, I.P., Kochina, I.N., 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]. J APPL MATH MEC. 24(5), 1286-1303. https://doi.org/10.1016/0021-8928(60)90107-6. Bathe, K.J., Khoshgoftaar, M.R., 1979. Finite element free surface seepage analysis without
mesh
iteration.
Int
J
Numer
Anal
Met.
3(1),
13–22.
https://doi.org/10.1002/nag.1610030103. Bertels, S.P., Dicarlo, D.A., Blunt, M.J., 2001. Measurement of aperture distribution, capillary pressure, relative permeability, and in situ saturation in a rock
33
fracture using computed tomography scanning. Water Resour. Res. 37 (3), 649–662. https://doi.org/10.1029/2000WR900316. Borsi, I., Farina, A., Primicerio, M., 2006. A rain water infiltration model with unilateral boundary condition: qualitative analysis and numerical simulations. Math. Methods Appl. Sci. 29, 2047–2077. https://doi.org/10.1002/mma.763. Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media. Hydrology Paper No. 3. Colorado State University, Fort Collins, CO. Burdine, N.T., 1953. Relative permeability calculations from pore-size distribution data.
Trans.
Am.
Inst.
Min.
Metall.
Pet.
Eng.
198,
71–77.
https://doi.org/10.2118/225-G. Carrera, J., Alcolea, A., Medina, A., Hidalgo, J., Slooten, L.J., 2005. Inverse problem in
hydrogeology.
Hydrogeol
J.
13
(1),
206–222.
https://doi.org/10.1007/s10040-004-0404-7. Chen, C., Horne, R.N., 2006. Two-phase flow in rough-walled fractures: Experiments and
a
flow
structure
model.
Water
Resour.
Res.
42,
W03430.
https://doi.org/10.1029/2004WR003837. Chen, Y.F., Hong, J.M., Zheng, H.K., Li, Y., Hu, R., Zhou, C.B., 2016. Evaluation of groundwater leakage into a drainage tunnel in Jinping-I arch dam foundation in Southwestern China: A case study. Rock Mech. Rock Eng. 49 (3), 961–979. https://doi.org/10.1007/s00603-015-0786-y. Chen, Y.F., Hu, S.H., Hu, R., Zhou, C.B., 2015. Estimating hydraulic conductivity of fractured rocks from high pressure packer tests with an Izbash’s law-based
34
empirical
model.
Water
Resour.
Res.
51,
2096–2118.
https://doi.org/10.1002/2014WR016458. Chen, Y.F., Ling, X.M., Liu, M.M., Hu, R., Yang, Z.B., 2018. Statistical distribution of hydraulic conductivity of rocks in deep-incised valleys, southwest China. J. Hydrol. 566, 216-226. https://doi.org/10.1016/j.jhydrol.2018.09.016. Chen, Y.F., Zhou, C.B., Jing, L., 2009. Modeling coupled THM processes of geological porous media with multiphase flow: theory and validation against laboratory and field scale experiments. Comput. Geotech. 36 (8), 1308–1329. https://doi.org/10.1016/j.compgeo.2009.06.001. Chen, Y.F., Zhou, C.B., Zheng, H., 2008. A numerical solution to seepage problems with complex drainage systems. Comput Geotech. 35(3), 383–393. https://doi.org/10.1016/j.compgeo.2007.08.005. Desai, C.S., Li, G.C., 1983. A residual flow procedure and application for free surface flow
in
porous
media.
Adv
Water
Resour.
6(1),
27–35.
https://doi.org/10.1016/0309-1708(83)90076-3. Doughty, C., 1999. Investigation of conceptual and numerical approaches for evaluating moisture, gas, chemical, and heat transport in fractured unsaturated rock.
J
Contam
Hydrol.
38(1-3),
69-106.
https://doi.org/10.1016/S0169-7722(99)00012-1. Eaton, R.R., McCord, J.T., 1995. Effective hydraulic conductivities in unsaturated heterogeneous media by Monte Carlo simulation. Transport Porous Med. 18(3), 203-216. https://xs.scihub.ltd/https://doi.org/10.1007/BF00616931.
35
Evans, D.D., Rasmussen, T.C., 1991. Unsaturated flow and transport through fractured rock related to high-level waste repositories, Final Report: Phase III. NUREG/CR-5581, 67 pp. https://www.osti.gov/servlets/purl/137966. Faybishenko, B., and Finsterle, S., 2000. Tensiometry in fractured rocks, in Zhang, D., Winter, C.L., eds., Theory, Modeling, and Field Investigation in Hydeogeology: A Special Volume in Honor of Shlomo P. Neuman’s 60th Birthday: Boulder, Colorado, Geological Society of America Special Paper 348, p. 161-174. Finsterle, S., 2000. Using the continuum approach to model unsaturated flow in fractured
rock.
Water
Resour.
Res.
36
(8),
2055–2066.
https://doi.org/10.1029/2000WR900122. Flint, L.E., 2003. Physical and hydraulic properties of volcanic rocks from Yucca Mountain,
Nevada.
Water
Resour.
Res.
39
(5).
https://doi.org/10.1029/2001WR001010. Fredlund, D.G., Xing, A.Q., 1994. Equations for the soil-water characteristic curve. Can. Geotech. J. 31, 521–532. https://doi.org/10.1139/t94-061. Gallipoli, D., 2012. A hysteretic soil-water retention model accounting for cyclic variations of suction and void ratio. Géotechnique 62 (7), 605–616. http://dx.doi.org/10.1680/geot.11.P.007. Garcia, L.A., Shigidi, A., 2006. Using neural networks for parameter estimation in ground
water.
J.
Hydrol.
https://doi.org/10.1016/j.jhydrol.2005.05.028.
36
318,
215–231.
Gong, W.Y., Cai, Z.H., Jiang, L.X., 2008. Enhancing the performance of differential evolution using orthogonal design method. Appl. Math. Comput. 206, 56–69. https://doi.org/10.1016/j.amc.2008.08.053. Gribb, M.M., 1996. Parameter estimation for determining hydraulic properties of a fine sand from transient flow measurements. Water Resour. Res. 32(7), 1965– 1974. https://doi.org/10.1029/96WR00894. Gribb, M.M., Šimůnek, J., Leonard, M.F., 1998. Development of cone penetrometer method to determine soil hydraulic properties. J Geotech Geoenviron. 124(9), 820–829. https://doi.org/10.1061/(ASCE)1090-0241(1998)124:9(820). Guarracino, L., 2006. A fractal constitutive model for unsaturated flow in fractured hard
rocks.
J.
Hydrol.
324,
154–162.
https://doi.org/10.1016/j.jhydrol.2005.10.004. Guo, T., Zhang, S., Zou, Y., Xiao, B., 2015. Numerical simulation of hydraulic fracture propagation in shale gas reservoir. J Nat Gas Sci Eng. 26, 847-856. https://doi.org/10.1016/j.jngse.2015.07.024. Hvorslev, M.L., 1951. Time lag and soil permeability in groundwater observations. Bull 36. Waterways Experiment Station Corps of Engineers. U. S. Army, Vicksburg, Mississippi. Hughson, D.L., Yeh, T.-C.J., 2000. An inverse model for three-dimensional flow in variably saturated porous media. Water Resour. Res. 36 (4), 829–839. https://doi.org/10.1029/2000WR900001. Hong, J.M., Chen, Y.F., Liu, M.M., Zhou, C.B., 2017. Inverse modelling of
37
groundwater flow around a large-scale underground cavern system considering the excavation-induced hydraulic conductivity variation. Comput Geotech. 81, 346-359. https://doi.org/10.1016/j.compgeo.2016.09.008. Hsieh, P.A., Neuman, S.P., 1985. Field determination of the three-dimensional hydraulic conductivity tensor of anisotropic media. Water Resour. Res. 21 (11), 1655–1665. https://doi.org/10.1029/WR021i011p01667. Hu, R., Chen, Y.F., Liu, H.H., Zhou, C.B., 2013. A water retention curve and unsaturated hydraulic conductivity model for deformable soils: consideration of the change in pore-size distribution. Géotechnique 63 (16), 13891405. https://doi.org/10.1680/geot.12.P.182. Hu, R., Chen, Y.F., Liu, H.H., Zhou, C.B., 2015. A relative permeability model for deformable soils and its impact on coupled unsaturated flow and elastoplastic deformation processes. Sci. China Tech. Sci. 58 (11), 1971–1982. https://xs.scihub.ltd/https://doi.org/10.1007/s11431-015-5948-3. Hu, R., Chen, Y.F., Zhou, C.B., Liu, H.H., 2017. A numerical formulation with unified unilateral boundary condition for unsaturated flow problems in porous media.
Acta
Geotech.
12
(2),
277291.
https://xs.scihub.ltd/https://doi.org/10.1007/s11440-016-0475-3. Illman, W.A., Hughson, D.L., 2005. Stochastic simulations of steady state unsaturated flow in a three-layer, heterogeneous, dual continuum model of fractured rock. J. Hydrol. 307, 17–37. https://doi.org/10.1016/j.jhydrol.2004.09.015. Illman, W.A., Neuman, S.P., 2000. Type-curve interpretation of multi-rate single-hole
38
pneumatic injection tests in unsaturated fractured rock. Groundwater 38(6), 899–911. https://doi.org/10.1111/j.1745-6584.2000.tb00690.x. Illman, W.A., Neuman, S.P., 2001. Type-curve interpretation of a cross-hole pneumatic test in unsaturated fractured tuff. Water Resour. Res. 37(3), 583– 604. https://doi.org/10.1029/2000WR900273. Illman, W.A., Neuman, S.P., 2003. Steady-state analyses of cross-hole pneumatic injection tests in unsaturated fractured tuff. J. Hydrol. 281(1-2), 36–54. https://doi.org/10.1016/S0022-1694(03)00199-9. Ireson, A.M., Wheater, H.S., Butler, A.P., Mathias, S.A., Finch, J., Cooper, J.D., 2006. Hydrological processes in the Chalk unsaturated zone–insights from an intensive field monitoring programme. J. Hydrol. 330(1-2), 29–43. https://doi.org/10.1016/j.jhydrol.2006.04.021. Jiang, Q., Feng, X.T., Hatzor, Y.H., Hao, X.J., Li, S.J., 2014. Mechanical anisotropy of columnar jointed basalts: An example from the Baihetan hydropower station,
China.
Eng.
Geol.
175,
35–45.
https://doi.org/10.1016/j.enggeo.2014.03.019. Karpouzos, D.K., Delay, F., Katsifarakis, K.L., de Marsily, G., 2001. A multipopulation genetic algorithm to solve the inverse problem in hydrogeology.
Water
Resour.
Res.
37,
2291–2302.
https://doi.org/10.1029/2000WR900411. Kool, J.B., Parker, J.C., 1988. Analysis of the inverse problem for transient unsaturated
flow.
Water
Resour.
39
Res.
24
(6),
817–830.
https://doi.org/10.1029/WR024i006p00817. Kueper, B.H., McWhorter, D.B., 1992. The use of macroscopic percolation theory to construct large-scale capillary pressure curves. Water Resour. Res. 28 (9), 2425–2436. https://doi.org/10.1029/92WR01176. Li, X., Chen, Y.F., Hu, R., Yang, Z.B., 2017a. Towards an optimization design of seepage control: A case study in dam engineering. Sci China Technol Sc. 60(12), 1903–1916. https://doi.org/10.1007/s11431-016-9160-y. Li, Y., Chen, Y.F., Zhang, G.J., Liu, Y., Zhou, C.B., 2017b. A numerical procedure for modeling the seepage field of water-sealed underground oil and gas storage caverns.
Tunn.
Undergr.
Space
Tech.
66,
5663.
https://doi.org/10.1016/j.tust.2017.04.002. Li, Y., Chen, Y.F., Zhou, C.B., 2014. Hydraulic properties of partially saturated rock fractures subjected to mechanical loading. Eng. Geol. 179, 2431. https://doi.org/10.1016/j.enggeo.2014.06.019. Liu, H.H., Bodvarsson, G.S., 2001. Constitutive relations for unsaturated flow in fracture
network.
J.
Hydrol.
252,
116–125.
https://doi.org/10.1016/S0022-1694(01)00449-8. Liu, H.H., Bodvarsson, G.S., 2003. Upscaling of constitutive relations in unsaturated heterogeneous
tuff
matrix.
J.
Hydrol.
276,
198–209.
https://doi.org/10.1016/S0022-1694(03)00071-4. Liu, H.H., Doughty, C., Bodvarsson, G.S., 1998. An active fracture model for unsaturated flow and transport in fractured rocks. Water Resour. Res. 34 (10),
40
2633-2646. https://doi.org/10.1029/98WR02040. Liu, H.H., Haukwa, C.B., Ahlers, C.F., Bodvarsson, G.S., Flint, A.L., Guertal, W.B., 2003. Modeling flow and transport in unsaturated fractured rock: an evaluation of
the
continuum
approach.
J.
Contam.
Hydrol.
62,
173–188.
https://doi.org/10.1016/S0169-7722(02)00170-5. Lu, Z., Kwicklis, E.M., 2012. Numerical evaluation of effective unsaturated hydraulic properties of fractured rocks using a stochastic continuum approach. Vadose Zone J. 11 (4). https://doi.org/10.2136/vzj2011.0164. Mantoglou, A., Gelhar, L.W., 1987a. Stochastic modeling of large-scale transient unsaturated
flow
systems.
Water
Resour.
Res.
23(1),
37-46.
https://doi.org/10.1029/WR023i001p00037. Mantoglou, A., Gelhar, L., 1987b. Capillary tension head variance, mean soil moisture content, and effective specific soil moisture capacity of transient unsaturated flow in stratified soils. Water Resour. Res. 23(1), 47-56. https://doi.org/10.1029/WR023i001p00047. Mantoglou, A., Gelhar, L.W., 1987c. Effective hydraulic conductivities of transient unsaturated flow in stratified soils. Water Resour. Res. 23(1), 57-67. https://doi.org/10.1029/WR023i001p00057. Mao, D., Wan, L., Yeh, T.-C.J., Lee, C.-H., Hsu, K.-C., Wen, J.-C., Lu, W., 2011. A revisit of drawdown behavior during pumping in unconfined aquifers. Water Resour. Res. 47(5). https://doi.org/10.1029/2010WR009326. Mathias, S.A., Butler, A.P., Jackson, B.M., Wheater, H.S., 2006. Transient
41
simulations of flow and transport in the Chalk unsaturated zone. J. Hydrol. 330, 10–28. https://doi.org/10.1016/j.jhydrol.2006.04.010. Montazer, P., Wilson, W.E., 1984. Conceptual hydrological model of flow in the unsaturated zone, Yucca Mountain, Nevada. USGS Water-Resources Investigation Report, 84-4345, 55 pp. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous
media.
Water
Resour.
Res.
12
(3),
513–522.
https://doi.org/10.1029/WR012i003p00513. National Research Council, 1996. Rock Fractures and Fluid Flow-Contemporary Understanding and Applications. National Academy Press, Washington D.C. Nuth, M., Laloui, L., 2008. Advances in modelling hysteretic water retention curve in deformable
soils,
Comput.
Geotech.,
35
(6),
835–844.
https://doi.org/10.1016/j.compgeo.2008.08.001. Perkins, K., Johnson, B.D., Mirus, B.B., 2014. Measurement of unsaturated hydraulic properties and evaluation of property-transfer models for deep sedimentary interbeds, Idaho National Laboratory, Idaho. Report 2328–0328. Persoff, P., Pruess, K., 1995. Two-phase flow visualization and relative permeability measurement in natural rough-walled rock fractures. Water Resour. Res. 31 (5), 1175–1186. https://doi.org/10.1029/95WR00171. Peters, R.P., Klavetter, E.A., 1988. A continuum model for water movement in an unsaturated fractured rock mass. Water Resour. Res. 24 (3), 416–430. https://doi.org/10.1029/WR024i003p00416.
42
Pruess, K., Tsang, Y.W., 1990. On two-phase relative permeability and capillary pressure of rough-walled rock fractures. Water Resour. Res. 26, 1915–1926. https://doi.org/10.1029/WR026i009p01915. Rafiai, H., Jafari, A., Mahmoudi, A., 2013. Application of ANN-based failure criteria to rocks under polyaxial stress conditions. Int. J. Rock Mech. Min. Sci. 59, 42–49. https://doi.org/10.1016/j.ijrmms.2012.12.003. Reitsma, S., Kueper, B. H., 1994. Laboratory measurement of capillary pressure-saturation relationships in a rock fracture. Water Resour. Res. 30, 865–878. https://doi.org/10.1029/93WR03451. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics 1, 318–333. https://doi.org/10.1063/1.1745010. Rulon, J.J., Rodway, R., Freeze, R.A., 1985. The development of multiple seepage faces
on
layered
slopes.
Water
Resour.
Res.
21,
1625–1636.
https://doi.org/10.1029/WR021i011p01625. Salem, M.B., Tomaso, L., 2018. Automatic selection for general surrogate models. Struct
Multidiscip
O.
58(2),
719-734.
https://xs.scihub.ltd/https://doi.org/10.1007/s00158-018-1925-3. Song, X., Shi, L., Ye, M., Yang, J., Navon, I.M., 2014. Numerical comparison of iterative ensemble Kalman filters for unsaturated flow inverse modeling. Vadose Zone J. 13(2). https://doi.org/10.2136/vzj2013.05.0083. Tarantino, A., 2009. A water retention model for deformable soils. Géotechnique 59 (9), 751–762. https://doi.org/10.1680/geot.7.00118.
43
van Genuchten, M., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soil. Soil Sci. Soc. Am. J. 44 (5), 892–898. https://doi.org/10.2136/sssaj1980.03615995004400050002x. Wang, J.S.Y., Narasimhan, T.N., 1985. Hydrologic mechanisms governing fluid flow in a partially saturated, fractured porous medium. Water Resour. Res. 20 (12), 1861–1874. https://doi.org/10.1029/WR021i012p01861. Wang, X.S., Jiang, X.W., Wan, L., Song, G., Xia, Q., 2009. Evaluation of depth-dependent
porosity
and
bulk
modulus
of
a
shear
using
permeability-depth trends. Int. J. Rock Mech. Min. Sci. 46, 1175–1181. https://doi.org/10.1016/j.ijrmms.2009.02.002. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE. J. 3(03), 245-255. https://doi.org/10.2118/426-PA. Weerakone, W., Wong, R.C.K., Mehrotra, A.K., 2012. Measurement of capillary pressure curve of DNAPL in a water-saturated sandstone fracture. J. Geotech. Geoenviron.
Eng.
138,
614–624.
https://doi.org/10.1061/(ASCE)GT.1943-5606.0000590. Wu, Y.S., Liu, H.H., Bodvarsson, G.S., 2004. A triple-continuum approach for modeling flow and transport processes in fractured rock. J. Contam. Hydrol. 73, 145–179. https://doi.org/10.1016/j.jconhyd.2004.01.002. Wu, Y.S., Ritcey, A.C., Bodvarsson, G.S., 1999. A modeling study of perched water phenomena in the unsaturated zone at Yucca Mountain. J. Contam. Hydrol. 38, 157–184. https://doi.org/10.1016/S0169-7722(99)00015-7.
44
Wu, Y.S., Zhang, K., Liu, H.H., 2006. Estimating large-scale fracture permeability of unsaturated rock using barometric pressure data. Vadose Zone J. 5(4), 1129– 1142. https://doi.org/10.2136/vzj2006.0015. Yan, C., Jiao, Y.Y., 2018. A 2D fully coupled hydro-mechanical finite-discrete element model with real pore seepage for simulating the deformation and fracture of porous medium driven by fluid. Comput Struct. 196, 311-326. https://doi.org/10.1016/j.compstruc.2017.10.005. Yeh, T.C.J., Gelhar, L.W., Gutjahr, A.L., 1985a. Stochastic analysis of unsaturated flow in heterogeneous soils: 1. Statistically isotropic media. Water Resour. Res. 21(4), 447-456. https://doi.org/10.1029/WR021i004p00447 Yeh, T.C.J., Gelhar, L.W., Gutjahr, A.L., 1985b. Stochastic analysis of unsaturated flow in heterogeneous soils: 2. Statistically anisotropic media with variable α, Water
Resour.
Res.
21(4),
457-464.
https://doi.org/10.1029/WR021i004p00457. Yeh, T.C.J., Gelhar, L.W., Gutjahr, A.L., 1985c. Stochastic analysis of unsaturated flow in heterogeneous soils: 3. Observations and applications. Water Resour. Res. 21(4), 465-471. https://doi.org/10.1029/WR021i004p00465. Younes, A., Mara, T.A., Voltz, M., Guellouz, L., Baalousha, H.M., Fahs, M., 2018. A new efficient Bayesian parameter inference strategy: Application to flow and pesticide transport through unsaturated porous media. J. Hydro. 563, 887–899. https://doi.org/10.1016/j.jhydrol.2018.06.043. Zhang, M., Hu, L., Yao, L., Yin, W., 2017. Surrogate models for sub-region
45
groundwater management in the beijing plain, China. Water-Sui. 9(10), 766. https://doi.org/10.3390/w9100766. Zhang, L.M., Fredlund, D.G., 2003. Characteristics of water retention curves for unsaturated fractured rocks. Second Asian Conference on Unsaturated Soils (UNSAT-ASIA 2003), Osaka, Japan, p.425–428. Zheng, H., Liu, D.F., Lee, C.F., Tham, L.G., 2005. A new formulation of Signorini's type for seepage problems with free surfaces. Int J Numer Meth Eng. 64(1), 1– 16. https://doi.org/10.1002/nme.1345. Zhou, A., Sheng, D., Carter, J., 2012. Modelling the effect of initial density on soil-water
characteristic
curves.
Géotechnique
62
(8),
669–680.
https://doi.org/10.1680/geot.10.P.120. Zhou, C.B., Liu, W., Chen, Y.F., Hu, R., Wei, K., 2015. Inverse modeling of leakage through a rockfill dam foundation during its construction stage using transient flow model, neural network and genetic algorithm. Eng. Geol. 187, 183–195. https://doi.org/10.1016/j.enggeo.2015.01.008. Zhou, C.B., Zhao, X.J., Chen, Y.F., Liao, L., Liu, M.M., 2018. Interpretation of high pressure pack tests for design of impervious barriers under high-head conditions.
Eng.
Geol.
234,
112–121.
https://doi.org/10.1016/j.enggeo.2018.01.006. Zhou, H., Gómez-Hernández, J.J., Li, L., 2014. Inverse methods in hydrogeology: Evolution
and
recent
trends.
Adv.
Water
https://doi.org/10.1016/j.advwatres.2013.10.014.
46
Resour.
63,
22–37.
Zimmerman, R.W., Chen, G., Hadgu, T., Bodvarsson, G.S., 1993. A numerical dual-porosity model with semianalytical treatment of fracture/matrix flow. Water Resour. Res. 29(7), 2127-2137. https://doi.org/10.1029/93WR00749.
47
48
List of Figures (12 figures in total)
Fig.1. Location and landscape of the Baihetan dam site.
Precipitation and evaporation rate (mm)
200 Precipitation Evaporation
150 100 50 0 -50 -100 1
2
3
4
5
6
7 Month
8
9
10
11
12
Fig.2. Mean monthly precipitation measured at Baihetan weather station (635 m a.s.l.) and mean monthly evaporation from the land surface in the study area.
49
Fig.3. Geologic cross-section I-I at the dam site (with its trace shown in Fig. 1): (a) regional view and (b) enlarged view for the near-bank region on the left bank.
50
Fig.4. Groundwater flow system and locations of boreholes and exploratory adits.
Fig.5. Schematic diagram of groundwater level observation at boreholes penetrating a tuff zone. The drilling of boreholes across the tuff zone resulted in drawdown of groundwater level above it. The groundwater levels were measured after a concrete plug was constructed at the intersection between the borehole and the tuff zone, and the vertical distance between the upper and lower groundwater levels were then calculated. 51
Fig.6. 2D finite element mesh for the study area: (a) overall view, and (b) enlarged view for the near-bank region.
52
Fig.7. Sketch of the two hidden-layer BPNN model. The input layer contains 19 neurons for input of the unsaturated hydraulic parameters P and the time-dependent precipitation/evaporation boundary conditions. The output layer contains 12 neurons for output of the time series data of groundwater level in 9 boreholes and discharge at 3 weirs. The numbers of neurons at the two hidden layers (22 and 35) were determined by minimizing an error function on a test data set with a trial-and-error method.
Fig.8. The estimated water retention curves for rocks in different rock zones. 53
(a) 770
Measured water level Predicted water level Precipitation
760
Water level (m)
750
760 750
ZK1123
740
740 730
730 720
720
ZK1119
710
710
700
700
ZK9319
690 680 640
Precipitation (mm)
770
690 680 640
ZK929
630
630
620
620
610
610
150
150
100
100
50
50
0 2009/12/15 2010/02/15
2010/05/15
2010/08/15
Time
2010/11/15
0 2011/02/15
(b)
750
ZK512
Water level (m)
740
Measured water level Predicted water level Precipitation
730
ZK412
740 730
720
720
710
710
700
700
ZK312
690
Precipitation (mm)
750
690
680
680
670
670
150
150
100
100
50
50
0
0
2005/12/15 2006/04/15 2006/08/15 2006/12/15 2007/04/15 2007/08/15
Time
54
(c)
Precipitation (mm)
Water level (m)
748
Water level (m) Precipitation (mm)
748
744
744
ZK37 726
726
724
724
722
722
720
720
150
150
100
100
50
50
0
2001/12/15
(d)
Measured water level Predicted water level Precipitation
2002/03/15
2002/06/15
2002/09/15
Time
0
2002/12/15 2003/02/15
Measured water level Predicted water level Precipitation
780 774
780 774
768
768
762
762
ZK112
756
756
710 700 690 680
710 700 690 680
150
150
100
100
50
50
0
2005/09/15 2006/03/15
2006/11/15
2007/07/15
Time
2010/05/15
0
2010/11/15
Fig.9. Comparison of the observed (dot points) and calculated (solid lines) groundwater levels at observation boreholes in the slope. The rainfall intensity is also plotted in each figure.
55
Fig. 10. Water pressure head contours (a) in wet season (on July 15, 2000) and (b) in dry season (on January 15, 2001).
56
Fig.11. Water head contours on May 15, 2011 after the exploratory adits were excavated.
300
200
2.0
May, 2011
1.5 150 109
100
0.64
50 0
Excavation began May, 2002
2002
2004
1.0 0.5
Flow Rate (L/min)
250
Flow Rate (L/min)
2.5
Prediction Measurement PD61 PD61-3 PD61-4 PD61-5
0.036
2006 2008 Time (year)
2010
0.0 2012
Fig.12. Variation of discharge into the exploratory adits during the exploration period. Also plotted are the available measurements in May 2011.
57
58
List of Tables (5 tables in total)
Table 1 Specific storage and hydraulic conductivity for rocks at the site. Statistics of hydraulic conductivity (K) Standard Rock No. of Range (m/s) Mean (m/s) deviation tests (m/s) Loose sediments 4.31×106 1.41×1048.20×104 7.00104 2.40104 63 6 6 5 5 High-moderate PZ 2.94×10 5.69×10 1.94×10 1.00×10 3.48×106 8 Moderate PZ 1.91×106 1.15×1071.29×105 3.00106 1.18106 128 Moderate-low PZ 1.53×106 1.55×1089.32×106 9.00107 4.96107 561 6 9 6 7 7 Low PZ 1.53×10 9.06×10 4.14×10 3.0010 3.2210 1603 Very low PZ 9.42×107 6.47×1091.26×106 8.00108 9.39108 694 Fault F17 5.78×106 5.18×1071.61×105 2.75106 1.95107 11 Faults F129, F133, F134* 4.22×106 / 2.00106 / / * 6 8 6 7 7 Tuff zone C2 5.55×10 1.8110 1.1010 2.2010 4.4310 31 Tuff zone C3(A) * 5.55×106 1.811083.40106 3.01×107 4.07107 34 Tuff zone C3(B) * 5.55×106 3.101071.30×105 2.98106 / / Tuff zone C3-1(A) * 5.55×106 1.681082.70106 4.50×107 9.07107 33 Tuff zone C3-1(B) * 5.55×106 1.601067.20106 1.68106 / / 6 8 5 Shear zone LS331 5.55×10 2.3310 5.5010 8.00107 2.35106 23 Shear zone LS3318 5.55×106 3.621086.00106 8.15107 2.26107 8 Shear zone LS3319 5.55×106 7.121081.20106 3.73107 3.14107 18 * The tuff zones C and C 3 3-1 are divided into an unweathered subzone A and a weathered subzone B. The hydraulic conductivity of C2, C3(A) and C3-1(A) is highly anisotropic. The listed K values of C2, C3(A) and C3-1(A) represent the in-plane hydraulic conductivity (K//), and the values of K Specific storage Ss (m1)
perpendicular to the plane (K) are about two orders of magnitude smaller than the in-plane values (i.e., K=K///100). The hydraulic conductivity of all other rock zones is assumed to be isotropic. The hydraulic conductivity of the tuff zones C3(B) and C3-1(B) and the small-scale faults F129, F133 and F134 were estimated by hydraulic tests other than the packer tests.
Table 2 Comparison of the observed and calculated vertical distance between the upper and lower groundwater levels in wet season (on June 18, 2000). Borehole ZK925 ZK1119 ZK1121 ZK412
Distance to section I-I (m)
Tuff zone penetrated
0 111.6 111.6 259.6
C3 C3-1 C3-1 C3 59
Vertical distance of groundwater levels (m) Observed
Calculated
13.20 24.40 38.90 22.99
9.68 16.26 23.74 19.68
Table 3 Unsaturated hydraulic parameters for rocks at the site. Rock
θs
θr
Loose sediments High-moderate/moderate PZ Moderate-low/low PZ Very low PZ Tuff and shear zones Faults
0.32 0.20 0.12 0.08 0.25 0.30
0.05 0.03 0.02 0.01 0.04 0.05
α (m1) Range Estimated 0.110 0.600 0.0010.1 0.034 0.0010.1 0.012 0.0010.1 0.010 0.110 3.033 0.110 4.646
m Range Estimated 01 0.437 01 0.430 01 0.425 01 0.305 01 0.689 01 0.634
β Range Estimated 03 0.533 03 0.580 03 0.611 03 0.756 03 0.782 03 0.792
Table 4 Groundwater level observations in nine boreholes used for inverse modelling. Borehole
Distance to section I-I (m)
Collar elevation (m)
Depth (m)
Groundwater level (m)
Water level fluctuation (m)
Observation period
ZK929 ZK9319 ZK1123 ZK1119 ZK312 ZK412 ZK37 ZK112 ZK512
0 49.9 111.6 111.6 205.8 259.6 309.2 309.2 361.9
695.45 712.63 853.04 840.96 726.81 877.65 873.71 859.50 883.23
130.72 181.10 200.89 200.00 140.28 377.81 200.90 120.20 165.40
624.3620.6 687.0685.3 750.1745.9 717.3712.8 686.5677.5 715.0707.1 745.4741.9 767.7752.0 740.7733.2
3.7 1.7 4.2 4.5 9.0 7.9 3.5 15.7 7.5
2010/022011/01 2010/052011/01 2010/052011/01 2010/052011/01 2006/062007/07 2006/082007/10 2004/012005/01 2010/022011/01 2005/122007/10
Table 5 Comparison of the measured and calculated discharges in May 2011 at three weirs installed in the exploratory adit system. No.
Location of weir
Excavation time of adit
1 2 3 4 5 6
PD61 PD61-1 PD61-2 PD61-3 PD61-4 PD61-5
2002/052006/05 2006/032006/05 2006/052006/06 2006/072006/08 2006/092006/10 2009/012009/05
60
Discharge in May 2011 (L/min) Measured
Calculated
109 N/A N/A 0.036 N/A 0.64
113.08 0 0 0.041 0.259 0.841
61
evaporation
C3
46% 54%
C3-1
faults
65% 35%
C3
tuff zones
C3-1
multiple wa ter tables
river 50 m
Graphical Abstract: We investigate the saturated-unsaturated flow in a large-scale rock slope at the Baihetan dam site, Southwest China, using the continuum approach and inverse modeling. The unsaturated hydraulic parameters in water retention models are determined based on the time-series observations of groundwater level and the available discharge measurements. The typical features of the groundwater flow mainly governed by gently sloping tuff zones of low vertical permeability, such as multiple water tables, wedge-shaped unsaturated zones and significant groundwater level drawdown during ten-year site characterization, are well simulated.
62
63
Highlights
A methodology was presented for field-scale inverse modeling of unsaturated flow (82/85)
The determined unsaturated parameters are well representative of rocks at the site (84/85)
Numerical results reproduce the main features of multiple water tables at the site (84/85)
64