Modeling unsaturated flow in fractured rocks with scaling relationships between hydraulic parameters

Yi-Feng Chen,Yuke Ye,Ran Hu*,Zhibing Yang,Chuang-Bing Zhou,c

a State Key Laboratory of Water Resources and Hydropower Engineering Science Wuhan University,Wuhan,430072,China

b Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering of the Ministry of Education,Wuhan,430072,China

c School of Civil Engineering and Architecture,Nanchang University,Nanchang,330031,China

Keywords:Unsaturated flow van genuchten model Hydraulic properties Fractured rocks Continuum approach

ABSTRACT Modeling unsaturated flow in fractured rocks is essential in various subsurface engineering applications,but it remains a great challenge due to the difficulties in determining the unsaturated hydraulic properties of rocks that contain various scales of fractures.It is generally accepted that the van Genuchten(VG) model can be applied to fractured rocks,provided that the hydraulic parameters could be representatively determined.In this study,scaling relationships between the VG parameters (α and n) and hydraulic conductivity (K) across 8 orders of magnitude,from 10-10 m/s to 10-2 m/s,were proposed by statistical analysis of data obtained from 1416 soil samples.The correlations were then generalized to predict the upper bounds of VG parameters for fractured rocks from the K data that could be obtained more easily under field conditions,and were validated against a limited set of data from cores,fractures and fractured rocks available in the literature.The upper bound estimates significantly narrow the ranges of VG parameters,and the representative values of α and n for fractured rocks at the field scale can then be determined with confidence by inverse modeling using groundwater observations in saturated zones.The proposed methodology was applied to saturated-unsaturated flow modeling in the right-bank slope at the Baihetan dam site with a continuum approach,showing that most of the flow behaviors in fractured rocks in this complex hydrogeological condition could be properly reproduced.The proposed method overcomes difficulties in suction measurement in fractured rocks with strong heterogeneity,and provides a feasible way for modeling of saturated-unsaturated flow in fractured rocks with acceptable engineering accuracy.

1.Introduction

Unsaturated flow is a process of great importance in vadose zones where water moves from land surface to aquifer.Unsaturated flow occurs not only in soil covers,but also commonly in rock slopes or even in the surrounding rocks of deep tunnels (Liu et al.,1998;Pirastru and Niedda,2010;Chen et al.,2020;Zhou et al.,2021).Over the past decades,there has been an increasing interest in modeling the unsaturated flow in fractured rocks for a wide range of subsurface applications such as landslide mitigation (Rulon et al.,1985; Li et al.,2014),tunnel inflow prediction (Zhou et al.,2021),nuclear waste disposal (Liu et al.,1998; Mukhopadhyay et al.,2006),contaminant transport (Farouk et al.,2015),and underground oil storage(Li et al.,2017).A variety of approaches have been developed for modeling the unsaturated flow behaviors through fractured media,including the continuum (Peters and Klavetter,1988;Finsterle,2000),dual-continuum (McLaren et al.,2000; Mathias et al.,2006),triple-continuum (Wu et al.,2004),discrete fracture network (Evans and Rasmussen,1991;Therrien and Sudicky,1996),percolation network (Kueper and McWhorter,1992),and stochastic representation approaches (Yang et al.,2004).In engineering practices,the continuum approach is still the dominant modeling approach due to its simplicity and high computational efficiency.This approach applies on conditions that two important constitutive relations,i.e.the water retention curve and the relative permeability function,could be properly determined for the media.

For soils,various techniques have been developed to obtain the water retention curve by measuring the matrix or total suction of soil samples,with some of them even applicable under field conditions.A number of water retention relation models,such as Gardner (1956) model,Brooks and Corey (1964) model,van Genuchten (VG) (1980) model and Fredlund and Xing (1994) model,were also developed.These models have been widely adopted for curve fitting to the test data,and to estimate the relative permeability of soils based on Burdine’s(1953)or Mualem’s(1976)model.Statistical and empirical correlations were proposed to predict the water retention curve or its model parameters from soil properties such as mineralogical composition or pore size distribution(Gupta and Larson,1979;Rawls et al.,1982;Saxton et al.,1986;Rawls and Brakensiek,1989),by synthesizing the abundant experimental data of soil moisture and suction that has been accumulated for decades.

Fractured rocks are typically characteristic of multi-scale,multidimensional void geometries with high heterogeneity,having a much larger representative elementary volume(REV)size than that of soils.The local equilibrium condition could hardly be guaranteed within the REVs; therefore,the suction measurement techniques originally developed for soils fail to apply for fractured rocks at the field scale.How to determine a proper water retention relation for fractured rocks remains a challenging issue in subsurface hydrology.As validated by numerical upscaling and inverse modeling,the VG model is believed to apply for macroscopically characterizing the water retention behaviors of fractured rocks,provided that the model parameters could be representatively determined (Peters and Klavetter,1988; Liu et al.,1998; Chen et al.,2020).

The goal of this study is to provide a feasible way to estimate the VG parameters for fractured rocks by statistical correlation and inverse modeling.The statistical correlations are proposed by synthesizing data from soils,rocks,fractures,and fractured rocks available in the literature,and are used to predict the upper bounds of VG parameters from the more easily-obtained permeability data for fractured rocks.Inverse modeling is then invoked to determine the representative values of the VG parameters for fractured rocks,using predicted bounds and groundwater observations in saturated zones as inputs.The saturated-unsaturated flow in the right bank slope at the Baihetan dam site is studied with the proposed method,showing that the complex flow behaviors such as multiple water tables,wedge-shaped unsaturated zones,and preferential flow in backbone structures could be satisfactorily described.The proposed method features a rather narrow bound for the estimation of VG parameters and a feasible way of data collection for inverse modeling,making it practical for unsaturated flow modeling in highly heterogeneous fractured rocks with acceptable engineering accuracy.

2.Scaling between VG parameters and hydraulic conductivity

2.1.The VG model

The VG model(van Genuchten,1980)is one of the most widely used relations for describing the water retention behavior in porous media,which reads:

where Seis the effective saturation; θ is the volumetric water content;θsis the saturated water content;θris the residual water content; s is the suction; and α,n,and m=1-1/n are the VG parameters.

Based on Mualem’s (1976) model,the VG permeability relation for the media at unsaturated state is known as

where kris the relative permeability.

The VG model was originally proposed for soils(van Genuchten,1980),and its applicability in fractured rocks was later examined by numerical upscaling or inverse modeling (Peters and Klavetter,1988; Liu et al.,1998; Chen et al.,2020).However,it remains extremely difficult in determining the VG parameters(α and n)for fractured rocks with strong heterogeneity,where local equilibrium could hardly be ensured by the available field techniques at present.Given the fact that both VG parameters (α and n) and hydraulic conductivity (K) are governed by the void structure of fractured rocks,empirical correlations could be established to obtain interval estimates for the parameters.As long as the confidence intervals are sufficiently narrow,the values of α and n representative of site conditions could be much easily determined,e.g.by inverse modeling.

2.2.Statistics of VG parameters

A set of data (α,n,and K) on 1416 soil samples was collected from unsaturated soil hydraulic database(UNSODA)and published papers(see Table A1),covering all the 12 kinds of soil textures from sand to clay in the United States Department of Agriculture(USDA)classification(Fig.1).The data were obtained both in laboratory and in field.The dimensions of laboratory samples varied from 5.4 cm × 3 cm (diameter × height) to 25 cm × 25 cm × 30 cm(length × width × height),and the dimensions of field tests spanned from 1 m×1 m to 10 m×10 m(length×width).The θ-s relationships of the soil samples were obtained by the techniques including tensiometer,pressure plate,Tempe cell,electrical conductivity sensor,and centrifuge permeameter.In the collected dataset,there are at least five pairs (θ,s) for each soil sample,and hence the VG parameters can be well determined by curve fitting to the test data with the nonlinear least-squares method.The saturated hydraulic conductivity (K) of the soil samples varies in eight orders of magnitude from 10-10m/s to 10-2m/s.The data sources are listed in Table A1.

Fig.1.Classification of four soil groups.

The α value of the soil samples varies from 0.001 kPa-1to 4.789 kPa-1,while the n value in the range of 1.01-9.631.Table 1 lists the statistics of the VG parameters for four classes of soils merged from the 12 soil textures,showing that the geometric means of α and n increase as the soil particles become coarser and the sand proportion becomes higher.This scenario can be explained by higher adsorptive force and capillary force existing in smaller pores and clays(Miller et al.,2002;Aubertin et al.,2003),or alternatively,by the increasing heterogeneity from clays to sands.The VG parameters of soils are typical of high variability,with the variation coefficient greater than 1 and 0.35 for α and n,respectively.

Table 1 Statistics of VG parameters for soils.

Carsel and Parrish (1988) proposed that the VG parameters are normally distributed after Johnson transformation.Russo and Bouton (1992) concluded that α and n follow a lognormal and a normal distribution,respectively.Fig.2 plots the histograms of the α and n values of the soil samples,where the α values are transformed by the Box-Cox method.Using hypothesis test and Anderson-Darling test,it is found that the Box-Cox transformed α value is normally distributed (Fig.2a),while n follows a threeparameter lognormal distribution (Fig.2b).This conclusion also holds for each class of soil textures(Figs.A1 and A2),as well as the data of α and n within every 0.2 interval of the log-transformed K data(log10K)in the range from 10-7m/s to 10-4m/s(Figs.A3 and A4).

2.3.Upscaling relationships between VG parameters and K

It is recognized that the water retention curve is strongly correlated to the particle size distribution or pore size distribution of soils (Arya and Paris,1981; Fredlund et al.,1997).For the VG model,the parameter α approximately takes the inverse of soil’s air-entry value (van Genuchten and Nielsen,1985; Lenhard et al.,1989; Morel-seytoux et al.,1996; Tinjum et al.,1997),indicating that α is roughly proportional to the pore radius according to Young-Laplace equation.The parameter n characterizes the slope of the transition zone on the water retention curve,and a broader distribution of particle sizes generally leads to a broader distribution of pore sizes and a smaller value of n (Sillers,2001; Benson et al.,2014).

The hydraulic conductivity(K)is related to the porosity and pore size distribution of the media as well as the interconnectivity of pores.Denote the characteristic size of pores in the media by l,one has K ∝l2.According to Young-Laplace equation,there holds approximately α ∝l,which immediately yields the relation α ∝K0.5that can also be inferred from the capillary model (Mishra and Parker,1990).However,this relation does not always hold due to oversimplification of flow geometries.For instance,the linear relation s ∝ K-0.5is not valid for soils with smaller sand proportions (Yang et al.,2013).This correlation can thus be generalized into the power law model:

Fig.2.Histogram and probability density functions (a) For Box-Cox transformed α;and (b) For n.

where λ is the exponent and ? is the scaling coefficient.

The collected soil data in the α ～K and n ～K planes are plotted in Fig.3.In spite of a large dispersion,there is a clear trend that as K increases,the values of α and n increase with increasing variability.The arithmetic mean of α values averaged at each 0.2-interval of log10K can be well correlated to K with R2=0.97 (Fig.3a):

Similarly,the arithmetic mean of n values averaged at each 0.2-interval of log10K can be best fitted with the following function with R2=0.9 (Fig.3b):

Fig.3.Variation of (a) α and (b) n of soils with K.

Eqs.(4)and(5)indicate that as K increases from 1×10-10m/s to 1 ×10-2m/s,the mean of α value increases from 0.024 kPa-1to 2.372 kPa-1,and the mean of n value increases from 1.5 to 6.Unless otherwise noted,K and n are in units of m/s and kPa-1,respectively.

Fig.4 shows the standard deviations of α and n values calculated at each 0.2-interval of log10K,which also increase with K and can be best described with the following power-law correlations:

where σαand σnare the standard deviations of the α and n values,respectively,with σαin units of kPa-1.The best fitted values of the exponent for the above two equations are both 0.15,with R2=0.76 and 0.81,respectively.

Fig.4.Variation of standard deviation (a) For α (σα); and (b) For n (σn) with K.

It is not surprising to see from Eq.(4)-(7) that both the arithmetic means and the standard deviations of the α and n values for soils increase with the enhancement of K across eight orders of magnitude from 10-10m/s to 10-2m/s.As K increases,both the pore size and the heterogeneity of soils increase,the former leads to an increase of the means (and),while the latter causes an increase of the standard deviations(σαand σn).

2.4.Upper bound estimates for VG parameters

According to Chebyshev’s theorem,the upper bounds of α and n can be predicted from K as

where αUBand nUBare the upper bound estimates of α and n.

The lower bounds of α and n are known as αLB=0 and nLB=1,hence Eqs.(8) and (9) give simple confidence intervals for VG parameters.Fig.3 plots the curves of αUBand nUBfor soils,showing that the majority of the data points fall below the upper bounds.Owing to higher heterogeneity of pore structures,the interval between the lower and upper bounds for both α and n increases with the increasing value of K for higher permeability media.

Fig.5.Variation of (a) α and (b) n of fractured media with K.

Although Eqs.(8) and (9) are established based on soil’s data,they are supposed to apply for fractured rocks.Fig.5 plots a limited set of (α,n,K) data on 115 rock cores,fractures,and fractured formations available in the literature(e.g.Liu et al.,2003;Pirastru and Niedda,2010;Chen et al.,2020),including the data of tuff cores and fractures at Yucca Mountain (BSC,2000).The lithology includes extrusive rocks(e.g.basalt),sedimentary rocks(e.g.limestone),and metamorphic rocks (e.g.tuff,mylonite and slate).The data was obtained by laboratory tests,field tests,numerical upscaling or inverse modeling (Table A1).

Table 2 The third-type boundary conditions.

Fig.5 shows again that the majority of data falls below the upper bounds,verifying the applicability of Eqs.(8) and (9) in fractured rocks.Within the diverse range of K between 1 ×10-10m/s and 1 ×10-2m/s,the values of α for rock cores are mostly below the mean curve for soils predicted by Eq.(4) while the values of n for rock cores are mostly around the mean curve (Eq.(5)).A higher variability of data is observed for fractures and fractured rocks,with the values of α and n mostly lying above the mean curves.This implies that the VG parameters become more variable as the heterogeneity increases from cores to fractures and to fractured rocks.

Eqs.(8) and (9) are of practical significance to predict a proper range of VG parameters from the value of K for highly heterogeneous media.This is because,even in unsaturated zones,the hydraulic conductivity can be easily and reliably characterized by a variety of techniques that have been well developed for decades.The predicted confidence intervals are sufficiently narrow,hence significantly reducing the uncertainties and difficulties in determining the representative values of α and n by,e.g.the inverse modeling technique or simply the trial and error method.

3.Unsaturated flow model in fractured formations

3.1.Saturated-unsaturated flow model

The continuum approach is the most efficient and widely used approach for groundwater flow modeling,which produces simulation results of sufficient accuracy for a variety of engineering applications.In this framework,the saturated-unsaturated flow through the formations is governed by Richards (1931) equation:

where h is the pressure head;Cw=-?θ/?h is the specific moisture capacity,i.e.the slope of θ-s curve(s in units of m);Ssis the specific storage; ω is an indicator 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-Buckingham’s law:

where K is the saturated hydraulic conductivity tensor,z is the positive upward vertical coordinate,and h+z is the hydraulic head.

The initial condition is given by:

where h0is the initial pressure head in the domain Ω.

The boundary conditions for Eq.(10) are given by:

(1) The Dirichlet boundary condition:

(2) The Neumann boundary condition:

(3) The unified unilateral conditions on the third-type boundaries including seepage face Γs,rain infiltration surface Γi,and evaporation surface Γe(Borsi et al.,2006; Hu et al.,2017):

Fig.6.Location and landform of the Baihetan dam site.

where δ is a sign indicator,and h* and q*nare the maximum allowable values of pressure head and flux,respectively.The values of δ,h*and q*non the boundaries are listed in Table 2.The unilateral conditions on the third-type boundaries give rise to strong nonlinearity in Eq.(10),which can be efficiently solved with the discretized parabolic variational inequality(PVI)method(Hu et al.,2017).The algorithm was implemented in the FE code THYME(Chen et al.,2009).

Fig.7.Mean monthly precipitation and evaporation from the land surface at the dam site.

3.2.Inverse modeling technique

The saturated-unsaturated flow model (Eq.(10) and (11)) contains two constitutive relations (i.e.the water retention curve,Eq.(1),and the relative permeability function,Eq.(2)) and a set of material parameters(i.e.K or its isotropic counterpart K,Ss,θs,θr,α and n).The magnitude of hydraulic conductivity (K) for fractured rocks can generally be reliably determined by field hydraulic testing(e.g.packer tests),with rich experiences having been accumulated in engineering practices.The specific storage(Ss)can be estimated by hydraulic testing,analogy or the following equation:

where βrand βware the bulk compressibilities of fractured rock and water,respectively; φ is the porosity or the volume fraction of connected fractures and pores;ρwis the density of water;and g is the gravitational acceleration.

The saturated water content(θs)theoretically takes θs=φ when the effect of air trapping is negligible.The residual water content(θr)depends on the size of voids and pore throats,and is considered to have a less impact on the seepage field in fractured media (Liu et al.,2003).In the most simplified case,one may take θr≈0.

It is technically impossible to directly determine the VG parameters(α and n)by field tests for fractured rocks of large REV size.Alternatively,a sufficiently narrow range of α and n could be first predicted by Eqs.(8)and(9)from the more easily-obtained value of K.Then,the values of α and n could be much easily determined by inverse modeling,preferably using the values predicted by Eqs.(4)and (5) as the initial guesses.The inverse modeling utilizes the groundwater observations that can be reliably collected in saturated zones,and it works as the moisture movement in unsaturated zones would eventually influence the distribution of groundwater in saturated zones.

The objective function could be defined as (Zhou et al.,2015;Chen et al.,2020):

Fig.8.Groundwater flow system and observed groundwater levels.

where P={α1,n1,α2,n2,… αr,nr}Tis a vector denoting the VG parameters of materials to be back-calculated,r is the number of materials,M is the number of groundwater observation boreholes,N is the number of weirs,and Hiare the measured and simulated time series of groundwater level in borehole i,and Qjare the measured and simulated time series of discharge at weir j,and w is a weight coefficient introduced to ensure a balance between the relative errors of groundwater level and discharge measurements.

Various techniques are available to seek the optimal set of parameters by minimizing the objective function(e.g.neural network,genetic algorithm,Bayesian optimization).For large-scale inverse problems under field conditions,the procedure combining orthogonal design,finite element modeling of flow,artificial neural network and genetic algorithm generally works well,with a good balance between computational efficiency and accuracy (Zhou et al.,2015; Chen et al.,2016,2021; Hong et al.,2017).Interested readers may refer to Chen et al.(2020)for details of the algorithm.

4.Application

4.1.Site description

Located on the lower Jinsha River in Southwest China,the Baihetan hydropower station is the second largest hydropower station in the world.The project mainly consists of a double curvature arch dam of 289 m high and two underground powerhouse cavern systems.The site is typical of a steep-sided V-shaped valley formed through erosion by the river,with three perennial streams (i.e.Baihetan,Dazhai and Haizi) on the right bank,as shown in Fig.6.The site is situated in the subtropical monsoon climate zone,where 80%of the annual precipitation occurs during the wet season from May to October.Fig.7 plots the mean monthly precipitation and evaporation at the site.

The bedrocks outcropping in the right bank at the dam site are mainly the basaltic formations of the Emei Mountain group (P2β),which can be divided in to 11 basalt flow layers(P2β1-P2β11)according to their eruption sequence(Fig.8).The rock layers dip gently towards the right bank,with an orientation of N30-50 E/SE ∠15-25.There is a tuff layer about 0.2-1.8 m thick developed on the topofeach rock layers of the basaltic formations except the layer P2β1.Each tuff layer contains a shear zone(hereafter referred to as tuff zone)of 5-40 cm thick.The basalt formations are overlain by sedimentary rocks of the Feixianguan group(T1f),mainly consisting of low-permeability sandstone,claystone and clayey siltstone.Loose deposits (Q4) are mainly distributed in riverbeds and on gentle hillslopes.

The main geological structures at the site consist of subvertical faults(e.g.F19,F20and f232),gently-sloping tuff zones developed in the tuff layers(e.g.C3-C10)and local shear zones developed within each basalt layer (e.g.RS331,RS336and RS621),as shown in Fig.8.Fault gouge layers of 1-5 cm thick are developed in each tuff zone,as a result of tectonic shear between basalt layers (Fig.A5).Consequently,the permeability of the tuff zones is highly anisotropic,with their in-plane permeability (K//) about 1-2 orders of magnitude higher than the permeability perpendicular to the structural planes (K⊥).

The hydraulic conductivity of fractured rocks,geological structures and loose deposits were characterized by a variety of hydraulic tests including packer tests,pumping tests,recovery tests,slug tests and in situ seepage tests.For the fractured basalt rocks,the hydraulic conductivity decreases with increasing depth,and can be well correlated to the degree of weathering,the density of fractures,and the level of in situ stress(Chen et al.,2018).In dam engineering,the permeability of rocks is commonly categorized into six classes based on the Lugeon value(q)or the magnitude of hydraulic conductivity(K)(in units of m/s):very high(K ≥10-2),high(10-2＞K ≥10-4),moderate (10-4＞K ≥10-6),low (10-6＞K ≥10-7),rather low(10-7＞K ≥10-8),and very low(K ＜10-8).According to this classification,the rocks downward from ground surface are classified into five zones of decreasing permeability for the basalt formations(i.e.moderate I,moderate II,low I,low II,and rather low)and four zones for the sedimentary formations(i.e.moderate,low,rather low and verylow).The representative values of the hydraulic parameters(K,Ss,θs,θr)in each rock zone are listed in Tables 3 and A2.Also listed inTable 3 are the ranges of the VG parameters(α and n)estimated by Eqs.(8)and(9).

Table 3 Hydraulic properties for rocks at the site.

The groundwater at the site is mainly recharged by precipitation and discharged to streams (e.g.the Dazhai stream,Fig.6) and the Jinsha River.The infiltrated water mostly stores and moves in the loose deposits(Q4),and only a limited quantity could further infiltrate downward to recharge the basalt formations(P2β),due to low permeability of the sedimentary formations (T1f) between the aquifers.Groundwater observations fromboreholes and exploratory adits in 2001-2011 showed that the groundwater level in the deposits on the Xiahongyan platform is relatively stable,typically fluctuating within the range of 2.2-7.4 m(see Fig.8).The groundwater in the basalt formations(P2β)moves mainly along the largescale faults and tuff zones,characteristic of slow vertical infiltration and multiple water tables(Fig.8).The change of groundwater level in boreholes betweenwet and dry seasons is mostly within 10-20 m,with larger changes occurring in near-bank boreholes.

4.2.Numerical setup

A two-dimensional (2D) model along section IX-IX (with its trace shown in Fig.6) was created for saturated-unsaturated flow modeling in the right-bank slope.As shown in Figs.6 and 9,the left boundary of the model was aligned with the center of the Jinsha River,while the right boundary was located at the divide between the Dazhai and Haizi Streams.The dimension of the model was 4470 m long and 1680 m high,with 29,659 four-node isoparametric elements and 60,090 nodes (Fig.9).The topographic and geologic features were well represented in the mesh,where faults and tuff zones were also treated as a continuum of proper hydraulic properties (Table 3).The exploratory adits parallel to the profile (e.g.PD62) could not be accurately represented in the 2D mesh,and were approximately modeled with equivalent spacing per unit length(Chen et al.,2020).The 2D model could be applied for flow modeling because the groundwater at the dam site mainly flows along the direction perpendicular to the Jinsha River,especially in the near-bank slope where the flow behavior is the major concern in this study.

Fig.9.2D finite element mesh for the study area: (a) Overall view; and (b) Enlarged view for the near-bank region.The boundary conditions are also described in the plots.

According to the hydrogeological settings at the site,the boundary conditions of the model were specified as follows(Fig.9):(i)The riverbed surface was prescribed with a water head boundary according to the fluctuation of river stage,and a water depth of about 1 m was prescribed on the bed of Dazhai stream;(ii)The left and right vertical boundaries were assumed to be water divides,and were imposed with no-flow boundary.The base of the model was also assumed to be impermeable; (iii) The surfaces of the exploratory adits were taken as a seepage face; and (iv) The slope surface was prescribed as a precipitation/evaporation boundary with hp=0 m and hd=-200 m.Precipitation was assumed to occur in the middle ten days of a month,and the mean monthly precipitation and evaporation rates were plotted in Fig.7.

The average groundwater level observations in boreholes were used to calibrate a steady-state flow model,hence the initial head distribution in the saturated zone of the slope was estimated(Fig.A6).The head distribution in the unsaturated zone was determined by assuming a vertically linear reduction of saturation from 100%at the water table to 60%at the ground surface.On this basis,a natural precipitation/evaporation process was modeled to eliminate the influence of initial head distribution (Chen et al.,2020),which lasted for 3 decades (from May 1972 to May 2002).The site exploration process was then modeled(from May 2002 to December 2009).The initial time step and maximum time step were taken as 7 h and 7 d,respectively.

The time series of groundwater level observations in 5 boreholes drilled in the near-bank basaltic formations and 3 boreholes drilled in loose deposits on the Xiahongyan plat form,in addition to the discharge data collected at the 350 m near-face section of adit PD62(Fig.8),were chosen to construct the objective function (Eq.(17)).The weight coefficient(w)in Eq.(17)was calibrated to be 0.1 in this study.To reduce the difficulty in inverse modeling,some rock zones of comparable permeability were combined,resulting in a number of 10 kinds of materials and a total of 20 parameters for inverse modeling(Table 4).

During inverse modeling,an orthogonal table L161(920)was first constructed via orthogonal design by taking 9 levels for each parameter from the estimated ranges,producing a set of 161 parameter combinations.These parameter combinations were then used as inputs for saturated-unsaturated flow modeling.The calculated time-series data of groundwater level and discharge at the measuring points,together with the parameter combinations and the time-dependent precipitation/evaporation boundary condition,were used to construct and train a four-layer back propagation neural network(BPNN)model,with a structure of 21-23-21-9 neurons from the input to output layer.An optimal estimate of the VG parameters was finally obtained via genetic algorithm by minimizing the objective function (Eq.(17)),in which the size of initial population was 100,and the probabilities for crossover and mutation operations were 0.9 and 0.1,respectively.

4.3.Results

The inversed values of the VG parameters in different rock zones are listed in Table 4,and the corresponding water retention curves and relative permeability functions are plotted in Fig.10.The estimated parameters are comparable to the measurements by Perkinset al.(2014) for sedimentary rocks and the data of tuff rocks and fractures at Yucca Mountain (BSC,2000).Fig.11 compares the training error of the BPNN model with that when the ranges of VG parameters are empirically determined,as done in Chen et al.(2020).This plot evidences the necessity of using the upperbound estimates (Eq.(8) and (9)) in neural network training and inverse modeling.

Fig.10.The estimated (a) Water retention curves and (b) Relative permeability functions for different rock zones.

Table 4 The estimated VG parameters for rocks at the site via inverse modeling.

Fig.12 compares the observed and calculated groundwater levels(or depths)at five boreholes in basaltic formations and three boreholes in loose deposits.It shows that even there are some obvious discrepancies between the observations and the numerical simulations,the numerical results roughly reproduce the trend and magnitude of groundwater level observations with engineering accuracy.The distributions of water pressure head (in long-term natural state) in wet and dry seasons are plotted in Fig.13.There are two major water tables in the slope,an upper one in the loose deposits (Q4) and a lower one in the basaltic formations (P2β),as separated by the aquitard (T1f) between them.Furthermore,the high anisotropy of permeability in the subvertical faults (e.g.F19,F20) and gently-sloping tuff zones (e.g.C3-C10) makes the groundwater store and move mainly in this backbone structure system.As an example,Fig.14 shows the preferential flow through the backbone structures in the near-bank slope during the second half of May 2002.The simulations show that the mean infiltration rate of rainfall into the loose deposits is 28%,and only 7% of the precipitation could be further infiltrated to recharge the underlying basaltic formations in the local area of Xiahongyan platform(Fig.8).All the above results are fairly consistent with the site observations.

Fig.11.Training error of the neural network model.

Fig.15 plots the distribution of pressure head in December 2009,two years after the exploratory adit PD62 and its four branches PD62-1-PD62-4 were excavated.The plot shows that the excavation of the adits led to a significant drop of groundwater in the basaltic formations,but the groundwater level in the overlying loose deposits was less influenced.The discharge out of the adit system,measured with a weir at the 350 m near-face section of adit PD62,is shown in Fig.16.The measurement generally decreased with time and was rarely influenced by rainfall,but there was a sudden increase in October 2006 when the branches PD62-2-PD62-4 were excavated.The discharge then decreased to a relatively stable magnitude around 2 L/min.The simulated results roughly capture this feature.

4.4.Discussion

This study is motivated by the difficulty in the suction measurement in fractured rocks,which renders the saturatedunsaturated flow modeling extremely difficult or even impossible.The proposed upscaling relationships (Eq.(4),(5),(8) and (9))provide an initial guess and the upper bound estimation for the VG parameters of fractured rocks.It enables their representative values to be more easily and efficiently determined via inverse modeling using groundwater observations in saturated zones.

However,it should be noted that there are some limitations in the proposed method,such as the use of continuum approach,the negligence of hydraulic hysteresis (Hu et al.,2013) and flow patterns(Chen et al.,2017;Hu et al.,2019),the limited data availability for fractured media,and the use of 2D mesh,etc.A sensitivity analysis shows that the differences between the simulations and measurements could be moderately reduced but not minimized by a refined mesh (Table A3 and Fig.A7),indicating that the differences (Figs.12 and 16) are mostly attributed to the above factors.For example,the complex flow patterns (e.g.capillary fingering,viscous fingering,compact displacement and channeling) in fractured media governed by capillary number,flow rate and flow geometry (Chen et al.,2017; Hu et al.,2019) could not be properly modeled by continuum approach,but can be roughly described in an average sense.

Fig.12.Comparisons of the observed and calculated groundwater levels(or depths)at observation boreholes in the slope.The rainfall intensity is also plotted in each figure:(a) 2001-2003; (b) 2005-2007; and (c) 2007-2009.

Fig.13.Pressure head contours (a) in wet season and (b) in dry season.

Fig.14.Groundwater movement through the main structures consisting of faults and tuff zones in the near-bank slope during the second half of May 2002 at different time (a)t=0 d; (b) t=7 d and (c) t=14 d.

Fig.15.Pressure head contours in December 2009,i.e.two years after the exploratory adits were excavated.

Despite of the above limitations,the case study demonstrates that the proposed method is able to reproduce macroscopically most of the important unsaturated flow behaviors in fractured rocks,and the numerical results are accurate enough for engineering applications such as the design of grout curtain and the estimate of tunnel discharge.On the basis of continuum approach and the VG model,the proposed method is simple and practical,and can be readily applied to addressing large-scale saturated-unsaturated flow problems under field conditions.

5.Conclusions

The study is devoted to determining the representative values of the VG parameters for fractured rocks via inverse modeling using the groundwater observations in saturated zones,thereby overcoming the difficulties in saturated-unsaturated flow modeling through fractured rocks.To improve the efficiency and accuracy of inverse modeling,the upper bounds of the VG parameters are estimated via scaling relationships between the VG parameters (α and n)and hydraulic conductivity(K)established from a dataset of 1416 soil samples.The applicability of the scaling relationships to fractured media is then validated against a limited set of available data on 115 rock cores,fractures,and fractured formations.A case study of the saturated-unsaturated flow in the right bank slope at the Baihetan dam site demonstrates that the unsaturated hydraulic properties of fractured rocks on the field scale could be effectively determined via the proposed approach.Moreover,most of the flow behaviors under complex hydrogeological conditions could be reasonably reproduced via continuum approach-based numerical modeling.The proposed method provides a simple and practical technique for unsaturated flow modeling and tunnel discharge prediction in fractured formations when groundwater observations are available.

Fig.16.Variation of the measured and predicted discharges out of the adit system during the exploration period.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The financial supports from the National Natural Science Foundation of China(Grant Nos.51925906 and 51988101) and the National Key R&D Program of China (Grant No.2018YFC0407001)are gratefully acknowledged.The authors gratefully thank POWERCHINA HUADONG Engineering Corporation Limited for the field data supporting this study.

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2022.02.008.

List of symbols

SeEffective saturation

SsSpecific storage

θ Volumetric water content

θs,θrSaturated water content and residual water content

s Suction

K,K Saturated hydraulic conductivity tensor and its isotropic counterpart

K⊥,K//Vertical and horizontal components of hydraulic conductivity

α,n,m Fitting parameters of van Genuchten model

σα,σnStandard deviations of α values and n values

αLB,nLBLower bounds of α values and n values

αUB,nUBUpper bounds of α values and n values

krRelative permeability

λ,? Fitting coefficients

CwSpecific moisture capacity

t Time

v Darcy velocity

n Outward unit normal vector to the boundary

z Vertical coordinate

Γh,ΓqWater head boundary and flux boundary

Γs,ΓI,ΓeBoundaries of seepage face,rain infiltration surface,and evaporation surface

h Pressure head

h0Initial pressure head

ω Indicator function of h

δ Sign indicator

h*,q*nMaximum allowable values of pressure head and flux

hpPonding depth of water on the infiltration boundary Γi

hdMinimum pressure head on the evaporation boundary Γe

I(t) Rainfall intensity

ReMaximum evaporation rate in the field

βr,βwBulk compressibilities of fractured rock and water

φ Porosity

ρwDensity of water

g Gravitational acceleration

P Vector denoting the VG parameters of materials to be back-calculated

r Number of materials

M,N Numbers of groundwater observation boreholes and weirs

w Weight coefficient