Seepage Characteristics and Its Control Mechanism of Rock Mass in High-Steep Slopes

In Southwest China large-scale hydropower projects, the hydraulic conductivity and fracture aperture within the rock mass of a reservoir bank slope has dramatically undergone a time series of evolution during dam abutment excavation, reservoir impounding and fluctuation operation, and discharge atomization. Accordingly, seepage control measures by hydro-structures such as drainage or water insulation curtains should be guided by scientific foundation with a dynamic process covering life-cycle performance. In this paper, the up-to-date status of studying the evolution mechanism of seepage characteristics relating to fractured rock hydraulics from experimental samples to the engineering scale of the rock mass is reviewed for the first time. Then, the experimental findings and improved practice method on nonlinear seepage flow under intensive pressure drives are introduced. Finally, the scientific progress made in fractured rock seepage control theory and optimization of the design technology of high-steep slope engineering is outlined. The undertaken studies summarized herewith are expected to contribute to laying a foundation to guide the further development of effective geophysical means and integrated monitoring systems in hydropower station construction fields.


Introduction
In Southwest China, from Sichuan to Yunnan province, large-scale hydro-power stations such as Xiaowan in Lantsang River, Jinping in Ya-lung (Nyag Chu in Tibetan) River (shown in Figure 1), Dagangshan in Dadu River, Xiluodu, Xiangjiaba, Wudongde and Baihetan successively in Chin-sha Chiang River have been constructed in the last 10 years.As shown in Table 1, taking four of them as examples, their complex rock slopes as dam abutments are high and steep.
walls on seepage flows in a fracture, the applicability of the cubic law and its modifications as well as stress sensitivity in fracture seepage [5][6][7][8].Nevertheless, in the case of unsaturated conditions, there is not yet plenty of in-depth findings on water and air two-phase flow through a rock fracture.The theoretical models on multiphase flow in fractured rocks are mostly founded on soil-water characteristic curves in porous media [9][10][11].More seriously, the study on flow regime as well as the hysteretic behavior of water and air two-phase flows through a rock fracture are rarely seen.To study the seepage characteristics and seepage law of motion in fractured rock mass, equivalent continuum and discrete network approaches have been mainly adopted.Between them, the highlight of the equivalent continuum approach is to study the scale effect and anisotropy of hydraulic conductivity and specific storage coefficient so as to establish the equivalent permeability tensor of rock mass and discover the stress sensitive evolution regularity of seepage characteristics [12][13][14].Although field testing methods represented by the cross-hole permeability test [15,16] have been made, progress to develop effective geophysical approaches such as microseismic monitoring is urgently needed.This facilitates understanding of how damage or failure occurring inside the rock mass influences its corresponding seepage characteristics while a dam abutment slope is being excavated or loaded [17,18].Many researchers [19][20][21] studied the features of permeability evolution along the process of rock deformation as well as progressive failure.Their findings embodied the intensive scaling effect of permeability.Nevertheless, for unsaturated seepage properties of fractured rock mass, mechanisms of air and water two-phase flow as well as rainfall infiltration, study progresses only slowly.
Equivalent continuum, discrete network, and dual porosity methods are usually adopted to numerically simulate the seepage field of slope rock mass.Concerning the steady and transient seepage flow related problems, in order to solve the numerical stability problem caused by nonlinearity of flux spillage boundary, formulations of elliptical and parabolic Signorini's type variational inequality have been successively established [22].Accordingly, incorporated with a drainage substructure method to numerically simulate its seepage flow, a precise simulation technique was proposed for analyzing large-scale drainage curtain structures [23].Nowadays, the seepage control technique for slope rock mass has developed from the early stage of leakage prevention to the stage of taking drainage to be as important as leakproof.The rationality and  Slope stability is a major issue during design, construction.and operation of hydraulic and hydro-power engineering [1].It is justified by the large quantity of research that the vast majority of landslides are related to water [2].Distinguishing rock hydraulics [3] from classical poromechanics was recognized by humans by learning tragic lessons from world-shaking catastrophes such as the French Malpasset arch dam break in 1959 and the Italian Vajont tremendous left bank landslide in 1963 [4].Hydraulic characteristics of rock discontinuities and the regularity of groundwater movement are the foundation and key to implement slope seepage analysis.In the case of saturated conditions, theoretical and experimental findings on the seepage in a single rock fracture are abundant.They reflect the influences of geometric configurations and stretch features of fracture walls on seepage flows in a fracture, the applicability of the cubic law and its modifications as well as stress sensitivity in fracture seepage [5][6][7][8].Nevertheless, in the case of unsaturated conditions, there is not yet plenty of in-depth findings on water and air two-phase flow through a rock fracture.The theoretical models on multiphase flow in fractured rocks are mostly founded on soil-water characteristic curves in porous media [9][10][11].More seriously, the study on flow regime as well as the hysteretic behavior of water and air two-phase flows through a rock fracture are rarely seen.
To study the seepage characteristics and seepage law of motion in fractured rock mass, equivalent continuum and discrete network approaches have been mainly adopted.Between them, the highlight of the equivalent continuum approach is to study the scale effect and anisotropy of hydraulic conductivity and specific storage coefficient so as to establish the equivalent permeability tensor of rock mass and discover the stress sensitive evolution regularity of seepage characteristics [12][13][14].Although field testing methods represented by the cross-hole permeability test [15,16] have been made, progress to develop effective geophysical approaches such as microseismic monitoring is urgently needed.This facilitates understanding of how damage or failure occurring inside the rock mass influences its corresponding seepage characteristics while a dam abutment slope is being excavated or loaded [17,18].Many researchers [19][20][21] studied the features of permeability evolution along the process of rock deformation as well as progressive failure.Their findings embodied the intensive scaling effect of permeability.Nevertheless, for unsaturated seepage properties of fractured rock mass, mechanisms of air and water two-phase flow as well as rainfall infiltration, study progresses only slowly.
Equivalent continuum, discrete network, and dual porosity methods are usually adopted to numerically simulate the seepage field of slope rock mass.Concerning the steady and transient seepage flow related problems, in order to solve the numerical stability problem caused by nonlinearity of flux spillage boundary, formulations of elliptical and parabolic Signorini's type variational inequality have been successively established [22].Accordingly, incorporated with a drainage substructure method to numerically simulate its seepage flow, a precise simulation technique was proposed for analyzing large-scale drainage curtain structures [23].Nowadays, the seepage control technique for slope rock mass has developed from the early stage of leakage prevention to the stage of taking drainage to be as important as leakproof.The rationality and effectiveness of slope seepage control design depends on not only correctly understanding the seepage control mechanism but also scientifically evaluating seepage control results.Fortunately, a number of researchers' achievements paved the way to evaluate the effect of seepage control in high-steep slopes [24][25][26][27][28].It has been summarized in recent years that the mechanism comes down to process, state, parameter, and boundary, i.e., the four categories of rock mass seepage control [23].Thus, laying a theoretical cornerstone for seepage control optimization design.
Therefore, it is badly needed in an intensive way to study the topics: evolution of seepage characteristics for intensively un-loading slope rock mass, unsaturated seepage and water-air two-phase flow motion features, rainfall infiltration mechanism as well as seepage control optimization design.
Encompassing three key aspects of the titled problem, i.e., seepage related properties, seepage motion regularities, and seepage control, it is the aim to systematically summarize the following latest findings in this review.These profound studies concern not only the macro-and meso-scopic mechanism, multiple scaling effect and anisotropic feature on the process evolution of the seepage characteristics of high-steep slope fractured rocks, but also the model solution and value selection of non-linear seepage flow parameters, as well as multi-purpose whole process dynamic inverse analysis of rock mass seepage fields.

Evolution Mechanism on Seepage Characteristics of High-Steep Slope Rock Mass
Natural evolution of deep valleys and large scale excavation of high-steep slopes both cause a strong unloading impact on the slope rock mass.Extensively developed discontinuous structures as well as rock damage play a dominant role in seepage flow properties.The following parts contain three aspects of the distinct findings: (1) the rapid laboratory technique on the permeability testing for low-permeable rocks and evolution models; (2) evolution mechanism and its scaling effect on the seepage characteristics of fractured rocks; (3) the multipurpose and whole process dynamic inverse analysis method on the seepage characteristics of fractured rock mass.They have been studied to deal with the challenge of predicting the performance of instantaneous seepage characteristics in the rock mass.
Utilizing a TAW-1000 deep water pore pressure tri-axial servo rock testing machine and a French tri-axial seepage rheological experimental system where the water pressure attains 70 MPa, the association between the seepage and deformation properties of fractured rock was studied [29].Furthermore, acoustic emission signals were recorded from the whole stress-strain process of rock deformation and failure [30].The inherent connection between deformation, permeability, and acoustic emission data was setup.In the meantime, sonic wave and radial testing using thick-wall cylindrical ring rock sample [31] was also attempted and developed, respectively.As well, simulators of various kinds of numerical methods such as Finite Element Method, Boundary Element Method, Distinct Element Method, Control Volume Finite Difference Method, and Particle Flow Code were adopted to analyze the experiments accordingly in both continuum and discontinuous, macro-and mesoscopic, parallel and large-scale high performance scientific computational ways.These studies obtained a large amount of fundamental systematic data sets on the relationship between hydraulic conductivities and the whole process of stress-strain behaviors, especially in (red) sandstone, limestone, granite, and marble.The fact that modelling computations displayed matching results with experimental data laid a solid foundation to further understand the evolution of seepage characteristics along with various loading and deformation states.
Blocked by high cost and long cycle, the existing rapid testing way of rock low permeability has become a technical difficulty.Through integrating an instantaneous pressure pulse device into the Temperature-Hydraulic-Mechanical coupling testing facility, shown in Figure 2, an experimental technique coupling water-air two-phase seepage flows and triaxial compression of low permeable rock samples was researched and developed.This system contains servo measuring of the axial, transverse surrounding and air-water partial pore pressures as well as automatic data acquisition and recording.As its range reaches 10 −12 -10 −22 m 2 wide with internationally equivalent level of precision, it was successfully applied in the Canadian URL (Underground Research Laboratory) deep underground laboratory [32,33].Through analyzing sonic waveforms which penetrate different categories of rock samples under various conditions of moisture and saturation, the association between the features in time, frequency, time-frequency domains of acoustic signals and the pore and microcrack development, water-bearing condition of rocks was explained [35].Through modeling acoustic emissions and energy dissipation of particle flows in rocks under water-mechanical coupling conditions, it was concluded that tensile cracks are always dominant under low confining pressure and shearing cracks are dominant when strain increases under high confining pressure [36].Under an environment of particle flow discrete element software PFC2D, a fluid-solid coupling numerical simulation analysis of the particle flow discrete element is carried out by using embedded FISH (An imbedded programming language to enable users to define variables and functions).The influence of porosity dynamic distribution on permeability evolution characteristics under different osmotic pressures was explained by the micro-mechanism [37], as illustrated in Figure 3.The stress-strain curves and acoustic emission data of sandstone were obtained through a hydraulic coupling triaxial test and realtime Acoustic Emission monitoring.In the meantime, a hydraulic coupling bi-axial model of discrete element and particle flow was established to study the laws of acoustic emission and energy evolution in sandstone.Through analyzing sonic waveforms which penetrate different categories of rock samples under various conditions of moisture and saturation, the association between the features in time, frequency, time-frequency domains of acoustic signals and the pore and microcrack development, water-bearing condition of rocks was explained [35].Through modeling acoustic emissions and energy dissipation of particle flows in rocks under water-mechanical coupling conditions, it was concluded that tensile cracks are always dominant under low confining pressure and shearing cracks are dominant when strain increases under high confining pressure [36].Under an environment of particle flow discrete element software PFC2D, a fluid-solid coupling numerical simulation analysis of the particle flow discrete element is carried out by using embedded FISH (An imbedded programming language to enable users to define variables and functions).The influence of porosity dynamic distribution on permeability evolution characteristics under different osmotic pressures was explained by the micro-mechanism [37], as illustrated in Figure 3.The stress-strain curves and acoustic emission data of sandstone were obtained through a hydraulic coupling triaxial test and real-time Acoustic Emission monitoring.In the meantime, a hydraulic coupling bi-axial model of discrete element and particle flow was established to study the laws of acoustic emission and energy evolution in sandstone.
dynamic distribution on permeability evolution characteristics under different osmotic pressures was explained by the micro-mechanism [37], as illustrated in Figure 3.The stress-strain curves and acoustic emission data of sandstone were obtained through a hydraulic coupling triaxial test and realtime Acoustic Emission monitoring.In the meantime, a hydraulic coupling bi-axial model of discrete element and particle flow was established to study the laws of acoustic emission and energy evolution in sandstone.Coupled Hydraulic-Mechanical-Damage model and efficient solution of large-scale Finite Element equations for fully coupled problems were also studied [38].The numerical solution to the two-part Hooke model (TPHM) was realized in TOUGH-FLAC3D.The conditions of its applicability were also defined.This method essentially reflects the nonlinear deforming behavior during the course of loading/unloading in a low stress stage so as to conform to the excavation unloading induced deforming features in the surrounding rocks.This model is more practical to design and analyze deformation characters of surrounding rocks in engineering with unloading disturbance [39].Through the numerical simulation using the statistical meso-damage finite element method [40], as shown in Figure 4 [41], the effects of inhomogeneity, confining pressure, and joint density on the size effects of mechanical and seepage characteristics were discussed.A quantitative estimation method was proposed for the mechanical and hydraulic properties of rock mass with meso-to macro-levels of multi-scale fractures through numerous systematic simulations.The parameters of deformation, strength, fracture porosity, and permeability of the rock mass were represented in the numerical models.Although computational resources and cost are in most cases rather expensive and the applicability of cubic law as well as its restriction needs further validation, it can provide a reference for reasonably acquiring and extrapolating computational engineering parameters for practical highsteep slope rock mass.Coupled Hydraulic-Mechanical-Damage model and efficient solution of large-scale Finite Element equations for fully coupled problems were also studied [38].The numerical solution to the two-part Hooke model (TPHM) was realized in TOUGH-FLAC3D.The conditions of its applicability were also defined.This method essentially reflects the nonlinear deforming behavior during the course of loading/unloading in a low stress stage so as to conform to the excavation unloading induced deforming features in the surrounding rocks.This model is more practical to design and analyze deformation characters of surrounding rocks in engineering with unloading disturbance [39].Through the numerical simulation using the statistical meso-damage finite element method [40], as shown in Figure 4 [41], the effects of inhomogeneity, confining pressure, and joint density on the size effects of mechanical and seepage characteristics were discussed.A quantitative estimation method was proposed for the mechanical and hydraulic properties of rock mass with meso-to macro-levels of multi-scale fractures through numerous systematic simulations.The parameters of deformation, strength, fracture porosity, and permeability of the rock mass were represented in the numerical models.Although computational resources and cost are in most cases rather expensive and the applicability of cubic law as well as its restriction needs further validation, it can provide a reference for reasonably acquiring and extrapolating computational engineering parameters for practical high-steep slope rock mass.
strength, fracture porosity, and permeability of the rock mass were represented in the numerical models.Although computational resources and cost are in most cases rather expensive and the applicability of cubic law as well as its restriction needs further validation, it can provide a reference for reasonably acquiring and extrapolating computational engineering parameters for practical highsteep slope rock mass.In a recent study, an energy conservative microscopic and anisotropic damage mechanics model was established for brittle rocks [34].The aforementioned study considered the normal stiffness recovery and sliding dilatancy of microcracks as follows [34]: where, E is the macropscopic strain tensor; S s is the flexibility for rock matrix; Σ is the macroscopic stress tensor of inhomogeneous rock; β, γ are internal variables characterizing crack opening and sliding, respectively; l 2 is the unit spherical area,  represents the symmetrical tensor formed by cross multiplication from two vectors.Also, considering crack initiation, propagation, changes in crack morphology, and volume fraction as well as the micro-structural feature of crack connectivity, a new damage-induced permeability evolution model was established [34].This model in addition characterizes crack connectivity as an important microscopic structural feature as follows [34]: In a recent study, an energy conservative microscopic and anisotropic damage mechanics model was established for brittle rocks [34].The aforementioned study considered the normal stiffness recovery and sliding dilatancy of microcracks as follows [34]: where, E is the macropscopic strain tensor; S s is the flexibility for rock matrix; Σ is the macroscopic stress tensor of inhomogeneous rock; β, γ are internal variables characterizing crack opening and sliding, respectively; l 2 is the unit spherical area, ⊗ represents the symmetrical tensor formed by cross multiplication from two vectors.Also, considering crack initiation, propagation, changes in crack morphology, and volume fraction as well as the micro-structural feature of crack connectivity, a new damage-induced permeability evolution model was established [34].This model in addition characterizes crack connectivity as an important microscopic structural feature as follows [34]: where, ϕ c is the density function of cracks; n is the unit normal vector of a crack; δ is the two-order unit tensor; k s is the permeability tensor of rock matrix; k 0 c is initial permeability of microcracks; β 0 , β are the initial and current internal variables of microcracks; d 0 , d are the initial and current density of crack distribution, respectively.It was proved that the traditional multi-scale homogenization based model has a lower bound estimate nature and overcomes the limitation due to lack of describing the microcracks' connectivity.The number of parameters is less by 1 to 4 that of other models [42,43].
Compared with experimental tests on Swedish Äspö diorite, Lac du Bonnet, Senones and Beishan granite, the correctness of the above models was verified due to the consistent agreement between prediction and measurement [44].Anisotropic damage of rock, structural plane dilatancy, and wear as well as statistical characteristics of fracture development is taken into account in this model.It reveals the macro-meso mechanism, multi-scale effect, and anisotropic characteristics of permeability evolution in rock mass, which was applied to study the excavation disturbance effect and permeability evolution law for Jinping I hydropower station left-bank slope.
Based on a multi-scale homogenization approach, the model of elasto-plastic damage and seepage characteristics evolution for fractured rock mass was established [45].Among the scales of rock, structural plane and rock mass, macro-/meso-scopic mechanism, the multiscale effect and anisotropy in rock mass seepage characteristics evolution was revealed.The above mentioned models count not only the dominant function of fracture development features on rock mass permeability but also the contribution to rock mass seepage characteristics from rock damage caused by intensive disturbance and post-peak deformation behavior of structural planes.The advantage brought by converging between the rock mass structure model and the damage model created a condition for coupling mechanical and hydraulic analysis on rocks during implementing large-scale excavation.The findings were successfully applied in studying the influence of excavation relaxation on seepage characteristics evolution as well as field seepage control effects of the high-steep slopes and underground powerhouse group in Jinping I and Kala hydropower stations [46,47].

Nonlinearity in Seepage Flow of High-Steep Slope Rock Mass and Its Analysis
A series of theories and control technologies on seepage flow analysis under the condition of high osmotic pressure was established.It not only manifests the seepage motion characters in slopes under rainfall or atomization conditions, but also reveals the unsaturation, multi-phase coupling, and non-Darcy flow features of seepage flows in deformable rough-wall fractured rocks.
First, experiment and modeling on the seepage mechanism of a rough fracture under high osmotic pressure were undertaken.As shown in Figure 5, using high-precision 3D laser scanning, an accurate technology was proposed to reconstruct the initial aperture distribution of a structural plane.
A series of theories and control technologies on seepage flow analysis under the condition of high osmotic pressure was established.It not only manifests the seepage motion characters in slopes under rainfall or atomization conditions, but also reveals the unsaturation, multi-phase coupling, and non-Darcy flow features of seepage flows in deformable rough-wall fractured rocks.
First, experiment and modeling on the seepage mechanism of a rough fracture under high osmotic pressure were undertaken.As shown in Figure 5, using high-precision 3D laser scanning, an accurate technology was proposed to reconstruct the initial aperture distribution of a structural plane.In order to study the seepage characteristics of rough structural planes in a visual environment, a multi-phase percolating flow testing device which is capable of replicating a transparent rough crack by epoxy resin material was developed.At the same time, as shown in Figure 6, it provides an important means for visual observation of coarse-fracture multiphase flows, solute migration, and viscous fingering [49,50].By laboratory experiments, the nonlinear coefficient for Forchheimer's law on the non-Darcy flow relationship between flow rate Q and pressure gradient P of a deformable rough fracture was determined.Furthermore, as shown in Figure 7  In order to study the seepage characteristics of rough structural planes in a visual environment, a multi-phase percolating flow testing device which is capable of replicating a transparent rough crack by epoxy resin material was developed.At the same time, as shown in Figure 6, it provides an important means for visual observation of coarse-fracture multiphase flows, solute migration, and viscous fingering [49,50].By laboratory experiments, the nonlinear coefficient for Forchheimer's law on the non-Darcy flow relationship between flow rate Q and pressure gradient P of a deformable rough fracture was determined.Furthermore, as shown in Figure 7  Second, an analytical method for nonlinear seepage parameters of fractured rock mass was established.Based on a borehole high pressure test, aided by borehole video recording and sonic detection, it yielded a P-Q curve type, its feature and the intrinsic relationship between water injection pressure, rock mass integrity, and geostress level.An analytical model for high water pressure testing under nonlinear flow conditions was given.The model overcomes the disadvantage that the permeability coefficient obviously turns out to be smaller when conventional pressure test Second, an analytical method for nonlinear seepage parameters of fractured rock mass was established.Based on a borehole high pressure test, aided by borehole video recording and sonic wave detection, it yielded a P-Q curve type, its feature and the intrinsic relationship between water injection pressure, rock mass integrity, and geostress level.An analytical model for high water pressure testing under nonlinear flow conditions was given.The model overcomes the disadvantage that the permeability coefficient obviously turns out to be smaller when conventional pressure test rules are applied to a high water pressure test which brings a huge safety risk to the engineering anti-percolation and drainage design.The test showed that the P-Q curve can be divided into three typical stages: Darcy flow segment (I), non-Darcy flow segment (II), and hydraulic splitting segment (III).As shown in Figure 8, it can also be divided into two basic types, A and B: Curve A-type only contains two sections I and II.The pressure rising curve and falling curve are basically coincident.The fracture irreversible deformation is not distinct; Curve B-type includes three sections I, II, and III.The rising and falling curves do not coincide and irreversible deformation or hydraulic cleavage occurs to the fracture [58].Under the spherical diffusion flow assumption, an analytical model for nonlinear seepage parameters of rock mass based on nonlinear Izbash law and Forchheimer's law (−∇P = AQ + BQ2) under high pressure condition was established.The method to estimate the permeability coefficient of the rock mass in three different stages, namely linear, nonlinear, and hydraulic fracturing sections, was recommended [59,60].It should be noted that the above model would soon be included in the revising test code after going through successful verifications in the field application.
fracture [58].Under the spherical diffusion flow assumption, an analytical model for nonlinear seepage parameters of rock mass based on nonlinear Izbash law and Forchheimer's law (−P = AQ + BQ2) under high pressure condition was established.The method to estimate the permeability coefficient of the rock mass in three different stages, namely linear, nonlinear, and hydraulic fracturing sections, was recommended [59,60].It should noted that the above model would soon be included in the revising test code after going through successful verifications in the field application.Third, as shown in Figure 9, a large-scale physical experiment was carried out to study the water coupling deformation features of a slope under strong atomization.The model scale is 8000 × 2000 × 4000 mm 3 , which realized water level fluctuation, simulated groundwater head varying in rock layers, and studied the evolution process of slope displacements under reservoir water fluctuation and rainfall conditions.The dip angle of the slope was adjusted to study the ultimately stable dip angle under different saturation.In addition, a contactable monitor system and a non-contactable digital speckle measurement system realized on-line continuous deformation measurement [61].On the other hand, the model for the soil-water characteristic curve containing the deformation and hysteresis effect was established.It revealed the coupling mechanism between nonlinear deformation of soil swelling, its collapse by moisture, and unsaturated seepage flow.Based on these, as shown in Figure 10, the water-holding characteristic curve of deformable rock mass and the calculation method of unsaturated hydraulic conductivity were formed [62,63].According to the proposed soil water characteristic curve model, a modified Mualem statistical model for the relative permeability coefficient of unsaturated soil, considering soil deformation effect, was established.The comparison with the experimental results showed that the model can estimate the relative permeability coefficient of the soil under deforming conditions, whose prediction ability is stronger than the Mualem statistical model.It provides a more reliable theoretical model for studying the seepage law of unsaturated soil and its numerical simulation of H-M coupling.In addition, according to the seepage characteristics of the slope rock mass in the atomized area, a coupled model of slope water-air twophase flow motion considering rock mass deformation was established to explain the hysteretic effect of the wet front advancing process, gas migration, rock mass deformation, and drainage on the infiltration evolution process under heavy rainfall [64].Third, as shown in Figure 9, a large-scale physical experiment was carried out to study the water coupling deformation features of a slope under strong atomization.The model scale is 8000 × 2000 × 4000 mm 3 , which realized water level fluctuation, simulated groundwater head varying in rock layers, and studied the evolution process of slope displacements under reservoir water fluctuation and rainfall conditions.The dip angle of the slope was adjusted to study the ultimately stable dip angle under different saturation.In addition, a contactable monitor system and a non-contactable digital speckle measurement system realized on-line continuous deformation measurement [61].On the other hand, the model for the soil-water characteristic curve containing the deformation and hysteresis effect was established.It revealed the coupling mechanism between nonlinear deformation of soil swelling, its collapse by moisture, and unsaturated seepage flow.Based on these, as shown in Figure 10, the water-holding characteristic curve of deformable rock mass and the calculation method of unsaturated hydraulic conductivity were formed [62,63].According to the proposed soil water characteristic curve model, a modified Mualem statistical model for the relative permeability coefficient of unsaturated soil, considering soil deformation effect, was established.The comparison with the experimental results showed that the model can estimate the relative permeability coefficient of the soil under deforming conditions, whose prediction ability is stronger than the Mualem statistical model.It provides a more reliable theoretical model for studying the seepage law of unsaturated soil and its numerical simulation of H-M coupling.In addition, according to the seepage characteristics of the slope rock mass in the atomized area, a coupled model of slope water-air two-phase flow motion considering rock mass deformation was established to explain the hysteretic effect of the wet front advancing process, gas migration, rock mass deformation, and drainage on the infiltration evolution process under heavy rainfall [64].Fourth, the thermo-hydro-mechanical coupling mechanism and modeling of unsaturated soils were established.A large number of studies have shown that the adhesion effect between soil particles in the unsaturated state has a significant impact on its deformation.Pore distribution as well as its evolution inside soil, play a decisive role on the soil-water characteristic curve.Based on the soil meso-structure and its evolution characters, an explicitly expressed factor which has definite mechanical meaning in characterizing bonding between unsaturated soil particles was proposed.Under the framework of the modified Cambridge model, an elastoplastic constitutive model of unsaturated soil was setup with a homogenized skeleton stress and particle bonding factor as basic variables.On the basis of the soil pore distributive evolution law during the deformation process, a soil-water characteristic curve model considering the deformation and hysteresis effect was expressed.As shown in Figure 11, the model of unsaturated seepage characteristics, the full-coupling model of elastoplastic deformation, and the capillary hysteresis of unsaturated soil were established under a thermodynamic framework, revealing the coupling mechanism of seepage flow on nonlinear deformation and unsaturation.Compared with the classic Barcelona (BBM) model in the field of unsaturated soil mechanics [65], the above fully coupled hydraulic model for unsaturated soil not only in a more fundamental way reflects two-way coupling characteristics of elastoplastic deformation and capillary hysteresis under the thermodynamic framework, but also has a more direct building method and clearer physical meaning.Only a single yield surface (BBMs has two) and fewer parameters (four less than BBM) are adopted.A mechanism of shear zone formation during a rain-induced landslide was revealed, laying a foundation for precisely simulating the hydraulic coupling process of the 300 m level of rockfill dam wall.Fourth, the thermo-hydro-mechanical coupling mechanism and modeling of unsaturated soils were established.A large number of studies have shown that the adhesion effect between soil particles in the unsaturated state has a significant impact on its deformation.Pore distribution as well as its evolution inside soil, play a decisive role on the soil-water characteristic curve.Based on the soil meso-structure and its evolution characters, an explicitly expressed factor which has definite mechanical meaning in characterizing bonding between unsaturated soil particles was proposed.Under the framework of the modified Cambridge model, an elastoplastic constitutive model of unsaturated soil was setup with a homogenized skeleton stress and particle bonding factor as basic variables.On the basis of the soil pore distributive evolution law during the deformation process, a soil-water characteristic curve model considering the deformation and hysteresis effect was expressed.As shown in Figure 11, the model of unsaturated seepage characteristics, the full-coupling model of elastoplastic deformation, and the capillary hysteresis of unsaturated soil were established under a thermodynamic framework, revealing the coupling mechanism of seepage flow on nonlinear deformation and unsaturation.Compared with the classic Barcelona (BBM) model in the field of unsaturated soil mechanics [65], the above fully coupled hydraulic model for unsaturated soil not only in a more fundamental way reflects two-way coupling characteristics of elastoplastic deformation and capillary hysteresis under the thermodynamic framework, but also has a more direct building method and clearer physical meaning.Only a single yield surface (BBMs has two) and fewer parameters (four less than BBM) are adopted.A mechanism of shear zone formation during a rain-induced landslide was revealed, laying a foundation for precisely simulating the hydraulic coupling process of the 300 m level of rockfill dam wall.
Lastly, the theory and method for seepage analysis of high-steep slope engineering were established.The seepage flow problem under high water head conditions involves strong physical and boundary nonlinearities.Conventional analysis methods generally have defects such as poor convergence and numerical instability.They are difficult to adapt to complex large engineering problems.The Signorini type unified complementary form for complex seepage flow boundary conditions such as infiltration, evaporation, and overflow were proposed.The Signorini parabolic variational inequality analysis method based on both continuum media and fracture network discrete models was established.The theory of seepage flow analysis in high-steep slope engineering was extended: from stationary/transit to saturated/unsaturated analyses; from linear Darcy to nonlinear non-Darcy flow analysis; from continuous medium to fracture network seepage flow analyses.In general, a high-efficiency numerical simulation method on the complex seepage flow process of high-steep slopes in hydropower projects with high water heads was solved under a unified framework.Lastly, the theory and method for seepage analysis of high-steep slope engineering were established.The seepage flow problem under high water head conditions involves strong physical and boundary nonlinearities.Conventional analysis methods generally have defects such as poor convergence and numerical instability.They are difficult to adapt to complex large engineering problems.The Signorini type unified complementary form for complex seepage flow boundary conditions such as infiltration, evaporation, and overflow were proposed.The Signorini parabolic variational inequality analysis method based on both continuum media and fracture network discrete models was established.The theory of seepage flow analysis in high-steep slope engineering was extended: from stationary/transit to saturated/unsaturated analyses; from linear Darcy to nonlinear non-Darcy flow analysis; from continuous medium to fracture network seepage flow analyses.In general, a high-efficiency numerical simulation method on the complex seepage flow process of highsteep slopes in hydropower projects with high water heads was solved under a unified framework.

Theory and Technology of Seepage Control in High-Steep Slope Engineering
This section mainly explains how to realize feedback of the seepage anomaly, the fine simulation of the seepage control structure, and the dynamic optimization design of leak-proof and drainage capability.
The standards of seepage control design in hydropower projects are mainly determined by hydraulic conductivity (q) of bedrock.Nevertheless, its hydrogeological significance and value acquisition under high water head conditions are not yet well understood.Thus a seepage control design criterion was proposed based on the maximum permeability of the rock mass under various pressures in a high-pressure packer permeability test.Based on the analytical model for non-linear seepage parameters of rock mass, this criterion has a rigorous theoretical foundation.It fits for laminar flow, turbulence, expansion, erosion, filling, and other P-Q curve types.Having been applied to the seepage control optimization design of the high-pressure diversion system in the Qiongzhong Pumped Storage Power Station, it was furthermore adopted by the energy industry standard 'The Rules for Drilling Packer Permeability Test for Hydropower Engineering' in China.A seepage control theory divided into medium property, boundary condition, initial state, and coupling process was proposed.Founded on seepage motion models, the theory broke through traditional empirical design methods such as engineering analogy, scheme comparison, and selection by quantifying leak-proof, drainage and filtration mechanisms.From the above, as shown in Figure 12, a leak-proof and drainage optimization design method on the basis of hydrogeological conditions, precise simulation of seepage control structure, and dynamic feedback was established.In the project design stage, it is beneficial to layout and optimize the seepage control system.In construction and initial impounding

Theory and Technology of Seepage Control in High-Steep Slope Engineering
This section mainly explains how to realize feedback of the seepage anomaly, the fine simulation of the seepage control structure, and the dynamic optimization design of leak-proof and drainage capability.
The standards of seepage control design in hydropower projects are mainly determined by hydraulic conductivity (q) of bedrock.Nevertheless, its hydrogeological significance and value acquisition under high water head conditions are not yet well understood.Thus a seepage control design criterion was proposed based on the maximum permeability of the rock mass under various pressures in a high-pressure packer permeability test.Based on the analytical model for non-linear seepage parameters of rock mass, this criterion has a rigorous theoretical foundation.It fits for laminar flow, turbulence, expansion, erosion, filling, and other P-Q curve types.Having been applied to the seepage control optimization design of the high-pressure diversion system in the Qiongzhong Pumped Storage Power Station, it was furthermore adopted by the energy industry standard 'The Rules for Drilling Packer Permeability Test for Hydropower Engineering' in China.A seepage control theory divided into medium property, boundary condition, initial state, and coupling process was proposed.Founded on seepage motion models, the theory broke through traditional empirical design methods such as engineering analogy, scheme comparison, and selection by quantifying leak-proof, drainage and filtration mechanisms.From the above, as shown in Figure 12, a leak-proof and drainage optimization design method on the basis of hydrogeological conditions, precise simulation of seepage control structure, and dynamic feedback was established.In the project design stage, it is beneficial to layout and optimize the seepage control system.In construction and initial impounding stages, it enables feeding abnormal seepage flow back and carrying out treatment.In the long-term operation stage, it realized seepage safety evaluation and regulation.
stages, it enables feeding abnormal seepage flow back and carrying out treatment.In the long-term operation stage, it realized seepage safety evaluation and regulation.As shown in Figure 13 [59,60], a multi-objective and whole-process dynamic inverse solution for analyzing seepage flow field was proposed based on orthogonal design, forwarding analysis of transient seepage, the BP neural network, and the non-inferior sequencing genetic algorithm.As to this approach: The fine precise simulation of the large drainage hole-curtain system and Signorinitype parabolic variational inequality method for transit seepage analysis are taken as the mean; the dynamic evolution of permeability and hydrogeological boundary conditions are taken as the parameters to be inversed; The optimal approximation of measured time series of data as seepage pressure, water level and flow rate are taken as the objectives as follows [69], where K,   are the parameter sets of permeability coefficient and water-level boundary conditions to be inverted; Qi, Qi m are the calculated and measured time series of data for flow amount at measuring point i, respectively; i, i m are the calculated and measured time series of data of osmotic pressure at measuring point i, respectively; M and N are the numbers of measuring points for osmotic pressure and flow amount, respectively.This method belongs to the positive feedback analysis method, which embodies the characteristics of multi-objective optimization and whole-process inversion.It realizes the uniqueness and reliability of inverse solution to seepage field and has good application to large-scale hydropower projects.
Based on multi-field coupling theory and precise and fine simulating technology for the engineering seepage control system, CoupledHM3D, a platform for both hydro-mechanical coupling As shown in Figure 13 [59,60], a multi-objective and whole-process dynamic inverse solution for analyzing seepage flow field was proposed based on orthogonal design, forwarding analysis of transient seepage, the BP neural network, and the non-inferior sequencing genetic algorithm.As to this approach: The fine precise simulation of the large drainage hole-curtain system and Signorini-type parabolic variational inequality method for transit seepage analysis are taken as the mean; the dynamic evolution of permeability and hydrogeological boundary conditions are taken as the parameters to be inversed; The optimal approximation of measured time series of data such as seepage pressure, water level and flow rate are taken as the objectives as follows [68], where K, φ L are the parameter sets of permeability coefficient and water-level boundary conditions to be inverted; Q i , Q i m are the calculated and measured time series of data for flow amount at measuring point i, respectively; φ i , φ i m are the calculated and measured time series of data of osmotic pressure at measuring point i, respectively; M and N are the numbers of measuring points for osmotic pressure and flow amount, respectively.
impoundment period.The Jinping first-class double-curvature arch dam on Yalong River has a maximum dam height of 305 m (shown in Figure 1) and complex geological conditions in the site area.Since the second stage of impoundment in June 2013, the water inflow of the drainage gallery in the base layer of the left bank (at elevation of 1595 m) had been large.The total water inflow exceeded 30 L/s and the water pressure inside the holes reached 0.3 MPa after closure of drainage holes at Pile No. 266 m.Based on analyzing geology, water quality, monitoring and testing data, by means of multiobjective and whole-process inversion analysis of seepage field, the dynamic inversion analysis of the seepage tensor and hydrogeological boundary conditions of the rock mass at the left bank in Jinping was carried out in stages.As illustrated in Figure 14, three important conclusions were drawn as follows: (1) It was proved that water in the reservoir is the source of water inflow from the dam foundation; (2) The leakage passage of the gushing water was identified to be T2(6-8) 2-3z, i.e., namely two groups of predominant steeply dipping fractures along the river flow direction developed in marble.The geological origin of leakage passage and the mechanism of water inflow were also revealed.(3) It demonstrated the effectiveness of curtains, revealing the permeability characteristics as well as their strong anisotropy in the rock mass.A treatment measure was put forward to directly close and seal the drainage holes which fall in part with large water inflows.This method belongs to the positive feedback analysis method, which embodies the characteristics of multi-objective optimization and whole-process inversion.It realizes the uniqueness and reliability of inverse solution to seepage field and has good application to large-scale hydropower projects.
Based on multi-field coupling theory and precise and fine simulating technology for the engineering seepage control system, CoupledHM3D, a platform for both hydro-mechanical coupling analysis and seepage control optimization design, was developed by integrating the multi-scale model of seepage characteristics evolution, the unsaturated soil-water-mechanics fully coupling model, and the non-Darcy flow model of Geotechnical media.Considering seepage control, drainage, and reverse filtration, a global optimization design method for seepage control of water conservancy and hydropower engineering was developed.Thus a complete seepage control system including medium property control, boundary condition control, initial state control, and multi-field coupling control was realized.The inverse analysis system has been successfully applied to dynamically analyze water gushing of a gallery in the dam foundation at the left bank in the Jinping hydropower project during the impoundment period.The Jinping first-class double-curvature arch dam on Yalong River has a maximum dam height of 305 m (shown in Figure 1) and complex geological conditions in the site area.Since the second stage of impoundment in June 2013, the water inflow of the drainage gallery in the base layer of the left bank (at elevation of 1595 m) had been large.The total water inflow exceeded 30 L/s and the water pressure inside the holes reached 0.3 MPa after closure of drainage holes at Pile No. 266 m.
Based on analyzing geology, water quality, monitoring and testing data, by means of multi-objective and whole-process inversion analysis of seepage field, the dynamic inversion analysis of the seepage tensor and hydrogeological boundary conditions of the rock mass at the left bank in Jinping was carried out in stages.As illustrated in Figure 14, three important conclusions were drawn as follows: (1) It was proved that water in the reservoir is the source of water inflow from the dam foundation; (2) The leakage passage of the gushing water was identified to be T2(6-8) 2-3z, i.e., namely two groups of predominant steeply dipping fractures along the river flow direction developed in marble.The geological origin of leakage passage and the mechanism of water inflow were also revealed.(3) It demonstrated the effectiveness of curtains, revealing the permeability characteristics as well as their strong anisotropy in the rock mass.A treatment measure was put forward to directly close and seal the drainage holes which fall in part with large water inflows.

Discussion and Conclusions
The up-to-date studies were reviewed on three main aspects from safety perspectives in Southwest China large-scale hydropower projects.They are the dynamic evolution of hydraulic properties of the slope rock mass fractured in different scales, regularities of non-Darcy flows under extreme pressure differences in lofty rock slopes, and field hydraulic controls based on systematical and feedback modelling analyses.
In this paper, the experimental tests on seepage characteristics of both rock samples and roughwalled fractures as well as the corresponding simulation analyses were summarized in a relatively systematic way.Not only the damage change defined by crack growth but also the evolution of seepage characteristics were predicted using the models established by researchers.Furthermore, the advances in studying the non-linear regularities of high-speed seepage flows under high pressure gradient and the unsaturated water-holding characteristics within deformable rock mass are introduced for the sake of broader application to hydraulic motion problems in high-steep slopes.Based on these latest findings, the significant achievements are also concisely elaborated in promoting seepage control theory such as a complete system covering medium property, boundary condition, initial state, and multi-field coupling as well as seepage control practice through inverse modeling analysis during construction of large-scale hydropower stations.
The hydraulic characteristics of fractured hard rocks of a reservoir bank slope undergo progressive evolution in the course of dam abutment excavation, reservoir impounding and fluctuation operation, and discharge atomization in hydropower engineering.Accordingly, seepage control measures by hydro-structures such as drainage or water insulation curtains should be guided

Discussion and Conclusions
The up-to-date studies were reviewed on three main aspects from safety perspectives in Southwest China large-scale hydropower projects.They are the dynamic evolution of hydraulic properties of the slope rock mass fractured in different scales, regularities of non-Darcy flows under extreme pressure differences in lofty rock slopes, and field hydraulic controls based on systematical and feedback modelling analyses.
In this paper, the experimental tests on seepage characteristics of both rock samples and rough-walled fractures as well as the corresponding simulation analyses were summarized in a relatively systematic way.Not only the damage change defined by crack growth but also the evolution of seepage characteristics were predicted using the models established by researchers.Furthermore, the advances in studying the non-linear regularities of high-speed seepage flows under high pressure gradient and the unsaturated water-holding characteristics within deformable rock mass are introduced for the sake of broader application to hydraulic motion problems in high-steep slopes.Based on these latest findings, the significant achievements are also concisely elaborated in promoting seepage control theory such as a complete system covering medium property, boundary condition, initial state, and multi-field coupling as well as seepage control practice through inverse modeling analysis during construction of large-scale hydropower stations.
The hydraulic characteristics of fractured hard rocks of a reservoir bank slope undergo progressive evolution in the course of dam abutment excavation, reservoir impounding and fluctuation operation, and discharge atomization in hydropower engineering.Accordingly, seepage control measures by hydro-structures such as drainage or water insulation curtains should be guided by an innovative scientific foundation which concerns dynamic transformation of parameters, laws, and models during life-cycle performance.
Although applying hydraulic structures to seepage flowing control in a large-scale fractured rock mass resulted in a distinct achievement in Chinese hydropower engineering, a profound progress on fundamental fractured rock hydraulics is still on the way.It is expected that competent geophysical means and integrated data from effective monitoring systems will play a more valuable role in discovering and understanding the science intertwining complex flows and fractured rocks in engineering fields.

Figure 1 .
Figure 1.A bird's eye view of Jinping I Hydropower Station by man-made satellite.

Figure 1 .
Figure 1.A bird's eye view of Jinping I Hydropower Station by man-made satellite.

Processes 2018, 6 ,Figure 2 .
Figure 2. Coupled triaxial compression and two-phase flow testing system for low permeable rocks [34]: (a) Experimental facility by Research & Development; (b) Tested curves as an example.

Figure 2 .
Figure 2. Coupled triaxial compression and two-phase flow testing system for low permeable rocks [34]: (a) Experimental facility by Research & Development; (b) Tested curves as an example.

Figure 3 .
Figure 3. Contour of simulated progressive permeability distribution along with stress-strain process [37]: (a) under maintained pore pressure 7 MPa; (b) eventual peak stress states under various maintained pore pressures.

Figure 3 .
Figure 3. Contour of simulated progressive permeability distribution along with stress-strain process [37]: (a) under maintained pore pressure 7 MPa; (b) eventual peak stress states under various maintained pore pressures.

Figure 4 .
Figure 4. Illustration of high-performance fine modelling on a rock sample stress-strain process using statistical meso-damage Finite Element Method [41]: (a) Simulated permeability evolution along with stress-strain process; (b) Progress cracking states on the middle cross-section at five typical moments; (c) Cracking states on three typical cross-sections at the peak-stress moment.

Figure 4 .
Figure 4. Illustration of high-performance fine modelling on a rock sample stress-strain process using statistical meso-damage Finite Element Method [41]: (a) Simulated permeability evolution along with stress-strain process; (b) Progress cracking states on the middle cross-section at five typical moments; (c) Cracking states on three typical cross-sections at the peak-stress moment.
, three mechanisms causing non-Darcy flow were revealed.They are namely inertia effect of water flow (a), pore swelling caused by dislocating rough fissure walls (b), and interaction between fluid-solid deformations (c).The physical parameter model for Forchheimer's law, critical Reynolds number model, and flow regime criterion were proposed and validated.Considering fracture morphology and fluid regime, the friction factor was given in a theoretical expression.Three flow zones covering fracture flow regimes were established.They are Darcy flow zone, Forchheimer flow zone, and complete turbulent zone [51,52].It should be noted that the existing practice in fields had followed the linear seepage rule of Darcy flow over a long history.The newly developed permeability estimate method for high-pressure packer testing based on the above findings achieved favorable effects in the 450 m deep practice at the central region of Hainan Island.

Figure 6 .
Figure 6.Experimental outcomes of solute migration and viscous fingering along a transparent rough fracture [52-56]: (a) Solute migration test and breakthrough curve; (b) developing process of viscous fingering.

Figure 7 .
Figure 7. Experimental regularity and theoretical study of non-Darcy flows along a deformable rough rock fracture: (a) non-linearity from inertia effect of water flow (b) non-linearity from pore swelling due to dislocating rough fissure walls (c) non-linearity from coupling fluid-solid deformations [50]; (d) zone division for flow regimes [57].

Figure 8 .
Figure 8. Types and features of P-Q curves in drilled borehole high-pressure Parker tests [59]: (a) Typical P-Q curve; (b) curve Type A; (c) curve Type B.

Figure 8 .
Figure 8. Types and features of P-Q curves in drilled borehole high-pressure Parker tests [59]: (a) Typical P-Q curve; (b) curve Type A; (c) curve Type B.

Figure 9 .
Figure 9. Large-scale physical experiment on deformation features of the reservoir bank under intensive atomization.

Figure 9 .
Figure 9. Large-scale physical experiment on deformation features of the reservoir bank under intensive atomization.

Figure 9 .Figure 10 .
Figure 9. Large-scale physical experiment on deformation features of the reservoir bank under intensive atomization.

Figure 12 .
Figure 12.Optimization of the design method for seepage flow control of high-steep slope engineering [67].

Figure 12 .
Figure 12.Optimization of the design method for seepage flow control of high-steep slope engineering [67].

Figure 14 .
Figure 14.Inverse modelling and enforcement to water gushing from a drainage gallery at the left bank with altitude 1995 m in project Jinping [68]: (a) the inverse analysis numerical model; (b) the measured flow rates out of the drainage holes in the drainage tunnel at 1595 m above sea level on July 23, 2013, when the reservoir water level was about 1800 m above sea level; (c) comparison of the measured and calculated hydraulic heads at piezometers installed on the downstream side of grouting curtain; (d) comparison of the measured and calculated discharges through the measuring weirs.

Figure 14 .
Figure 14.Inverse modelling and enforcement to water gushing from a drainage gallery at the left bank with altitude 1995 m in project Jinping [68]: (a) the inverse analysis numerical model; (b) the measured flow rates out of the drainage holes in the drainage tunnel at 1595 m above sea level on July 23, 2013, when the reservoir water level was about 1800 m above sea level; (c) comparison of the measured and calculated hydraulic heads at piezometers installed on the downstream side of grouting curtain; (d) comparison of the measured and calculated discharges through the measuring weirs.

Table 1 .
Height of the dam abutment slopes for some hydro-power stations in southwest China 1 .