Groundwater flow through fractured rocks and seepage control in geotechnical engineering: Theories and practices

Chuang-Bing Zhou,Yi-Feng Chen,Ran Hu,Zhibing Yang

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 Infrastructure Engineering,Nanchang University,Nanchang,330031,China

Keywords: Fractured rock Groundwater flow Flow visualization Hydraulic property Hydromechanical coupling Groundwater flow modeling Seepage control

ABSTRACT Groundwater flow through fractured rocks has been recognized as an important issue in many geotechnical engineering practices.Several key aspects of fundamental mechanisms,numerical modeling and engineering applications of flow in fractured rocks are discussed.First,the microscopic mechanisms of fluid flow in fractured rocks,especially under the complex conditions of non-Darcian flow,multiphase flow,rock dissolution,and particle transport,have been revealed through a combined effort of visualized experiments and theoretical analysis.Then,laboratory and field methods of characterizing hydraulic properties (e.g.intrinsic permeability,inertial permeability,and unsaturated flow parameters) of fractured rocks in different flow regimes have been proposed.Subsequently,highperformance numerical simulation approaches for large-scale modeling of groundwater flow in fractured rocks and aquifers have been developed.Numerical procedures for optimization design of seepage control systems in various settings have also been proposed.Mechanisms of coupled hydro-mechanical processes and control of flow-induced deformation have been discussed.Finally,three case studies are presented to illustrate the applications of the improved theoretical understanding,characterization methods,modeling approaches,and seepage and deformation control strategies to geotechnical engineering projects.

1.Introduction

Fluid flow through the subsurface rocks is an interdisciplinary topic closely related to a number of scientific and engineering fields including hydrogeology,rock mechanics,geotechnical engineering,and earth resources engineering(Berkowitz,2002;Neuman,2005;Meakin and Tartakovsky,2009;Chen et al.,2014a;Lei et al.,2017).Fractures can form and change in response to many factors such as lithostatic,tectonic,and thermal stresses,high fluid pressures,and fluid erosion (Brown,1995;Chen et al.,2018a;Frash et al.,2019;McBeck et al.,2019;Yang et al.,2019a).Thus,fractures are ubiquitous in subsurface rocks,occurring at different spatial scales,from microscopic to continental(National Research Council,1996).Groundwater flow through fractured rocks plays an important role in many geological hazards and engineering failures,ranging from landslides to damage in foundation structures,and to water inrush in tunnels and underground caverns (Li et al.,2014;Zhou et al.,2015a,2018).A fundamental understanding of mechanisms of groundwater flow in fractured rocks across scales is a prerequisite for predicting and controlling of seepage as well as ensuring engineering safety in many related geotechnical applications.

The study of groundwater flow through fractured rocks faces several challenges.First,fractured rocks are often characterized by multi-scale heterogeneity due to the geometrical features that exhibit variations in surface roughness,fracture apertures,fracture density and connectivity (Neuman,2005;Pyrak-Nolte and Nolte,2016;Yang et al.,2016).The hydraulic and geomechanical properties of fractured rocks are strongly dependent on these spatially varying geometrical parameters (Barton et al.,1985;Olsson and Barton,2001;Rutqvist and Stephansson,2003;Hyman et al.,2016;da Silva et al.,2019;Hu et al.,2019b).Second,the hydraulic properties of fractured rocks and the seepage field can undergo complex,long-term evolutionary changes due to disturbance caused by engineering activities during construction and operation such as excavation,grouting,and reservoir impoundment (Chen et al.,2015c,2021;Liu et al.,2016b).Third,predicative modeling of groundwater flow in fractured rocks at the field scale can be challenging because of the difficulties,e.g.in obtaining accurate hydraulic parameters of fractured rocks and in high-fidelity representation of engineering and geological structures given the computational resources restrictions (Chen,2022).

To address the above challenges and to better understand the mechanisms of groundwater flow in fractured rocks,many theoretical,experimental and simulation studies have been conducted(Zimmerman and Yeo,2000;Bonnet et al.,2001;Meakin and Tartakovsky,2009;Rutqvist and Tsang,2012;Hyman et al.,2016;Antonellini et al.,2017;Maldaner et al.,2019;Babadagli,2020).Laboratory techniques and field methods have been proposed for better characterization of hydraulic properties of fractured rocks(Rutqvist,2015;Massiot et al.,2017;Arshadi et al.,2018;Maldaner et al.,2019;Becker et al.,2020).Various numerical models and simulators have been developed and applied to studying groundwater flow processes in engineering practices,e.g.hydropower projects,geothermal exploitation,geological storage of CO2,and nuclear waste disposal (Rutqvist et al.,2002;Strack,2003;Chen et al.,2009;Taron et al.,2009;Hyman et al.,2015;Yang et al.,2015;Lei et al.,2017;Birkholzer et al.,2019;Hu and Rutqvist,2020;Toller,2022).These advances have deepened our understanding of groundwater flow in fractured rocks.However,there remain some key issues to be addressed:

(1) How does the flow regime change in response to factors such as fracture geometry,hydrodynamic condition,multiple fluid phases,solid phase dissolution,and presence of suspended particles?What are their microscopic mechanisms?How can we visually observe the flow behaviors?

(2) What theoretical models best characterize the hydraulic properties (especially for non-Darcian and unsaturated flows) of fractured rocks at different scales and their evolution with time? How can the inertial permeability be accurately determined? How to reliably interpret field data in various hydraulic tests under complex geological and engineering conditions?

(3) What are the efficient approaches to numerically solve the transient,saturated-unsaturated flow problems at the field scale? How to perform efficient inverse modeling? How to treat large conduits or underground rivers in modeling groundwater flow in karst aquifers?

(4) What are the efficient procedures to optimize the design of seepage control systems in geotechnical engineering?

(5) What are the mechanisms of groundwater flow-induced deformation? What are the impacts of reservoir impoundment on large-scale deformation and how can the deformation be mitigated with engineering measures?

To address the above issues,this paper first summarizes different types of groundwater flow models in geological media.Then microscopic mechanisms of fluid flow through fractured media are investigated,and theories and methodologies for determining hydraulic properties of fractured rocks are proposed.Approaches of field-scale groundwater flow modeling and procedures of optimization design for seepage control systems are developed.Mechanisms of groundwater flow-induced deformation and its control are also discussed.Finally,three cases studies of engineering application of the theories and models are demonstrated.

2.Flow models in geological media

Geological media contain a complex void system ranging from pores and fissures to fractures and conduits of distinct geneses,scales and geometries,and are hence characteristic of inherent heterogeneity and complex flow geometry.According to their pore structures,geological media are commonly classified into three categories: porous media,fractured media and karst media.From porous to fractured to karst media,the pore structure becomes more heterogeneous and the flow behavior becomes more complex.A variety of groundwater flow models are established to describe the different flow behaviors through geological media(Chen,2022),as shown in Fig.1.The continuum model applies in most cases to the groundwater flow through porous media,with sufficient accuracy and high computational efficiency.For fractured media,the fracture network model and dual porosity model are often used to describe the preferential/channeling flow in the fracture system and the interaction between the pore and fracture systems.For karst media,the dual and triple porosity models are preferable for characterizing the base flow through porous matrix,intermediate flow through fractures and quick flow through karst conduits(Shoemaker et al.,2008);but when only the groundwater balance and the recharge-discharge relationship are of concern,the lumped-parameter model becomes the primary choice for its simplicity and less data requirements(Fleury et al.,2009;Schmidt et al.,2014;Zhou et al.,2021).

The linear Darcy’s law is frequently assumed for groundwater flow through geological media.But with the increase of flow velocity,the relationship between flow velocity and hydraulic gradient deviates gradually from linearity,and the flow regime transitions from laminar to non-Darcian.An inertial term should be added to Darcy’s law to account for the nonlinearity of flow,which gives the Forchheimer’s law.The groundwater is a linear viscous fluid,but when grouting process is considered,the fluid becomes non-Newtonian.When two and multiple phases of miscible or immiscible fluids are involved,the flow patterns become much more complicated,as influenced by various factors including flow geometry,wettability,flow rate,etc.The unsaturated flow is a simplification of air-water two-phase flow,which describes the moisture movement in vadose zones and follows the Darcy-Buckingham’s law.Furthermore,the fluid flow in geological media is often impacted by and coupled with other physical and chemical processes,such as deformation,heat transfer,particle transport and chemical dissolution (Rutqvist et al.,2002;Chen et al.,2009).

The selection of a proper flow model requires that the site conditions be well represented and the problem of concern be addressed with expected accuracy.For example,one may find that the fracture network model is best suited for flow simulation when the fractures are sparse and the rock matrix is of extremely low porosity and low permeability.But with the increases of fracture density and porosity,the continuum approach gradually becomes a proper choice.It should be noted that other factors,such as data availability,computational cost and the modeler’s preference,could also impact the choice of flow model.Among the flow models shown in Fig.1,the saturated(single-phase)steady-state/transient flow model and saturated-unsaturated flow model are most commonly used to address the groundwater flow problems in geotechnical engineering,as listed in Table 1.It can be inferred from Table 1 that from steady-state to transient to saturated-unsaturated flow,and from Darcian to non-Darcian flow,the nonlinearity of the governing equations and the number of hydraulic properties that needs to be determined increase,hence increasing the difficulties in solving the partially differential equations and in determining the hydraulic properties of the media.

Table 1 Groundwater flow models commonly used in geotechnical engineering.

Fig.1.Classification of groundwater flow models (modified from Chen,2022).

3.Microscopic mechanisms of fluid flow through fractured media

Understanding the microscopic mechanisms is fundamental to developing constitutive models for hydraulic properties of fractured rocks.Flow-visualization experiment is an essential methodology to identify and quantify the microscopic mechanisms.This section starts with the flow-visualization experiments and then presents the observation and quantification of flow behaviors,including non-Darcian flow,two-phase flow,fracture dissolution,and particle transport.

3.1.Visualization experiments

Flow-visualization technologies are an effective method that has been employed to investigate microscopic mechanisms of complex flow in fractured media.These visualization technologies include optical imaging (Detwiler et al.,1999;Hu et al.,2017b),computed tomography (CT) imaging (Karpyn et al.,2007),and nuclear magnetic resonance (NMR) method (Dijk et al.,2002).Compared to the CT imaging and NMR method,optical imaging can permit real-time imaging of flow behavior both at the microscopic and the sample scales with relatively low cost,and thereby has become a widely-used technology to study mechanisms of fluid flow in fractured media.Fig.2a shows the flow-visualization experimental setup based on optical imaging technologies (Chen et al.,2017;Hu et al.,2017b,2018a,2019a;Wu et al.,2021a;Wang et al.,2022),consisting of a flow-rate control system,an imaging system,a transparent fracture flow cell,and a temperature-control system.

In the flow-rate control system,two syringe pumps are used to inject dyed liquid into the fracture with precise flow rates from 0.001 mL/min to 100 mL/min.In the imaging system,an inverted microscope is installed below the fracture to record images of fluid flow at the micro-scale,and another two 14-bit charge-coupled device (CCD) cameras are installed over the fracture that records images at the sample scale.A micro-PIV (particle image velocimetry) system is also included in the imaging system to capture the flow velocity fields.Pressure sensors and a concentration sensor are connected to the fracture,providing a real-time recording of the pressure difference and concentration at a sampling frequency of 100 Hz.The entire experimental system is contained in an insulated enclosure to maintain a temperature ranging from 20°C to 45°C.

One of the challenges in the flow-visualization experiments is the fabrication of fracture.Different from the fracture samples used in the CT and NMR experiments,the fracture sample in the optical imaging experiments should be transparent and well represent the characteristics of natural rough surfaces.To meet this condition,the fracture is fabricated using the replica cast method through the following five steps (Chen et al.,2017): (1) creating a naturalfracture surface,(2) creating silicone rubber mould,(3) molding fracture replicas from the mould,(4) assembling the two replicas,and (5) modifying surface wettability as desired.Following this procedure,the transparent fracture can well represent the aperture field and the rough surfaces of the natural rock fractures (Fig.2b).Through the visualization system,the fluid flow processes through rough fractures can be visualized in real-time,including non-Darcian flow (Fig.2c),two-phase flow through a single fracture(Fig.2d) and a fracture network (Fig.2e),fracture dissolution(Fig.2f),and fluid-driven particle transport (Fig.2g).

3.2.Non-Darcian flow through rough fractures

Non-Darcian flow behaviors through rough rock fractures have been widely observed in both natural groundwater flow systems and forced flows (such as hydraulic fracturing,pumping or injection wells,dams,and tunnels).This nonlinear flow characteristic has been ascribed to the formation and growth of eddies caused by significant inertial forces.Accurate understanding of eddy evolution is of vital importance to characterizing non-Darcian flow in rock fractures.This section introduces quantitative characterization of eddy growth and its evolution and discusses the mechanism controlling the onset of non-Darcian flow.

Fig.2.(a)Schematic of the flow-visualization system for fluid flow through rough fractures,(b)The transparent rough fracture,(c-g)Observation of fluid flow processes for(c)Non-Darcian flow,(d) Two-phase flow,(e) Unsaturated flow in a fracture network,(f) Fracture dissolution,and (g) Particle transport and deposition.

3.2.1.Eddy development and growth

The occurrence of eddy greatly narrows the main flow channel,which is responsible for the occurrence of non-Darcian flow behaviors.The development and growth of eddies in rock fractures are shaped by both geometric configurations and hydrodynamic conditions.Cardenas et al.(2009) made a comparative investigation on the flow velocity field of natural rock fracture obtained by solving Navier-Stokes equations and Stokes equations.It was found that the eddy size obtained by solving Navier-Stokes equations increased with the Reynolds number.However,for the result of the Stokes equation without considering the inertial term,the eddy did not change with the Reynolds number and seemed to be determined only by the fracture geometry.Lee et al.(2014) directly observed the eddy growth with increasing Reynolds number in rock fracture with the help of a high-speed camera and fluorescent particle imaging technology.In their experiment,the correlation between the eddy and the inertial term of Navier-Stokes equation was confirmed.To quantitatively characterize the eddy in rough fractures,an eddy detection technique by Zhou et al.(2019a) was adopted to analyze the solved velocity fields.This technique was developed based on a zero-integral-flux method,which can precisely detect the eddy boundary and further calculate the total eddy volume.Eddies near the upper and lower walls of regions with large aperture expansions were evident visually(Fig.3a),showing that there was a large difference in the velocity magnitude between the eddy and the main channel(Fig.3b).With increasing Reynolds number,the eddy initially developed as viscous Moffatt eddies(Moffatt,1964;Cardenas et al.,2009),and then grew because of inertial effects(Fig.3c).Besides,zone shifting and sudden change in eddy area can be observed with increasing Reynolds number(Fig.3c).These diverse eddy shapes and complex growth patterns were controlled by geometric and hydraulic factors under the joint influence of viscous and inertial forces.

3.2.2.Onset of non-Darcianflow

Non-Darcian flow behaviors in rock fractures have been widely observed in direct numerical simulations (Cardenas et al.,2009;Zhou et al.,2019a),laboratory experiments(Javadi et al.,2014;Chen et al.,2015d),and in situ tests (Quinn et al.,2011;Chen et al.,2015b).The critical Reynolds number (Rec) has been widely used for characterizing the onset of flow transition to non-Darcian,and a wide range ofRecranging from 0.001 to 2300 has been suggested for rock fractures,as summarized by Zhou et al.(2015b).This larger difference inReccan be attributed to the non-unique criteria of flow regime transition,inconsistent definitions ofRe,and variability of fracture geometries and flow channels.This makes the prediction of non-Darcian flow onset inconvenient for practical applications.Zeng and Grigg (2006) recommended a Forchheimer number representing the ratio of pressure drop caused by liquid-solid interactions to that caused by viscous resistance as a metric for characterizing non-Darcian flow.This strategy was implemented by Javadi et al.(2014),who developed a quantitative criterion that combines the Forchheimer’s law (-?P=AQ+BQ2) andReto quantify the onset of non-Darcian flow in rock fractures.This criterion quantitatively defined the critical point where the role of nonlinear pressure drop(BQ2)reaches to α percentage of the overall pressure drop(AQ+BQ2),as illustrated in Fig.4a α=10%has been widely accepted to determine this critical point (Zhang and Nemcik,2013;Javadi et al.,2014;Zhou et al.,2015b,2019b).Based on a series of laboratory experiments,Zhou et al.(2015a)proposed an empirical power-law function that relatesRecto the hydraulic aperturebhand α,as shown in Fig.4b.This proposed model has clear physical meanings,offers an explicit quantitative criterion for determiningRecand provides a method for quantitatively evaluating the non-Darcian flow effect.

3.3.Two-phase flow patterns through rough fractures

The displacement of one fluid by another immiscible one in rough fractures is an important process in many subsurface engineering applications.Different from the single-phase flow,the fluid-fluid interface involved in this two-phase flow exhibits various patterns.Although two-phase flow patterns through porous media have been investigated for more than 30 years(Lenormand et al.,1988),few studies focus on the flow patterns through rough fractures.Different from porous media,the fractures are inherently rough to various degrees,and the roughness induces irregular flow passages,which play a fundamental role in the formation and transition of the two-phase flow patterns.This section presents the experimental observation of two-phase flow patterns through rough fractures and then quantifies the roughness effects on the two-phase flow patterns.

Fig.3.Eddy development and growth in rock fractures:(a)A rock fracture with a series of eddies,(b)Velocity field at maximum Reynolds number with detected eddy boundaries(red lines),and (c) Eddy evolution with increasing Reynolds number.Modified from Zhou et al.(2019c).

Fig.4.Definition of the critical point of non-Darcian flow and the critical Reynolds number(CRN)model:(a)Quantitative criterion to quantify the onset of non-Darcian flow in rock fractures,and (b) CRN model where the onset of non-Darian flow can be predicted by the hydraulic aperture bh.Modified from Zhou et al.(2015b).

3.3.1.Observation of theflow patterns

Governed by the interplay between the capillary and viscous forces (e.g.Hu et al.,2018b),the instability of the fluid-fluid interface at the micro-scale forms different flow patterns ranging from capillary fingering,viscous fingering,and compact displacement(Yang et al.,2019b).This interplay can be quantified by two dimensionless numbers,i.e.the capillary numberCadefined asCa=μivi/σ and the viscosity ratioMdefined asM=μi/μd.Here μiand μdare,respectively,the viscosity of the invading and defending fluids;viis the Darcy velocity of the invading fluid;and σ is the interfacial tension.Based on the experimental observation,a phase diagram of two-phase flow patterns for a rough fracture is established(Chen et al.,2017,2018b),as shown in Fig.5.All of the typical flow patterns formed in porous media are also observed in the rough fracture.In the case ofM<1,the crossover zone from capillary fingering to viscous fingering is for the first time observed in the rough fracture,with lower displacement efficiency compared to that in the regimes of viscous and capillary fingering.The crossover zone in the rough fracture is much narrower than that in porous media (Lenormand et al.,1988),which is attributed to the different flow geometry,especially the fracture roughness.

In the case ofM>1,according to the two-phase flow theory,the flow pattern should be stable at relatively largeCa.However,the localized flow channel,a previously unidentified flow pattern,is observed in the rough fracture.The underlying mechanism for the onset of the localized flow channel in the rough fracture is that the heterogeneity of flow geometry,i.e.the roughness,triggers fingers under favorable conditions (M>1).Due to the rough walls of the fracture,the invading fluid did not fully occupy the fracture void space,i.e.incomplete displacement occurred (Chen et al.,2018b).This localized flow channel is observed in rough fractures but not in smooth-walled fractures,highlighting the destabilizing effect of the roughness on two-phase flow patterns through rough fractures.

3.3.2.Roughness effects on theflow patterns

Fig.5.Phase diagram of two-phase flow patterns in a rough fracture,showing the domains of flow patterns,including capillary fingering,viscous fingering,compact displacement,and the crossover and localized flow channel.The displacement efficiency (indicated by the color bar) as a function of capillary number Ca and viscous ratio M is also shown (Hu et al.,2018b).

Fig.5 implies the destabilizing effect of roughness on two-phase flow patterns.This subsection would establish a theoretical model to quantify this effect and then present a phase diagram of the flow patterns,including the roughness index and the capillary number.The fracture roughness index is defined as the coefficient of variation of the aperture field λb=σb/,where and σbare,respectively,the mean and the standard deviation of the aperture field(Glass et al.,2003).The experimental results showed that the two-phase flow exhibits nearly compact displacement for λb<0.05,indicating that the roughness effects can be neglected in this range(Fig.6a and e).In this case (λb<0.05),the rough-walled fractures can be simplified as the Hele-Shaw channel (parallel-plate).For λb>0.05,the fracture roughness affects the formation and transition of flow patterns,which is controlled by the competition between viscous and capillary forces(as discussed above).Under slow flow conditions (smallCa),the capillary force dominates the flow,resulting in wider fingers and a more stable interface.AsCagradually increases,the influence of viscous force becomes important,destabilizing the interface.AsCareaches a critical value at which the two effects of capillary and viscous forces are comparable,there exists the transition of flow patterns.On this basis,a theoretical model for the transitions of flow patterns was established(Hu et al.,2019b).Fig.6a shows that the proposed theoretical model well describes the two transitions of the flow patterns from capillary fingering (Fig.6d) to the crossover,denoted byCac,and from the crossover (Fig.6c) to viscous fingering (Fig.6b),denoted byCaf.It also demonstrates the significant effects of roughness in the capillary-fingering regime.In this regime,the irregular flow channels induced by the roughness would increase the fluctuation of capillary force that destabilizes the displacement front,and thusCacdecreases with λb.

Fig.6.(a)Phase diagram of the two-phase flow patterns through rock fractures(Hu et al.,2019).The experimental data are denoted by dots.The theoretical model is able to predict the boundary between capillary fingering(CF)and crossover zone(CZ),labeled with Cac,and the boundary between crossover zone(CZ)and viscous fingering(VF),labeled with Caf.The representative flow patterns (labeled with open circle) are shown in (b) viscous fingering,(c) crossover zone,(d) capillary fingering,and (e) nearly compact flow in smooth fracture.

Fig.7.(a) Evolution of the transient splitting ratio η with dimensionless time t* which is defined as the volume ratio of water passing the junction to the initial volume,and (b)Relationship between the cumulative splitting ratio η* and initial liquid slug length L (Xue et al.,2020).The dashed lines denote the critical liquid slug length.

3.4.Two-phase flow behavior through fracture networks

Based on the understanding of the two-phase flow patterns through rough fractures,this section presents the flow behaviors through fracture networks.Different from the single rough fracture,the fracture intersections control the flow splitting behavior and play an important role in the flow patterns through fracture networks.

3.4.1.Flow-splitting behavior in fracture intersections

The flow-splitting behavior in fracture intersections has a considerable influence on the macroscopic flow characteristics of the fracture network.Fracture intersections can act as capillary barriers,leading to unstable and intermittent flow behaviors(Wood et al.,2005;Yang et al.,2019c),which is the main reason for complex spatial and temporal dynamics in fracture networks.On the basis of observations of the droplet splitting experiments,a two-dimensional (2D) liquid dynamic splitting model at fracture intersections was proposed by combining the single fracture seepage theories and the solid-liquid-air three-phase contact line motion theory (Xue et al.,2020).The proposed model can completely describe the splitting process and accurately predict the splitting ratio of liquid at the fracture intersection,as shown in Fig.7a.It was found that there is a critical initial length that divides flow pattern into type I(dominated by flow in the branch channel)and type II(dominated by flow in the main channel)behaviors(Fig.7b),as suggested by both the model and the experimental results.Based on the mechanistic model,the effects of several key physical parameters such as fracture inclination,aperture,and dynamic contact angle on the liquid infiltration behavior were systematically studied.

3.4.2.Unsaturatedflow in discrete fracture networks

Field and laboratory experiments have indicated the complexity of temporal and spatial characteristics of unsaturated flow in fracture networks (Dahan et al.,1999).Through unsaturated flow experiments with a fracture network fabricated by threedimensional (3D) printing technology,liquid breakthrough path,breakthrough time,and spatial distribution of local flow states were investigated (Zhou et al.,2022).The experimental results showed obvious gravity-driven preferential flow pathways(Fig.8).Based on analysis of water droplet migration in a single fracture and capillary barrier effect-induced liquid accumulation at the fracture intersections,a theoretical model of liquid breakthrough time through fracture networks was proposed.The flow status in fracture networks can be classified into six categories: saturated section,film flow section,dry section,trapped air section,pulse section,and wet-surface section.The improved understanding and the proposed model provide new physical insights into the complex unsaturated flow dynamics and are of significance for predicting the time for water/contaminant to travel to a certain location in unsaturated fractured rocks.

3.5.Dissolution patterns through rough fractures

The dissolution of soluble minerals within a fracture expands the fracture aperture and produces distinct dissolution patterns,which eventually control the enhancement of hydraulic properties.Understanding and quantifying the dissolution patterns through rough fractures is fundamental to the seepage control for tunnels and dam foundations in karst areas (Romanov et al.,2003;Zheng et al.,2021).This section presents experimental observations and a theoretical model for the dissolution patterns and then reports a new dissolution phenomenon,i.e.dissolution hotspots.

3.5.1.Observation and transition of dissolution patterns

From the micro-scale perspective,the dissolution dynamics are controlled by the interplay among advection,diffusion and reaction.The dissolution dynamics at the sample scale exhibit three patterns: compact pattern,wormhole,and uniform pattern.The Péclet numberPeand the Damk?hler numberDaare often employed to quantify the interplay.Flow-through dissolution dissolution experiments in transparent radial fractures were performed (Fig.2),and three distinct dissolution patterns were observed.Two length scales(lfandlt)were employed to represent the geometry of dissolution channels,wherelfis the penetration length of the dissolution channel in the flow direction,andltis the penetration length in the transverse direction,as shown in Fig.9b.Iflfis smaller thanlt,the solid phase is mostly dissolved near the inlet,leading to a compact dissolution pattern.Conversely,whenltis much smaller thanlf,orlfis larger than the system size,the dissolution exhibits a uniform pattern.When the two length scales(lfandlt) are on the same order of magnitude,wormhole patterns form.Based on the analysis of penetration lengths of dissolution channels,a theoretical model of critical Damk?hler numbers (DacandDau)was developed(Fig.9a),and then a phase diagram in the 1/Da-Peplane was proposed for the dissolution pattern through the radial rough fracture (Wang et al.,2022).The phase diagram provides insights into the transition of dissolution patterns in rough fractures and is important in the management of reactive flow in geological systems.

3.5.2.Dissolution hotspots and the underlying mechanism

Fig.8.Spatial distribution of water at the steady-state(modified from Zhou et al.,2022).Water is injected at a single inlet(a-c)or multiple inlets(d-f)at a total flow rate Q,and the flow direction is from top to bottom: (a) Q=0.03 mL/min,(b) Q=1 mL/min,(c) Q=5 mL/min,(d) Q=0.1 mL/min,(e) Q=4 mL/min,and (f) Q=20 mL/min.

Fig.9.(a)Phase diagram of dissolution patterns in the 1/Da-Pe plane(Wang et al.,2022).The transition from compact to wormhole patterns are denoted by Dac,and the transition from wormhole to uniform pattern is denoted by Dau.(b)Two length scales(lf and lt)represent the shape of dissolution patterns.If lf is smaller than lt,compact dissolution pattern forms.If lt is much smaller than lf,or lf is larger than the system size,the dissolution exhibits a uniform pattern.The red color means a higher γ (= lf/lt),whereas the blue color represents a lower γ.The dots denote the experimental results.

Fig.10.Dissolution hotspots (modified from Hu et al.,2021):(a)Compact dissolution pattern,for which the aperture growth<Δb>r along radial profiles decreases with r(a′),(b)Uniform dissolution pattern where<Δb>r nearly keeps unchanged,(c)dissolution hotspot where<Δb>r first increases with radius r and then reaches its maximum at rhot(c′),and(d) Theoretical model for predicting the locations of dissolution hotspots.

Under radial flow conditions,the velocity decreases as the fluid flows away from the inlet (the center of the fracture),which significantly affects the increase of fracture aperture.Fig.10a′,b′and c′shows the radial profiles of aperture growth <Δb>rfor the three dissolution patterns.It can be seen that<Δb>rdecreases with radiusrin the compact dissolution pattern (Fig.10a and a′),but nearly keeps unchanged in the uniform dissolution pattern(Fig.10b and b′).Surprisingly,a previously unidentified phenomenon,dissolution hotspots,was observed in the wormhole dissolution pattern (Fig.10c).In this case,<Δb>rfirst increases with radiusrand then reaches its maximum at the locationrhot.The peak of the curve is the hotspot of dissolution because the local dissolution rate reaches the maximum (Fig.10c′).

The occurrence of the dissolution hotspot corresponds to the state where the effects of mass transport in the horizontal and vertical directions are comparable.The timescale in the horizontal directiontadvis determined by the advection of reactive fluid from the inlet to the location at radiusr,and the timescaletveris controlled by the buoyancy-driven instability in the vertical direction.The fracture dissolution is limited by the transport in the vertical direction iftadvtver.The dissolution hotspot occurs whentadv=tver.On this basis,a theoretical model for the dissolution hotspot was established (Hu et al.,2021).The model can well predict the locations of dissolution hotspots and reveal the underlying mechanism for the occurrence of hotspots (Fig.10d).The experimental observations and quantification of the dissolution hotspots highlight the key role of gravity in fracture dissolution,even in horizontal flows.

3.6.Fluid-driven particle transport patterns in fractures

Particle transport in fractures is of paramount importance for hydropower engineering and environmental restoration.Although considerable efforts from laboratory to field experiments have been devoted to the investigation of flow-driven particle migration in fractures,the coupled effects of various factors remain unclear,which hinders a full understanding of particle transport and retention in fracture networks.

3.6.1.Mechanism for particle retention

Previous studies (e.g.Wan and Tokunaga,1997) have revealed that the dominant mechanisms for particle retention include attachment at the solid-liquid interface,air-water interfacial capture,film straining,and pore exclusion.During unsaturated flow,deformation of the air-liquid interface induces an additional force on the particles,which can be decomposed into a lateral component that tends to push the particles toward the bulk liquid,and a vertical component that pushes the particles against the wall of the channels,as shown in Fig.11.The particle entrapment can only be observed when the total resistance (Fr) acting on the particles equals the sum of the lateral capillary force and hydrodynamic drag force (Fd) in the opposite direction (Li et al.,2022).Based on lubrication theory and force balance analysis,three regimes of particle retention were characterized,and a modified critical capillary number for particle retention under different conditions was obtained (Wu et al.,2021b).The results indicate that the particle retention largely depends on the flow rate,particle volume fraction,and geometry of the channel.This work elucidated the particle transport mechanism during unsaturated flow and therefore is of theoretical and practical importance to understanding the contaminant and sediment migration in many natural and engineered systems.

3.6.2.Particle deposition patterns

Particle deposition and fracture clogging depend on several factors,including the interaction forces,the fracture aperture-toparticle size ratio,the volume fraction of migrating particles,path tortuosity and ensuing inertial retardation (e.g.Valdes and Santamarina,2007).Through flow-visualization experiments,the evolution of particle deposition patterns through rough fractures with different geometries was observed.Due to particle deposition and clogging,the effective aperture and the permeability decrease.The distribution of clogging and the resulting aperture change are presented in Fig.12,showing the significant effect of the geometric ratioNgdefined as the ratio of the mechanical aperturebmto the median particle sized50.The probability of clogging reduces with increasingNg.

Fig.11.(a)Schematic of particle retention during air-suspension displacement;(b)Forces acting on a particle at the air-liquid interface,where Fc represents the capillary force,Fd is the hydrodynamic drag force,and Fr is the resistance due to,e.g.friction;and(c)Probability of particles retention as a function of volume fraction at different flow rates(Wu et al.,2021b).

Fig.12.Observation of deposition patterns in a set of fractures: (a-c) Evolution of deposition-induced aperture reduction Δb under various geometric ratios of Ng=11.1,24.2 and 37.4,where b0 is the initial local aperture;and (d) Variation of hydraulic conductivity during fracture deposition process,which can be best-fitted with a negative exponential model,where K0 denotes the initial hydraulic conductivity.

For a small geometric ratioNg=11.1 (Fig.12c),the particle deposition first occurs near the entrance and then propagates towards the outlet.During the fracture clogging process,the permeability is reduced by two orders of magnitude(Fig.12d).At a large geometric ratio (Ng=37.4),the particle deposition occurs throughout the fracture plane (Fig.12a),and the variation of the local aperture is affected by the combined effect of deposition and traction.The aperture reduction first increases more rapidly in small local aperture regions,and then particles are mainly deposited in large local aperture regions as time evolves and the permeability only reduces slightly (Fig.12d).For the intermediate geometric ratio(Ng=24.2),the local aperture reduction is affected by complex particle transport processes(Fig.12b),and permeability is generally reduced within one order of magnitude.These deposition patterns enrich our understanding of particle transport and fracture clogging in natural fractures,and provide a basis for further theoretical analysis and explanation for the time-dependent effect of permeability under field conditions (Section 4.3.3).

3.7.Upscaling of fluid flow from micro-scale to Darcy scale

Sections 3.2-3.6 discuss the complex flow patterns through rough fractures and fracture networks and explain the mechanisms that control the flow behaviors at the micro-scale.Upscaling methodologies are essential to determine the macroscopic-scale properties of the fluid flow from the experimental results.For single-phase flow,the methodologies of volume-averaging,thermodynamically constrained averaging (Hassanizadeh and Gray,1993;Battiato et al.,2019),and energy-based averaging (Zhou et al.,2018) have been widely used for porous and fractured rocks.However,developing upscaling methodologies for twophase flow in fractured media remains challenging because of the fluid-fluid interface involved.In the widely used volume-averaging method,it is ambiguous to upscale the pore-scale pressure to the Darcy-scale viscous pressure drop,mainly due to the capillary force at the fluid-fluid interface.Nevertheless,previous studies (e.g.Raeini et al.,2014) have shown that the Darcy-scale viscous pressure drops can be unambiguously defined using the energy dissipation rate.Therefore,it is of great potential to develop an upscaling methodology for two-phase flow in fractured media from the perspective of energy dissipation,and even it remains an open issue at present.

4.Hydraulic properties of fractured rocks

Characterizing the hydraulic properties of fractured rocks(Table 1) is the primary step to modeling of groundwater flow and design of impervious systems.A comprehensive methodology including laboratory tests,field tests,statistical correlations and inverse modeling is presented for determining or estimating the hydraulic properties of fractured rocks representative of site conditions.

4.1.Permeability of intact rocks

Compared to fractures,intact rocks are usually of low permeability and contribute less to the effective permeability of fractured rocks.However,the permeability of intact rocks offers a solid lower-bound estimate for the permeability of fractured rocks,and it may become significant when fractures are sparse or tightly closed,or rocks are intensively damaged.The transient pulse technique plays an important role in measuring the permeability of low-porosity rocks(Brace et al.,1968;Chen et al.,2014a;Liu et al.,2017a),which can be conveniently integrated into a triaxial cell for permeability measurements under confining pressure and various stress paths.The permeability of rocks is determined by the nature of pore and/or crack systems,and hence the permeability of low-porosity rocks typically evolves with damage growth (Chen et al.,2014a) and the cumulative acoustic emission counts (Li et al.,2020b) under conventional triaxial compression.Fig.13 shows the permeability variation of a granitic rock under increasing deviatoric stress,with a clear permeability decrease in the initial crack closure region,a constant permeability in the elastic region and a dramatic permeability increase in the crack growth region.The permeability may increase by 2 or even higher orders of magnitude as the deviatoric stress is increased up to sample failure.

This permeability variation pattern of intact rocks can be well formulated with damage models (Shao et al.,2005;Arson and Pereira,2013;Chen et al.,2014b) or simply with curve fitting models (Li et al.,2020b) by taking into account the microstructural alteration during loading.These models can then be applied to assessing the excavation-induced damage growth and permeability variation in the surrounding rocks of tunnels in massive rocks,as shown in Fig.14.It shows that the permeability of the massive rocks in the vicinity of the excavation surface could be enhanced by 2-4 orders of magnitude,depending on the lithology,geostress level,tunnel geometry and alignment.Therefore,this permeability enhancement is important for various applications such as radioactive waste disposal,underground oil/gas storage and tunnel sealing treatment.Furthermore,it is worthy of noting that the damage zone evaluated by permeability/transmissivity enhancement is always larger than that assessed by seismic velocities,microseismic events or deformation measurements.

Fig.13.Variation in permeability of a granitic rock with deviatoric stress at 10 MPa confining pressure(Chen et al.,2014a):(a)Stress-strain curve,and(b)Permeability-deviatoric stress curve.

Fig.14.Damage growth and permeability variation in the surrounding rocks of the TSX tunnel:(a)The damage zones assessed by the MVP(microvelocity probe)method and the SEPPI (in situ pulse test) method (Martino and Chandler,2004),(b) The damage density assessed by a micromechanics-based damage model (Chen et al.,2014b),and (c,d) The permeability variation along P2 and P3 directions.

4.2.Permeability of rough fractures

The permeability of single fractures is fundamental to the understanding of water conductivity in fractured rocks,which is influenced by various factors including the fracture aperture,surface roughness,infillings,and the lithology and weathering degree of rocks on both walls.The cubic law,originally proposed for the laminar flow of fluids through smooth-walled parallel plates,is modified to describe the fluid flow through rough fractures(Witherspoon et al.,1980;Barton et al.,1985;Zimmerman and Bodvarsson,1996).A hydraulic apertureehis commonly introduced to account for the deviation of rough morphologies from the ideal condition assumed in the parallel plate model.The hydraulic aperture simplifies the description of flow behavior through rough fractures and avoids the difficulty in aperture measurement under field conditions.The magnitude ofehis always smaller than the mean mechanical apertureeof fractures,and the correlation betweenehandedepends on the statistics of surface morphologies.The cubic law indicates that the permeability of rough fractures is proportional to.

The permeability of fractures varies when they undergo compaction or dilation.Under shear displacement,the fractures in hard rocks start to dilate as the shear stress approaches the peak strength,and the permeability will be significantly enhanced by about two orders of magnitude as shear deformation proceeds(Barton et al.,1985;Esaki et al.,1999),as shown in Fig.15a.The shear dilation and permeability enhancement of fractures are suppressed by increasing normal stress (Bandis et al.,1983).These phenomena could be well characterized with the concept of mobilized dilatancy (Chen et al.,2007).Another scenario of permeability change occurs when sediments migrate and deposit in fractures,which causes clogging of fracture aperture and decrease in permeability (Chen et al.,2021),as shown in Fig.15b.Furthermore,under increasing pressure gradient,the flow regime in fractures first transitions from laminar to non-Darcian,and then the fractures start to dilate (Zhou et al.,2019a).Correspondingly,the cubic law predicts a decreasing apparent permeability in the non-Darcian flow regime and then an enhancing permeability with further increasing pressure gradient (Fig.15c).

4.3.Permeability of fractured rocks

Given the heterogeneous nature of fractured rocks,the permeability is of high spatial variability and varies as the rocks undergo deformation,damage,fracturing and clogging.This section presents the spatial distribution of the permeability of fractured rocks in deep-incised valleys,and the patterns of permeability changes induced by excavation and fracture clogging.

4.3.1.Statistical distribution of permeability of fractured rocks in deep-incised valleys

Over the last three decades,a considerable number of high dams have been built in southwest China,where the river valleys are deep-incised and the bedrocks of various types outcrop as a result of fast river incision and active tectonics.The bedrocks have experienced distinct degrees of weathering and unloading downwards from the ground surface,which leads to stronger spatial variability of rock permeability compared to fractured aquifers in other areas.Numerous packer tests were performed at the dam sites to characterize the permeability of foundation rocks and for the design of impervious barriers.A permeability database RocKDB was developed,which up to now contains the data of 21,081 packer tests performed at 18 high dam sites in southwest China.The database records the depth,length,lithology,structural feature,weathering degree,unloading degree,P-Qcurve,Lugeon value (q),and hydraulic conductivity (K) of each test interval.

Fig.15.Permeability variation in fractures induced by (a) shear deformation (Esaki et al.,1999;Chen et al.,2007),(b) sediment deposition (see Section 3.6),and (c)increasing pressure gradient (Chen,2022).

Fig.16 shows a subset of the data containing 13,397 pack tests performed in 614 boreholes at 12 high dam sites (Chen et al.,2018a).The variability and spatial distribution of hydraulic conductivity for fractured rocks in deep-incised valleys are summarized as follows:

(1) The hydraulic conductivity (K) follows the log-normal distribution,as verified by the normal quantile plot and the Mann-WhitneyUtest(Fig.16a).Furthermore,the log-normal distribution also approximately describes the distribution ofKat each depth within 600 m from ground surface and the distribution ofKat each site.

Fig.16.Distribution of hydraulic conductivity(K)data for fractured rocks in deep-incised valleys(Chen et al.,2018a):(a)Histogram and normal probability density function(PDF)of log-transformed K values,(b) Variation of fractured rock’s K with depth,(c) Variation in fault’s K with depth,and (d) Variation of standard deviation for log-transformed K data(σlog10 K) with depth. N denotes the number of tests.

(2) The vertical distribution of hydraulic conductivity follows a power-law correlation with depth,which is valid not only for fractured rocks (Fig.16b),but also for structural planes such as faults and shear zones (Fig.16c).The extensional faults overall display higher hydraulic conductivity,while the compressional faults and shear faults have comparable magnitude of hydraulic conductivity.The power-law correlation comprehensively reflects the influences of fracture density,weathering,unloading and geostress level on the hydraulic conductivity of fractured rocks,which reads

whereKsis the hydraulic conductivity near the ground surface(m/s),dis the vertical depth from the valley flanks (m),and α is an exponent.

(3) The standard deviation of log-transformedKdata (σlog10K)follows the logistic function of depth (Fig.16d),which reads

where σ0,d0(m) and β are the fitting parameters.Eqs.(1) and (2)describe the distribution ofKat each depth in deep-incised valley areas.

4.3.2.Excavation-induced permeability variation in fractured rocks

Starting from its initial spatial distribution given in Fig.16 and Eqs.(1) and (2),the permeability of fractured rocks may change drastically due to excavation of foundations,slopes,tunnels and underground caverns.The excavation-induced permeability change is determined by the structural alterations of fractured rocks(e.g.fracture deformation and microcracking)and influenced by a series of factors such as lithology,fracture patterns,geostress level,and the excavation method and procedure.Given the multiscale heterogeneity of fractured rocks,the upscaling/homogenization technique is indispensable for characterizing the excavation or damage-induced permeability change.This technique could be applied in the framework of elasto-plasticity by taking into account the closure,opening and shear dilation of fractures (Oda,1986;Chen et al.,2007,2014b,2015c;Zhou et al.,2008b),or in the framework of damage mechanics by considering the microcracking of rock blocks and surface degradation of rough fractures(Liu et al.,2016b).Typically,an upper-bound estimate for the excavationinduced permeability change in fractured rocks withnfamilies of preferential fractures could be written as (Chen et al.,2007;Liu et al.,2016b):

where K is the effective hydraulic conductivity tensor of fractured rock;Kris the effective hydraulic conductivity tensor of rock blocks,which could again be formulated as a function of damage variables(Chen et al.,2014b);Kf0,ef0,sfand nfare the initial hydraulic conductivity,the initial mean aperture,the mean spacing and the unit normal vector of thefth set of fractures,respectively;Δεzfis the increment of the normal strain of thefth set of fractures during loading or excavation;and δ is the Kronecker delta tensor.

Eq.(3)implies that the excavation-induced permeability change has the following important features:

(1) K is a cubic function of Δεzf,and any variation in εzfunder excavation will trigger the change in K,even in orders of magnitude.

(2) K depends on incremental strains,rather than on stresses,which makes it possible to account for the impact of nonlinear or irreversible deformation of fractures on permeability variation.

(3) The orientation nfof fractures possibly renders K highly anisotropic,even if initial isotropy is assumed for K.

Fig.17 shows the excavation-induced relaxation and permeability change in the surrounding rocks of underground powerhouse caverns at the Jinping-I Hydropower Station (Chen et al.,2015c).The surrounding rocks mainly consist of marble and breccia marble,typically containing four groups of criticallyoriented fractures and mostly rated as “fairly good” (class III1).The ratio of the uniaxial compressive strength (UCS) to the maximum stress is in the range of 1.5-4,ranking in high to extremely high in situ geostress levels.Therefore,the depth of excavation-induced relaxation is rather large,reaching up to 2-8 m in the strong relaxation zones and 6-17 m in the weak zones,as detected by acoustic wave data and borehole TV images (Fig.17a).Correspondingly,a drastic permeability enhancement by up to 2-3 orders of magnitude is predicted by Eq.(3)on the side walls of the caverns (Fig.17b).The depth where the permeability changes by a factor of 2 reaches up to 30-35 m on the upstream side of the machine hall,which is much deeper than the weak relaxation zone and hence has a significant impact on the groundwater flow around the cavern system.

Fig.17.Excavation-induced relaxation detected by acoustic wave data and borehole TV images (a) and permeability change (b) in the surrounding rocks of underground powerhouse caverns at Jinping-I hydropower station (Chen et al.,2015c).The numbers in (b) represent the ratio of enhanced hydraulic conductivity to its initial value before excavation.

4.3.3.Time-dependent effect of permeability in the foundation rocks of high dams

Another scenario of permeability change becomes important when high dams are built on sediment-laden rivers.Under sustained hydraulic gradient,fine-grained sediments (including suspension and colloid) are transported by groundwater flow and deposited in the fracture system,causing physical,chemical or biological clogging of flow channels in the foundation rocks.The fracture clogging leads to a gradual decrease in the permeability of foundation rocks as time evolves (see Figs.12d and 15b).This phenomenon is evidenced by the annually decreasing amount of discharge out of the drainage system at quite a number of largescale dam sites during the operation period,occurring even at the weirs in the upper drainage tunnels much higher than the riverbeds (Chen et al.,2021).The fracture clogging-induced permeability decrease could be best characterized with a time-dependent model.Fig.18 shows the permeability change in the foundation rocks upstream the grout curtains at the Xiluodu arch dam site,obtained by inverse modeling of the transient flow.Combined with the observations in a single fracture(Fig.15b),the clogging-induced permeability variation is found to follow the negative exponential function (Chen et al.,2021):

whereKandK0are the current and initial hydraulic conductivity,respectively;tis the time in year;andaandbare the fitting coefficients.

Eq.(4) can also be used to predict when the dynamically equilibrated state could be attained for the transient flow through dam foundations.This equation is of great significance for longterm performance assessment of high dams on sediment-laden rivers because the permeability decay reduces the pressure gradient on the impervious barriers.Furthermore,as the clogging proceeds,the permeability anisotropy becomes slightly lower(Fig.18b).

Fig.18.Variation of horizontal hydraulic conductivity (a) and the anisotropy ratio of vertical to horizontal hydraulic conductivity (b) in the foundation rocks of Xiluodu arch dam upstream the grout curtains (Chen et al.,2021).

4.4.Inertial permeability of fractured rocks

Pioneered by Forchheimer(1901),the non-Darcian flow through porous/fractured media has been widely investigated over the past century,but this theory is rarely applied to geotechnical engineering problems.The underlying reason is that the Forchheimer’s law involves two hydraulic properties (i.e.the viscous hydraulic conductivityKvand inertial hydraulic conductivityKi,Table 1),and how to properly determineKiremains an open issue.If Darcy’s law is applied in the non-Darcian flow regime,the calculated apparent permeability (Ka) is often confused with the intrinsic permeability(Kv=K),hence significantly underestimating the permeability of the media even by orders of magnitude and increasing the risk of leakage in the design of seepage-proof systems(Chen et al.,2015a;Chen,2022).This section presents three techniques to determineKiof fractured rocks,aiming at promoting the application of the non-Darcian flow theories in geotechnical engineering practices.

4.4.1.High pressure pack test and its data interpretation

High-pressure packer test (HPPT) is an enhanced packer test technique in which the maximum injection pressure(Pmax)is much higher than the magnitude of 1 MPa adopted in the conventional test.This technique has been widely used for examining the flow behaviors,hydromechanical responses and hydrofracturing risks in rock formations under high water pressure and deeply buried environment (Rutqvist et al.,1998;Guglielmi et al.,2014;Chen et al.,2015a).The HPPT is performed at each isolated interval of borehole after inflation of the packers by a stepwise rise of the injection pressure untilPmaxis attained,which is followed by a stepwise release of the injection pressure.The test at each step of constant pressure (P) proceeds until the injected flow rate (Q) becomes stabilized.The most important raw data for each test interval are theP-Qcurve,but how to properly interpret the hydraulic properties of rock formations from this curve remains unsolved for decades.

The characteristics of the HPPTP-Qcurves are closely related to the intactness of the tested rocks,the flow conditions in the fractures and the hydraulic fracturing phenomenon occurring in the test rocks(Chen et al.,2015a),and can be classified into three types(A-C;Zhou et al.,2018),as shown in Fig.19.The most prominent feature of the curves is its nonlinearity due to the occurrence of non-Darcian flow or hydraulic fracturing.The former is characteristic of a less increase in volumetric flow rate than the proportional increase in pressure (Fig.19b and e),while the latter is characteristic of an inflection point defined as the critical pressurePcat which hydraulic fracturing or fracture dilation occurs.ByPc,theP-Qcurves can generally be divided into a laminar/non-Darcian flow phase and a hydraulic fracturing phase(Chen et al.,2015a,b;Zhou et al.,2018).

Type B is the most common type of the HPPTP-Qcurves,which occurs in a wide range of rock conditions from poor to perfect intactness (Chen et al.,2015a).With the assumption of 3D spherical,steady-state flow from the point sources on the axis of a test interval,a Forchheimer’s law-based analytical model was established to interpret the hydraulic properties(Kv,Ki)of the test rocks(Chen et al.,2015b):

whererwis the radius of the test borehole;Lis the length of the test interval;andc1andc2are two coefficients determined by quadratic fit to the measuredH-Qcurves in the non-Darcian flow phase:H=c1Q+c2Q2,in whichHis the injection pressure head.

The practical significance of Eq.(5) is summarized as follows:

(1) Eq.(5) is simple and concise,and can be reduced to Hvorslev’s(1951)equation as the inertial effect vanishes(i.e.1/Ki→0 for the type A curve).It is therefore included in the industrial standard ‘Specification for Water Pressure Test in Borehole of Hydropower Projects’ (NB/T 35113-2018).

(2) Eq.(5)implies that the maximum Lugeon value in the range of injection pressure between the first-step pressure and the operating pressure should be taken for design of impervious barriers (Zhou et al.,2018).This important criterion is also included in the above industrial standard.

(3) The hydraulic properties(Kv,Ki)become pressure-dependent when hydraulic fracturing occurs(i.e.P>Pc)or for the type C curves.Formulae for this behavior can be found in Chen et al.(2015b),Chen (2022) and Zhou et al.(2019a).

Fig.19.Types and the typical cases of HPPT P-Q curves: (a,d) Type A,(b,e) Type B,and (c,f) Type C (Chen et al.,2015a;Zhou et al.,2018).

(4) When HPPTs are performed in multiple test intervals or boreholes,the values ofKvandKiof a rock formation could be more representatively determined by using the meanP-Qcurve over all of the tests,though they may also be directly estimated by the means ofKvandKi(Chen et al.,2020a;Chen,2022).

(5) Eq.(5) applies to HPPTs in vertical boreholes,and a generalized formula was also proposed for inclined and horizontal boreholes (Li et al.,2021).

4.4.2.Constant-rate pumping test and its data interpretation

Constant rate test (CRT) is an important tool to assess the hydraulic properties of aquifers,where water is pumped from a well at a constant flow rate(Q)and the drawdown(s)is measured in the observation wells.In the vicinity of the pumping well,the flow velocity increases and the non-Darcian effect may become significant,especially at higher pumping rate (Mathias et al.,2008).Furthermore,the accuracy of CRT data interpretation may be reduced by the assumed one-dimensional (1D) linear or 2D radial flow geometry for the heterogeneous aquifers (Theis,1935;Miller,1962).On the basis of Barker’s (1988) generalized laminar radial flow model that integrates the non-integer dimension(1 ≤n≤3)of flow geometry,a generalized non-Darcian radial flow (GNRF)model was developed for data interpretation of CRTs (Liu et al.,2016a,2017b).The proposed model is based on Forcheimer’s law,and hence takes into account both the non-Darcian effect of flow(byKi)and the fractional flow dimension of aquifer(byn).The noninteger flow dimensionnis of hydrogeological significance,because it represents the fractal dimension of flow channels or multi-scale heterogeneity of aquifers (Le Borgne et al.,2004).The GNRF model can be simply written as (Liu et al.,2017b):

wheresis the drawdown;nis the flow dimension;ris the distance from the pumping well;tis the time;γ is the Euler’s constant;anda1-a3andb1-b3are the coefficients related to the hydraulic properties (Kv,Ki,Ssandn),the pumping flow rateQand the thickness of the aquiferT,respectively.

Fig.20.Type curves for the GNRF model (Liu et al.,2017b),where Γ denotes the complementary incomplete gamma function.

Fig.21.A universal power-law relationship between the inertial permeability and the viscous permeability of geological media: (a) The collected dataset,and(b) A subset of data obtained under deformation(Zhou et al.,2019b).κv and κi are the viscous and inertial permeabilities,respectively,with Kv=(ρg/μ)κv and Ki=gκi,in which ρ and μ are the density and dynamic viscosity of fluid,respectively;and g is the gravitational acceleration.

Fig.20 shows the type curves of the GNRF model,and the hydraulic properties (Kv,Ki,Ssandn) could be identified by curve fitting to the late-time drawdown data(Liu et al.,2017b).Numerical verification shows that the GNRF model behaves well forn>1.5,but it becomes less reliable forn≤1.5.As a remedy for this defect,a Forchheimer’s law-based flow model with 1D linear flow pattern was further proposed (Chen et al.,2019).With these proposed models,not only the non-Darcian flow properties (Kv,Ki,Ss),but also the possible transition of flow dimension (n) in the aquifers during a hydraulic test could be demonstrated (Chen et al.,2019;see Section 8.2).An important remark is that when the non-Darcian effect becomes significant,the intrinsic permeability (Kv) will be underestimated and the specific storage(Ss)will be overestimated by Darcy’s law-based models.

4.4.3.A universal relationship for estimating inertial permeability

The non-Darcian flow properties (Kv,Ki) are in nature determined by the geometries of the geological media.Due to the highly complex and diverse pore and fracture geometries as well as their variations (e.g.due to hydromechanical coupling),it is difficult to establish a purely theoretical relationship betweenKvandKithat is applicable across all types of fractured media.When the in situ hydraulic test data are absent,is there any correlation between the viscous hydraulic conductivityKvand the inertial hydraulic conductivityKisuch thatKican be statistically predicted from the more easily and cheaply obtained value ofKv? By compilation and statistical analysis of over 4000 pairs of(Kv,Ki)data on single fractures,porous media,fracture networks and fractured porous media(Fig.21a),a universal power-law upscaling betweenKiandKvwas established (Zhou et al.,2019b):

where ω is a scaling coefficient.

This empirical relationship is considered‘universal’,because it is statistically valid for a large variety of more than 100 types of porous and fractured media,from single pores to site-scale fractured rocks,withKvvarying over 12 orders of magnitude.It has also been found that Eq.(7)remains applicable for deformable media as long as the void structures maintain some kinds of similarity during deformation(Zhou et al.,2019b),as shown in Fig.21b.The practical significance of Eq.(7)is that once a medium has a knownKvand the scaling coefficient ω is determined by a limited number of tests or site conditions,one can estimateKiwith reasonable statistical uncertainty(Chen et al.,2020a).Besides,Eq.(7)can be further used to characterize the permeability variation of pressurized fractured rocks where the inertial effects and fracture dilation are intertwined (Zhou et al.,2019a).

4.5.Unsaturated hydraulic properties of fractured rocks

Unsaturated flow is an important process in the vadose zone,which occurs not only in soil covers,but also in fractured rocks.But unlike soils,fractured rocks are characteristic of much stronger heterogeneity and larger representative elementary volume (REV)size.Consequently,the suction measurement techniques originally developed for soils fail to apply for fractured rocks at the field scale,and how to properly determine the water retention and relative permeability relations for fractured rocks remains an extremely challenging issue.As validated by numerical upscaling and inverse modeling,it is now widely accepted that the VG (van Genuchten,1980) model also applies to fractured rocks (Liu et al.,1998;Peters and Klavetter,1988;Chen et al.,2020b).

By statistically analyzing a set of data on more than 1500 soils,rock cores,fractures and fractured rocks,upscaling relationships between the VG parameters (αVG,nVG) and the hydraulic conductivity (K) were proposed to predict the upper bounds of the VG parameters from the more easily-obtained value ofK(Chen et al.,2022),as shown in Fig.22.The proposed upper bounds of αVGandnVGare sufficiently narrow,hence their representative values can be determined more easily by inverse modeling with the groundwater observations in saturated zones.Combined with the continuum approach,the complex saturated-unsaturated flow behaviors,such as multiple water tables,wedge-shaped unsaturated zones and preferential flow in backbone structures,could be satisfactorily modeled (Fig.23).

5.Groundwater flow modeling through fractured rocks

Numerical modeling of groundwater flow through porous/fractured media is of great importance in geotechnical engineering for optimization design and performance/safety assessment.The nonlinear nature of the groundwater flow models,however,makes it difficult to be efficiently solved,especially in large-scale problems under difficult geological conditions.This section presents the variational inequality approach that effectively overcomes the nonlinearity originated from boundary conditions and can be applied for modeling of steady-state,transient and saturatedunsaturated flows in a unified framework.

5.1.The unified variational inequality method for groundwater flow simulations

There are two sources of nonlinearity encountered in groundwater flow modeling.One is originated from constitutive relations,such as the deformation-or time-dependent permeability variation (Section 4.3),Forchheimer’s law for non-Darcian flow (Section 4.4),and water retention and relative permeability curves for unsaturated flow (Section 4.5).The other is originated from the third-type boundary conditions on the rainfall infiltration,evaporation and seepage boundaries.These two sources of nonlinearity are intertwined,resulting in severe numerical instability of the solutions.It has been proved that on the thirdtype boundaries,the pressure head and flux uniformly satisfy the unilateral condition of Signorini’s type (Zheng et al.,2005;Borsi et al.,2006;Chen et al.,2011b;Hu et al.,2017a).Fig.24 shows the unilateral boundary conditions for the saturatedunsaturated flow,which are not only highly nonlinear,but also have complicated transitions of the types of boundaries during the rainfall-evaporation cycles (Chen,2022).

The groundwater flow models with unilateral boundary conditions can be rigorously formulated and solved with the variational inequality (VI) method.The VI method decomposes the unilateral conditions into a flux part and a pressure head part,with the former being treated as a natural boundary condition and the latter as a coercive boundary condition.Hence,the difficulty in selecting the trial functions is significantly reduced during numerical discretization (Zheng et al.,2005).The unified VI formulation of the steadystate,transient and saturated-unsaturated flow models (Table 1) is given by(Zheng et al.,2005;Chen et al.,2011b;Hu et al.,2017a):

Fig.22.Upscaling relationships for the upper bounds of VG parameters (Chen et al.,2022): (a) αVG,and (b) nVG.

Fig.23.Saturated-unsaturated flow in the right bank slope at Baihetan arch dam site (Chen et al.,2022): (a) Geological cross-section,and (b) The multi-layer distribution of groundwater in natural state simulated with a continuum model.

Fig.24.The unilateral boundary conditions and the transitions of boundary types for saturated-unsaturated flow (Chen,2022).

wherehisthe solution;? isthetrial function thatsatisfies the coercive boundary conditions(including the Dirichlet boundary condition);Ω isthe domain;Γq,ΓfandΓTare the flux,freesurfaceand the third-type boundaries,respectively;μ*is the specific yield;is the maximum allowable value of fluxon ΓT,with=0 on the seepage face,=Ion the infiltration boundary,and=Reon the evaporation boundary.Other symbols can be found in Table 1 and Fig.24.

For unconfined steady-state flow,Eq.(9a) and the second and third terms on the right hand side of Eq.(9c) vanish;and for unconfined transient flow,Cw=0 and the third term on the right hand side of Eq.(9c) vanishes.In these two flow models,the VI formulation is established by introducing a Heaviside functionHof pressure head and extending the flow laws (Darcian or Forchheimer) from wet to dry domains(Zheng et al.,2005;Chen et al.,2011b).The free surface between the wet and dry domains is hence transitioned to an inner boundary,and once the unilateral boundary conditions are iteratively determined,the location of the free surface is also determined by the condition on it.For saturated-unsaturated flow,the term 1-Hin Eq.(9b)and the second term on the right hand side of Eq.(9c)vanish.With the VI formulation,the strong nonlinearity and the complex transitions of the unilateral boundaries can then be rigorously and efficiently tackled (Hu et al.,2017a).Details of the algorithms can be referred to Chen(2022).

5.2.Inverse modeling of groundwater flow through fractured rocks

Inverse modeling plays an indispensable role in determining the hydraulic properties of rock formations on a scale representative of the site conditions.This technique can also be used to determine the boundary conditions when the flux on the boundaries of the domain is unknown (Chen et al.,2016c;Hong et al.,2017).It is performed by minimizing an objective function that measures the fitness between in situ measurements and numerical simulations.However,the inverse modeling is inevitably plagued by the problem of non-uniqueness,especially when only the groundwater level (or piezometric head) observations are used or the steadystate flow conditions are assumed (Zhou et al.,2015a).To improve the accuracy of the inverse solution and to account for the time-dependent nature of hydraulic properties and boundary conditions,a weighted objective function or double objective functions could be defined using the time-series measurements of both piezometric head and discharge (Zhou et al.,2015a;Chen et al.,2016b,c,2020b;Hong et al.,2017).

The use of multiple types of time-series observations increases the difficulties in multi-objective optimization,especially in largescale inverse problems under complex site conditions.The computational cost of inverse modeling mostly depends on the complexity of groundwater flow models and how many times the numerical simulations are invoked.A practical procedure requires that the number of forward simulations is sufficiently small.For this purpose,an inverse modeling procedure combining orthogonal design(OD),finite element analysis(FEA),artificial neural network(ANN)and genetic algorithm(GA)was proposed(Zhou et al.,2015a;Chen et al.,2016b).In this procedure (Fig.25),OD is employed to construct a small set of representative values for the hydraulic parameters from a pre-defined parameter space.The hydraulic parameters to be back-calculated may include the hydraulic properties (e.g.K,Ss,αVGandnVG) and/or the parameters preset for boundary conditions.FEA is then invoked for groundwater flow modeling on each parameter combination,which yields the simulated time-series data at the observation points.An implicit mapping from the parameter space to the seepage fields is obtained by constructing and training a back propagation neural network(BPNN) model.Finally,GA is used to obtain a globally optimal estimate of hydraulic parameters such that a best fit between the field measurements and the numerical results calculated by the trained BPNN model is ensured.

This procedure works well with balanced computational cost and accuracy,and hence has been widely applied to addressing the field-scale problems(Zhou et al.,2015a;Chen et al.,2016b,c,2020b,2022b;Hong et al.,2017).Furthermore,when this procedure is performed annually or periodically,the variations in hydraulic properties and boundary conditions could be back-calculated as time evolves(Chen et al.,2021).

5.3.Regional-scale modeling of groundwater flow in aquiferaquitard systems

To save computational cost,a site-scale model is commonly created for numerical simulation of groundwater flow in engineering practices,in which the geological settings and geotechnical structures could be well represented.The dimension of the sitescale model,however,is generally much smaller than the dimension of the hydrogeological unit (HU) where the site is located,as shown in Fig.26.A difficult problem arises about how to determine the conditions on the lateral boundaries of the computational domain,across which groundwater exchange occurs but the flux is unknown.This problem is further complicated when complex aquifer-aquitard systems are present in the region,where the flux or water head on the lateral boundaries likely becomes discontinuous.

Four strategies were proposed to roughly determine the lateral boundary conditions (Chen,2022).The first strategy is to locate the lateral boundaries on local watershed or local impervious structures,so that a no-flux boundary condition can be assumed.The second one is to prescribe a water head boundary condition,provided that borehole groundwater level observations are available near the boundaries.The third one is to assume a distribution of groundwater level on the lateral boundaries,and determine the coefficients of the function by inverse modeling (Chen et al.,2016c;see Section 5.2).The last and the best strategy is to construct a site-scale and a regional-scale model,so that these two models of hierarchy could be coupled and the flux or water head conditions on the lateral boundaries of the site-scale model could be rigorously determined (see Section 8.1).By doing so,the groundwater flow behaviors both on the regional scale and on the local site scale can be better understood at an acceptable computational cost.

Fig.25.The inverse modeling procedure for groundwater flow through rock formations (Zhou et al.,2015a;Chen et al.,2016c).

Fig.26.The unknown conditions on the lateral boundaries of a numerical model taken from the HU (Chen,2022).

5.4.Modeling of groundwater flow in karst aquifer systems

Karst aquifers are widely distributed in the world.In karst areas,many geotechnical engineering applications,especially tunnel constructions,have a significant impact on the hydrogeological environment.Karst aquifers are characterized by strong hydrogeological heterogeneity and complex flow mechanisms (coexistence of slow flow in the fractured rock matrix and fast,and often turbulent flow in the large conduits),which makes prediction of groundwater flow in karst area a challenging task(Hartmann et al.,2014).Numerical approaches for modeling groundwater flow in karst aquifers by solving the partial differential equation (PDE) for Darcian and/or non-Darcian flow include the equivalent porous medium(EPM)model,the double porosity model,and the coupled discrete-continuum model (Shoemaker et al.,2008;see Fig.1).Among the different approaches,the coupled discrete-continuum model explicitly considers the large conduits or underground rivers,and thus can better represent the hydrogeological features in karst aquifers.From a hydrological and water balance modeling perspective,the numerical simulation approach can be classified as distributed models.Another type of approach is the dimensionless lumped parameter model(Fleury et al.,2009;Schmidt et al.,2014;Zhou et al.,2021).In this approach,the spatial variability in the aquifer layers is not considered and no spatial discretization is needed.Rather,linear or nonlinear relationships are used in the lumped parameter model to describe the storage and discharge in separate model compartments according to the physical processes affecting water balance.Thus,the lumped parameter model is suitable for groundwater simulation with limited information on field characterization of hydrogeological conditions (Zhou et al.,2021).This section presents the lumped parameter approach and the coupled discrete-continuum approach.

5.4.1.Lumped parameter models for karst groundwater balance modeling

As shown conceptually in Fig.27,in a karst HU,the water from precipitation infiltrates to the karst system through diffuse flow in the soil cover and quick flow in sinkholes and vertical shafts.The diffuse flow recharges the epikarst zone (with a low degree of karstification),while the quick flow provides direct recharge to the underground river (Bailly-Comte et al.,2010).Then,the groundwater flows vertically into the more karstified vadose and phreatic zones.This can be in the form of base flow through the fissured matrix,quick flow through conduits and vertical shafts,and/or intermediate flow through enlarged karstified fractures(Zhou et al.,2021).The groundwater finally discharges either at the spring (or outlet of the underground river) or through the hyporheic zone of the surface rivers.Construction of engineering structures such as tunnels adds a new form of discharge.

Based on the conceptualized model,a modified lumped model(MLM) can be developed for karst systems,with a special concern on the discharge from the karst system to the engineering structures.As shown in Fig.27,the karst system is divided into four compartments,i.e.the soil moisture compartment (SMC),the epikarst compartment (EPC),the groundwater compartment (GWC),and the channel drainage compartment (CDC).The two compartments,SMC and EPC,account for the hydrological functioning,storage and production of the soil cover and the epikarst zone,respectively(Hosseini et al.,2017).The GWC includes three parallel sub-compartments for modeling of base flow,intermediate flow,and quick-flow,respectively.The CDC accounts for both the flow in the underground river and the influence of tunnel engineering on karst aquifer discharge.Water balance equations can be written for each compartment/sub-compartment,and these equations are linked by the flow or discharge terms with model parameters characterizing the temporal dynamics.The model parameters can be calibrated based on time-series data of e.g.precipitation and spring discharge.Detailed model equations and parameter calibration procedures can be found in Zhou et al.(2021).The MLM considers the exchange between tunnels and the aquifer system,and thus provides a useful tool for fast evaluation of karst groundwater budget changes induced by engineering activities such as tunneling and mining.

5.4.2.Coupled discrete-continuum modeling for groundwaterflow in karst aquifers

The steady-state and transient groundwater flow in continuum matrix for fractured rocks can be described by the equations presented in Section 5.1.However,these equations are inadequate when the flow in large conduits or underground rivers needs to be simulated.To address this issue,the coupled discrete-continuum model can be adopted.A commonly used coupled discretecontinuum model for groundwater flow in karst aquifer is the conduit flow process (CFP) model (Shoemaker et al.,2008).In this approach,both karst underground rivers and tunnels can be treated as discrete conduit elements with idealized cross-sectional shapes of circular pipes (Zheng et al.,2021).The conduits or conduit networks are composed of discrete segments and nodes embedded in the matrix grid (see Fig.28).The matrix is used to represent the fractured rock as an equivalent continuum.The water exchange between conduits and matrix occurs at the nodes,and for each given node,a water mass balance equation can be used to couple the flows in the conduit and in the matrix.Depending on the Reynolds number,laminar or turbulent flows can be modeled in the conduits.Both cases of fully saturated and partially filled conduits can be considered.The water exchange between the fractured rock matrix and the conduit can be computed based on the conduit conductance (which is a function of wall rock permeability and filling status of the conduit) at the embedded element and the hydraulic head difference between the conduit node and the matrix element.More details can be found in Zheng et al.(2021).

Fig.27.(a) Conceptualization of hydrological processes in a karst HU,and (b) Structure of the modified lumped parameter model for groundwater balance in a karst aquifer with tunnels.Modified from Zhou et al.(2021).

Fig.28.Schematic diagram of the coupled discrete-continuum model for modeling groundwater flow in karst aquifers.

6.Optimization design of seepage control systems

Seepage control is critical to maintaining the stability and safety of geotechnical structures.The optimization design of seepage control measures is enabled on the premise that the performance could be accurately quantified.This section summarizes the control mechanisms,numerical simulation and optimization design of seepage control systems.

6.1.Mechanisms of seepage control in geotechnical engineering

The goals of seepage control are to limit the amount of leakage through,reduce the pore pressure in and increase the stability of geological media and geotechnical structures.They are achieved in geotechnical engineering using independently or cooperatively the impervious barriers,drains and filters.Other techniques,such as the control on the rate of surface water level rise or drop,are also used to regulate the groundwater level and the stability of slopes and foundations.Based on the groundwater flow model consisting of the governing equation and the initial and boundary conditions,it has been demonstrated that the seepage control mechanisms can be classified into the following four categories (Chen et al.,2010;Chen,2022):

(1) Control by hydraulic properties.The impervious barriers and filters work by changing the permeability and seepage resistance of materials in the flow domain;

(2) Control by boundaryconditions.The drains and the control on the rising/dropping rate of surface water level function by creating new flux or seepage boundaries inside the domain or changing the water head boundary condition of the domain;

(3) Control by coupled processes.The stress/deformation,thermal and chemical processes have an impact on the storage terms or hydraulic properties of the continuity equation,and hence the control on these coupled processes would to some extent regulate the flow process;

(4) Control by initial conditions.The initial conditions have an important impact on the early evolution of the flow process,and they can be regulated in some particular situations,e.g.by changing the initial water content in the clay barriers.

The above mechanisms show that each engineering measure can be rigorously represented in the groundwater flow model,hence laying a solid basis for performance assessment and optimization design of the seepage control systems.It should be noted that the first two mechanisms are much more sensitive and effective for seepage control,as evidenced by their wide uses in geotechnical engineering (Chen et al.,2010).

6.2.Refined modeling and performance assessment of seepage control systems

A seepage control system typically consists of impervious barriers(e.g.grout curtains,clay cores,concrete face slabs),drains (e.g.drainage hole arrays,tunnels and galleries)and filter zones.Seepage control systems can be extremely complex in a large-scale hydropower project.A typical example of a complex seepage control system that contains 7980 drainage holes and grout curtains covering an area of 8.7×105m2in the foundation of Xiluodu arch dam(Chen et al.,2021)is shown in Section 8.1.The impervious barriers and filter zones can be relatively easily modeled with solid elements as long as their hydraulic properties can be properly determined.The drainage hole arrays,however,are characteristic of great number,dense spacing,small diameter and large length,which makes the generation of numerical grids quite difficult.The modeling is further complicated by their complex boundary conditions,which are summarized as follows(Chen et al.,2008;Chen,2022):

(1) Water head boundary condition (Fig.29a).The drainage holes deployed in the lowest level of drainage tunnels satisfy this boundary condition when overflow occurs.

(2) No-flux boundary condition(Fig.29b).If no overflow occurs,the above drainage holes should be treated as no-flux boundary.

(3) Unilateral boundary condition of Signorini’s type (Fig.29c).The vertical or inclined drainage holes installed between two horizontal drainage tunnels or drilled from the ceilings of drainage tunnels satisfy this boundary condition.

(4) Mixed boundary condition (Fig.29d).This boundary condition is associated with a draining well with known pumping rate.

Fortunately,the above two difficulties in the modeling of drainage hole arrays,i.e.the grid generation and boundary condition treatment,can be rigorously and effectively addressed by the unified VI method,together with a substructure technique(Chen et al.,2008,2011b;Hu et al.,2017a).The substructure technique was originally proposed by Wang et al.(1992)and then improved by Zhu and Zhang(1997).This technique automatically refines the grids to model the geometry of drainage holes and condenses the internal degrees of freedom,hence enabling the performance assessment of complex drainage systems containing tens of thousands of drainage holes.Fig.30 shows examples of the performance of drainage holes deployed in the Guangzhao roller compacted concrete dam and around the Jinping-I underground powerhouse caverns.

6.3.The optimization design procedure for seepage control systems

The optimization design of seepage control systems becomes feasible as long as their performance can be rigorously quantified.As shown in Fig.31,the goal is to select a set of design parameters for grout curtains (e.g.extension,depth,the rows and spacing of grouting holes,and permeability)and drainage holes(e.g.extension,depth,the number of rows,and spacing),such that a good balance between the cost of construction and the desired performance(e.g.a limited amount of leakage,lowered uplift pressure or pore pressure,and enhanced seepage resistance) is achieved.Fig.32 shows the procedure for optimization design of a seepage control system,which is realized on the premise that the site conditions are sufficiently characterized (Li et al.,2017).Particularly,the distribution,thickness and permeability of aquitards and the development patterns of concentrated leakage paths(e.g.faults,fracture zones,shear zones and bedding planes)play a critical role in figuring out an initial layout of the seepage control system.The layout can then be refined according to the design criteria and engineering experiences,and iteratively updated with an optimization rule until the desired performance is obtained with a balanced construction cost.With the aid of the refined modeling technique (Section 6.2),the optimization procedure can improve the performance of seepage control at much reduced cost,even under the site conditions with uncertainties and complexities.

7.Groundwater flow-induced deformation and its control

Groundwater flow can change the effective stress and induce mechanical deformation of fractured rocks,and the deformation can in turn affect the flow field through modifying the hydraulic properties of the medium (Chen et al.,2009,2014a,2015c),resulting in coupled hydro-mechanical processes.Large-scale hydraulic and hydropower engineering activities amplify the coupling effects.Under the influences of reservoir impoundment and alternating water level,the rocks and soils in the reservoir area can experience rapid and significant changes in the flow field,which can accelerate the deformation.This section mainly focuses on reservoir impoundment-induced deformation phenomena and their mechanisms in dam engineering as well as how to mitigate the deformation,while the discussion is also highly relevant for other similar applications in geotechnical engineering.

Fig.29.The boundary conditions of drains: (a) Water head boundary condition,(b) Impervious boundary condition,(c) Seepage boundary condition,and (d) Mixed boundary condition (Chen et al.,2008;Chen,2022).

Fig.30.The effects of drains on groundwater flow:(a-c)Through a gravity dam(Chen et al.,2008)and(d)around an underground powerhouse cavern system(Chen et al.,2015c).

Fig.31.Layout of a grout curtain and drainage holes in dam foundation (Li et al.,2017).

Fig.32.Flowchart of optimization design for seepage control systems (Li et al.,2017).

Fig.33.(a) Valley contraction monitoring lines,(b) Schematic diagram of the rock layers and vertical locations of the monitoring lines,and (c) Evolution of valley contraction at different monitoring lines after reservoir impoundment.

7.1.Reservoir impoundment-induced deformation phenomena in dam engineering

Typical reservoir impoundment-induced deformation phenomena include reservoir base subsidence or uplift,valley contractive deformation,wetting and seepage induced deformation in rockfill dams and foundations,and bank slope deformation and landslide.Among these phenomena,the latter two are mostly resulting from deformations at the local scale within the dam area or a specific unstable slope.Reservoir-induced landslides have been extensively studied within the fields of engineering geology and geohazards (Tang et al.,2019),and the seepage-induced deformation of rockfill dams and foundations is relatively well understood and can be controlled in engineering practice (Chen et al.,2011a;Wang et al.,2016).Reservoir impoundment after construction of high dams results in sharp rise in groundwater table in the aquifers,which can exert a wide-spread impact on the regional hydrogeological conditions.Unlike the latter two phenomena,the former two present relatively new issues in high dam engineering and has gained significant attention due to their direct threat to the safety of high dams(Cheng et al.,2017;Sun et al.,2018;Jiang et al.,2020;Yang et al.,2022).

Notable examples of deformation in forms of valley contraction and dam base uplift/subsidence include Xiluodu,Jinping-I,Jiangya,and Tongjiezi hydropower projects (Wu and Qi,2003;Liu et al.,2011;Cheng et al.,2017;Jiang et al.,2020).The dam base at Jiangya experienced uplift deformation both at the left and right banks after reservoir impoundment,while the reservoir base deformation pattern at the Tongjiezi dam exhibited subsidence at the left and uplift at the right bank.The reservoir impoundment of the Jinping-I arch dam (305 m) and the Xiluodu arch dam (285.5 m) induced uplift/subsidence deformation and more importantly valley contraction (Cheng et al.,2017),as revealed by data from continuous monitoring of the deformation at the reservoir banks.Field data show that the contractive deformation at the Xiluodu arch dam has reached up to 100 mm within 7 years after initial impoundment (see Fig.33).Whether the contractive trend at the Xiluodu dam site will converge in the near future and what are the exact mechanisms behind this large deformation are still under debate(Cheng et al.,2017;Jiang et al.,2020).

Particularly,the deformation at the Xiluodu dam site shows some typical characteristics.First,the deformation exhibits strong correlation with the reservoir impoundment cycles and water levels.Second,the deformation occurs both in the upstream and the downstream of the dam.Third,the deformation extends in a wide area and no clear evidence of tensile failure in the rocks can be observed.Mechanistically,the valley deformation and dam base uplifting/subsidence are associated with hydromechanical coupling effects induced by hydrogeological condition variations at the regional scale due to reservoir impoundment of high water levels.In particular,factors such as the occurrence of rock layers,lithology,surface terrain,and groundwater recharge-discharge patterns can all play a role.However,the exact underlying mechanisms,which are strongly dependent on the regional geological/hydrogeological conditions for each specific dam site,remain not fully understood.Due to the complex geological and hydrogeological conditions at the regional scale and intricacies in the coupled processes,more dedicated investigations are still needed to fully resolve the valley contraction problem.Special attention should be paid to the largescale unfavorable geological structures or layers that are mechanically weak and prone to physical/mechanical alteration and timedependent variations of rock properties due to the reservoir impoundment effects.

7.2.Mechanisms and modeling of groundwater flow-induced deformation in fractured rocks

By extending Boit’s poroelasticity theory (Biot,1941),the momentum conservation equation of continuum mechanics for coupled hydro-mechanical processes can be written as(Zhou et al.,2008a;Hu et al.,2016;Chen et al.,2022):

where u is the displacement vector;α is the Biot coefficient (Biot and Willis,1957);Depis the fourth-order tangential elastoplastic modulus tensor;δ is the Kronecker delta tensor;g is the gravitational acceleration;ρeff=(1 -φ)ρs+φρwis the effective density,with φ being the porosity;and ρsand ρware the densities of the solid and water,respectively.

According to Eq.(10),groundwater flow induces deformation of the fractured medium mainly through the following mechanisms(Chen,2022).First,the groundwater flow can cause pore pressure(p) changes and thus can change the effective stress field;also,water pressure can be applied on the boundaries of the fractured medium.Both effects can lead to deformation of the fractured rock.Second,through softening and internal erosion,groundwater can also act to deteriorate the rock mechanical properties (e.g.elastic moduli,cohesion,internal friction coefficient),which leads to decreased Depand thus increased deformation with possible consequences of yield,creep and ultimate failure.Third,groundwater flow can also induce changes in the effective density of ρeff,affecting the stress field and deformation.With the aid of coupled modeling,these mechanisms can provide guidance for the control and mitigation of flow-induced deformation.

Due to the complex mechanisms of coupled processes and the highly nonlinear characteristics of the constitutive relations in coupled models,it is necessary to develop efficient numerical models to perform prediction and assessment of groundwater flow-induced deformation in fractured rocks in various geotechnical applications.Numerical simulators for modeling coupled processes can be categorized into continuum models and discontinuum models.Large-scale geotechnical application problems often rely on continuum model simulations which have the advantage of computational efficiency since the hydraulic and mechanical properties due to the existence of fractures are represented by equivalent parameters.The finite element method(FEM)can be used to numerically solve the fully coupled governing equations with appropriate boundary and initial conditions as well as constitutive models for the tangential modulus tensor and the stress-dependent hydraulic conductivity tensor(Zhou et al.,2008a;Chen et al.,2009).Detailed information of the geological/hydrogeological conditions and field observations of pore pressure,flow rate and deformation are required to calibrate the numerous parameters in the coupled processes and to ensure reliability of the modeling results.

7.3.Engineering measures for mitigation of groundwater flowinduced deformation

How to guaranty stability of reservoir banks/foundations under the threat from groundwater flow-induced deformation is thus a primary task to ensure safety of engineering structures.For localscale bank slope deformation,mechanical measures such as reprofiling and reinforcement systems (e.g.anti-slide piles and anchors)are often adopted and accompanied by surface and internal drainage systems.For stability of underground engineering structures and foundations,seepage control is a key issue.Particularly,due to the areal extent or scale of the reservoir impoundmentinduced deformation,it is difficult to control the contraction through stabilizing or reinforcement systems at the bank slopes.Effective engineering measures to control valley contraction require a fundamental understanding of the underlying physical mechanisms which are site-specific and can be complex to pinpoint.To gain such an understanding and to predict the contractive deformation,one needs a refined coupled hydromechanical model with representative parameters,which faithfully considers the geological and hydrogeological conditions,the dam engineering activities,and the reservoir operating procedures.The evolution of regional hydrogeological environment after reservoir impoundment should be carefully studied (Chen et al.,2020b).More attention should be given to the weak zones/layers that are sensitive to water-induced softening and pore pressure changes and are hence able to undergo large deformation.On the basis of refined numerical modeling,strategies of constraining valley deformation may be developed by coordinated efforts of site selection,dam engineering design,reservoir operation,and seepage control systems.More field measurements(e.g.additional boreholes),laboratory experiments,and more comprehensive monitoring programs (e.g.combining ground measurements and satellite images) may be useful to further constrain the physical models and to establish early-warning systems for effective mitigation.

8.Case studies

The proposed methods,models and techniques (Sections 3-7)have been applied to optimization design of seepage control systems (e.g.Li et al.,2014,2017),inverse modeling of leakage events(e.g.Zhou et al.,2015a;Chen et al.,2016b),long-term performance assessment of impervious barriers(e.g.Chen et al.,2021)in tens of large-scale dam,underground powerhouse cavern and deep-buried hydraulic tunnel projects.This section presents three typical examples of applications to Xiluodu high arch dam(Chen et al.,2021),Changheba clay-core rockfill dam (Zhou et al.,2015a;Chen et al.,2019) and Jiayan water diversion tunnels (Zhou et al.,2021;Zheng et al.,2021).

8.1.Reservoir operation-induced hydrogeological changes at the site of Xiluodu arch dam

As the fourth highest dam of the same type in the world,the Xiluodu arch dam is situated in a deep-incised valley where the Jinsha River turns 90°from northeast to southeast (Fig.34).The arch dam is roughly located at the center of the Yongsheng synclinorium basin (Fig.34a),which constitutes an independent HU with groundwater recharged by precipitation and snow melt and discharging to the Jinsha River.The catchment area and the total capacity of the Xiluodu Reservoir are 4.54 × 105km2and 1.27×1010m3,respectively.For such a giant reservoir,the reservoir filling and operation significantly impact the groundwater flow systems and consequently induce deformation.After the initial reservoir impoundment,a remarkable valley contraction was observed unexpectedly at both the upstream and the downstream of the dam site.The deformation grows continuously and reaches about 108 mm,posing potential threats to the dam safety.As discussed in Section 7,the deformation of valley contraction is mainly induced by the hydrogeological changes at the basin scale (Jiang et al.,2020;Li et al.,2020a).Therefore,it is critical to characterize the hydrogeological changes induced by reservoir filling and operation,which is fundamental to predicting and controlling valley contraction deformation.

Fig.34.Hierarchical modeling of groundwater flow in the Yongsheng basin where the Xiluodu arch dam is located:(a)The independent HU,(b)The basin-scale model,(c)The sitescale model,and(d)The seepage-proof and drainage systems represented in the site-scale model(Chen et al.,2021).The basin-scale model(b)and site-scale mode(c)are coupled by the flux conditions on the lateral boundaries,and the site-scale model (c) and the drainage system are linked by the boundary conditions on the drain walls (see Section 6.2).

Fig.35.Hierarchical modeling results of the hydrogeological changes induced by reservoir filling and operation at the Xiluodu dam site.(a)Comparison between the simulated and measured piezometric head at the observation boreholes in the confined aquifer.Reservoir operation-induced changes in piezometric head at the top of the confined P1y aquifer on:(b) 1 January 2013,(c) 1 January 2016,(d) 1 January 2019,and (e) 1 January 2027.The dam site is marked with a purple square.

The hierarchical numerical modeling approach,as depicted in Section 5.3,is essential to evaluate groundwater flow behaviors.This is because the hydrogeological changes are not only influenced by the lateral recharge from the basin,but also controlled by the local geological settings,the rise of reservoir level,and the seepage control system consisting of grout curtains and numerous drains(Fig.34d).To fully model the groundwater flow behaviors,a basinscale model (～34 km,Fig.34a and b) and a site-scale model(～4 km,Fig.34c and d) were created.The basin-scale simulation was performed with MODFLOW to capture the overall flow behaviors in the aquifer-aquitard systems from 1 January 2012 to 31 December 2030.In this simulation,the base and lateral boundaries were prescribed with no-flow conditions,and the ground surface was set as the infiltration condition.Meanwhile,the site-scale simulation was performed with the FEM to model the groundwater flow near the dam foundation.These two models of distinct scales were linked to each other by the lateral boundary conditions prescribed with the water head time-series in the basin-scale model output (see Section 5.3).The hydraulic properties of aquifers and aquitards in the study area were determined by hydraulic tests and inverse calibrations (Chen et al.,2018;2020b,2021).

Fig.35a presents the simulated and observed groundwater levels in the confined limestone aquifer.The results confirmed that the hierarchical simulation reasonably captured the groundwater dynamics around the dam site in recent three years.The mean absolute errors between the simulation and measurements for the boreholes HZK04,HZK02 and XZK03 were 10.7 m,14.86 m and 14.64 m,respectively.The small discrepancy is possibly due to the local hydrogeological settings.To further quantify the hydrogeological changes,Fig.35b-e shows the spatial distribution and the evolution of the piezometric head at the top of the confined P1y aquifer.The results demonstrated a continuous increase of the piezometric head and a drastic expansion of the influence zone until the end of 2019 (Fig.35d).The dimension of the influence zone was predicted to be about 7-9 km from the reservoir.After then,the rate of the hydrogeological changes seems to decrease,and the flow system is predicted to reach a dynamically stable state at the end of 2027 (Fig.35e).In summary,the application of the hierarchical modeling approach provided a thorough evaluation of the hydrogeological changes across the basin scale.The results further imply that the convergence of valley deformation at the Xiluodu site may be reached by around 2027,and therefore continuous monitoring and further investigations are still needed.

8.2.The leakage through Changheba rockfill dam during construction

The Changheba dam,located in Kangding,Sichuan Province,China,is a gravelly clay-cored rockfill dam with a height of 240 m(Fig.36).The bedrocks at the site consist of igneous rocks(granite and quartz diorite),ranging from moderately weathered to fresh below theground.Thefoundationrocksare blocky,withfoursetsofcritically oriented fractures.The rocks are covered by four layers of unconsolidated deposits of alluvial or fluvioglacial origins,with a maximum thickness of 79.3 m(Fig.36b).During excavation of the loose deposits(with a maximum depth of 32 m) in the dam foundation,an unexpected leakage occurred,with the total amount of leakage up to 6000 m3/h in the dry seasons and 12,000 m3/h in the wet seasons(Fig.36e).For successful construction of the dam,there was an urgent needtoclarifyhowthe leakage occurred(i.e.theoriginand thepathof the leakage),what treatments should be taken,and how the leakage might influence the long-term safety of the dam.To answer these questions,the key was to re-examine the permeability of the loose deposits and foundation rocks.Two techniques were taken for this purpose:onewaspumping test(Chen et al.,2016a),and the other was inverse modeling(Zhou et al.,2015a).

Fig.36.(a) Picture of the Changheba Dam site,(b) Longitudinal geological cross-section of the dam and foundation,(c) 3D finite element mesh of the numerical model,(d)Comparison of the measured and simulated discharge into the dam foundation,and (e) Comparison of the measured and simulated hydraulic heads at the piezometers.Modified from Zhou et al.(2015a).

Fig.37.Pumping tests at the Changheba dam foundation (Chen et al.,2019):(a)Layout of the foundation,(b)Layout of the pumping and observation wells,and (c)Typical crosssection (A-A) of the aquifers.

Fig.38.Time-series of drawdown from pumping tests at the Changheba dam foundation (modified from Chen et al.,2019):(a)pumping test 1,(b)pumping test 2.The early-time data can be well fitted by the 1D flow model.The insets in (a) and (b) show that the late-time flow behavior is in agreement with the 2D radial flow model.

Fig.37 shows the layout of two pumping tests performed in the foundation bounded by two cut-off walls.The drawdown data were interpreted by the GNRF model (Section 4.4.2),showing interestingly that the well flow in the foundation involved a transition of flow dimension from 1D linear flow to 2D radial flow (Chen et al.,2019),as plotted in Fig.38.This transition exactly demonstrates the underlying mechanism of the leakage.During the early time of the tests(<3000 s),only the pore water in the deposits responded to the pumping;but as the tests proceeded (>7000 s),the lateral flow through the moderately-weathered bedrock layer played an important role because when the tests were performed,the grout curtain below the main cut-off wall had not been constructed yet(Fig.37c).The interpretation of the early-time and late-time drawdown data with the 1D and 2D flow models yielded hydraulic conductivity values of (1.18-2.52) × 10-3m/s and (1.61-2.32)×10-5m/s,respectively(Chen et al.,2019).These two ranges ofKvalues well represent the permeability of loose deposits and moderately-weathered bedrocks,as validated by a large number of packer tests and pumping tests during site characterization (Zhou et al.,2015a).Therefore,the pumping tests indicated that the leakage was originated from the Dadu river and took place along the moderately-weathered bedrock layer.

Based on the time-series data of discharge and piezometric head,the weighted objective inverse modeling technique(Section 5.2)was adopted to back calculate the permeability of the foundation rocks(Zhou et al.,2015a),as shown in Fig.36.It shows that the moderately-weathered rocks are of moderate permeability and have a representativeKvalue of 1.77×10-5m/s,which is almost the same as the result interpreted from pumping tests.The results confirm that the unexpected leakage occurred in the early construction stage through the moderately-weathered rock layer that was not cut off by the temporary barriers of cofferdams for seepage control in construction(Fig.36b),with the river water as the main source of supply.The inverse modeling further confirms that the seepage-proof performance of the cut-off walls meets the design criterion,and as long as the quality of the main grout curtain was well controlled in the subsequent construction,the leakage had little effect on the operation safety of the dam (Zhou et al.,2015a).By the end of 2017,the normal pool level of the reservoir was attained,and the measured discharge through the dam foundation was about 30 L/s,which justifies the above results of the leakage assessment.

8.3.Groundwater flow in karst areas in response to deep tunnel construction

Large-scale underground engineering projects in karst areas are not only often affected by geological disasters such as water inrush but also can result in strong disturbance of the karst groundwater system with significant ecological and/or environmental consequences.The Jiayan Water Diversion Project(JWDP)is a large-scale hydraulic engineering project in Guizhou Province,China.The project implements construction of deep tunnels with a total length of 648.19 km in the karst region.The rocks of the study area are mainly limestone and dolomite.Underground river networks develop in the area.During construction of the project,serious water inrush problems were encountered at two long tunnels(Maochang and Shuidaqiao tunnels,～36 km long) and the tunnel discharge caused significant groundwater drawdown in the area.This raised important issues of(1)how to accurately predict water inrush caused by progressive excavation to ensure engineering safety and (2) how to assess the impact on the karst groundwater system.To address these questions,the key was to comprehensively model the regional groundwater flow with special consideration on the karst hydrogeological characteristics.Two approaches were adopted to tackle the issues: the lumped parameter model (Zhou et al.,2021) and the coupled discretecontinuum modeling (Zheng et al.,2021).

Fig.39.Daily discharge with 95% confidence intervals simulated by the lumped parameters model.

The lumped parameter model (Section 5.4) was first used to characterize the relationship between groundwater recharge and discharge of the Shuchang HU in the karst system.The model includes a series of water balance equations describing four interacting flow compartments.The model parameters were calibrated using the precipitation and spring discharge data.It was found that the surface runoff,storage and recession in the epikarst zone,and the quick flow in the phreatic zones are the dominant processes that control the recharge-discharge relationship in the study area(Zhou et al.,2021).It was further shown that the tunnel construction captured an average of 19%of the total discharge(Fig.39).This model application demonstrated that as a screening approximation,the lumped models can provide a useful tool for identifying the dominant processes that govern the groundwater system response and for assessing the groundwater budget changes induced by engineering activities including tunneling and underground mining in karst areas.

3D numerical modeling of groundwater flow with detailed consideration of complex heterogeneity and the large-scale tunnel construction progression in karst aquifers at a regional scale is a challenging task.The coupled discrete-continuum modeling approach was adopted to simulate the groundwater flow for JWDP with a focus on the water inrush to tunnels and the associated hydrogeological impact at the large scale.In the 3D transient numerical model,both the karst conduits and tunnels are treated as discrete channels and the tunnel lengths and properties are varying with time according to the construction schedule.After model calibration based on site characterization and observation data,high-resolution simulations were performed to quantitatively evaluate the influence of tunnel construction on the karst groundwater system (Fig.40).It was shown that the simulation results accurately reproduce the time-series of tunnel discharges.The water inflow to the tunnels during the construction period led to significant lowering of the groundwater level (Fig.40a).The surface impact area(with drawdown ≥5 m)covered over 60%of the study area after 4 years of construction (Fig.40b).The simulation results showed that after lining of the tunnels,the groundwater system will slowly recover to a groundwater level 5-10 m below the undisturbed level (Fig.40d).When there is discharge through the drainage systems in the tunnel after lining,the recovery time will increase significantly.This application demonstrates the use of 3D numerical models to simulate regional groundwater flow in karst areas with large-scale underground engineering activities and to evaluate the spatial and temporal response of the groundwater system to the disturbance produced by tunnel construction.Therefore,the simulation methodology can be effectively utilized to optimize construction plans and ensure construction safety as well as to minimize regional scale environmental impact.

Fig.40.(a)Evolution of simulated groundwater level during tunnel construction at a point along the tunnel axis,(b)Spatial distribution of tunnel-construction induced drawdown at different times from the start of the project (top) to the end of excavation (bottom),(c) Water influx to karst conduits during the construction period,and (d) Predicted groundwater level recovery at the point as in (a).Modified from Zheng et al.(2021).

9.Concluding remarks and future directions

Significant advances both in fundamental theories and engineering applications of groundwater flow in fractured rocks have accumulated over the last couple of decades.Microscopic mechanisms of fluid flow in fractured rocks,especially under the influence of non-Darcian conditions,phase interference,rock dissolution,and particle transport,have been investigated through a combination of visualized experiments and theoretical analysis.Methods of characterizing hydraulic properties of fractured rocks in different flow regimes and at the field scale have been proposed.Highperformance numerical simulation approaches for large-scale modeling of groundwater flow in fractured rocks and aquifers have been developed,which enables simulation-based optimization design of seepage control systems in various settings.Mechanisms of coupled hydro-mechanical processes and mitigation of flow-induced deformation have also been discussed.The improved theoretical understanding,field characterization methods,numerical modeling approaches,and seepage control techniques have been successfully applied to engineering practices.The major contributions of this study are summarized as follows:

(1) A high-resolution but low-cost flow visualization system was developed.Some novel flow patterns of non-Darcian flow,two-phase flow,dissolution and fluid-driven particle transport through rough fractures or fracture networks were discovered.Theoretical modes were proposed to characterize the transitions of the flow,dissolution and particle transport patterns.

(2) A series of analytical models and empirical methods was proposed to characterize the hydraulic properties of fractured rocks at the field scale,including the permeability variations induced by excavation and fracture clogging,the identification of inertial permeability based on the techniques of highpressure packer test,pumping test and statistical correlation,and the estimation of unsaturated flow properties based on statistical correlation and inverse modeling.

(3) A unified variational inequality method was proposed for modeling of steady-state,transient and saturatedunsaturated flows with complex unilateral boundary conditions of Signorini’s type.Also presented were an efficient multi-objective inverse modeling technique,a regional-scale modeling strategy,and a combined lumped-parameter and discrete-continuum modeling strategy for groundwater flow under diverse difficult field conditions.

(4) The mechanisms of seepage control in geotechnical engineering were clarified and mathematically represented.A refined modeling technique was developed for large-scale drainage systems containing tens of thousands of drains with complex boundary conditions.An optimization design procedure was proposed for seepage control systems,aiming to maximize the seepage control performance,minimizing the construction cost and mitigating the groundwater flowinduced deformation.

(5) The proposed theories,models,methods and techniques were proved efficient and effective by a great number of applications to large-scale hydraulic,hydropower and geotechnical projects.Three typical case studies were shown to illustrate the applications to optimization design,performance assessment,leakage treatment and environmental impact evaluation,etc.

Despite the above advances,new frontiers constantly emerge due to the complexities of geological environment and physical processes as well as the evolving engineering practices,which calls for further dedicated research of this important topic.An outlook of future research directions is given below:

(1) Fundamental mechanisms of flow through fractured rocks merit further consideration with special focuses on the impacts of non-Darcian and multi-phase flows,fluid rheological properties,chemical reactions,fluid-rock interactions,etc.

(2) Effective upscaling methods need to be developed to incorporate small-scale flow mechanisms into processes/properties of flow in fractured rocks at larger scales.

(3) Accuracy in characterization of hydraulic parameters for fractured rocks at the field scale can be improved by e.g.high-resolution monitoring with integrated hydrogeological and geophysical tools and technologies.

(4) High-performance,high-fidelity,versatile numerical modeling approaches of flow and coupled processes in fractured rock systems await further development.

(5) More practical applications of the developed theories and models are warranted to address emergent challenges in assessment,optimization and safety control of large-scale geotechnical engineering problems.

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.51988101,51925906 and 52122905) are gratefully acknowledged.