# Seepage Evolution Model of the Fractured Rock Mass under High Seepage Pressure in Dam Foundation

Authors: Erfeng Zhao; Yufeng Jiang
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The seepage of the fractured rock mass in dam foundations involves complex fluid-structure coupling behavior, due to practical hydrogeological conditions. In this work, the seepage characteristics of the fractured rock mass and their correlations with the structural permeable mediums are experimentally explored to reveal the cracking effect on the hydromechanical properties firstly. Subsequently, the tangential and the compression creep damage constitutive models are, respectively, established by introducing a nonlinear viscoplastic body to improve the Nishihara model. Afterwards, an innovative evolution equation of the permeability coefficient considering the creep damage is proposed. It can indicate the time effect of the porosity, the permeability, and damage variables of the fractured rock mass under the long-term infiltration action of the hydraulic pressures. Ultimately, the proposed methods are applied to the seepage simulation on the dam foundation of the Longyangxia hydropower station and the significantly increased leakage is in good agreement with the measured values during the storage period. It was further confirmed that the crack expansion and penetration in the rock masses can be constantly intensified by the seepage pressures. The research results can provide a reference for engineering repair and supervision through controlling the permeability performance for long-term operations.


---

## Body

## 1. Introduction

Dam foundations are often located in high geostress and high hydraulic pressure environment with extremely complex hydrogeological conditions [1–3]. A large part of the deformation of the rock mass is primarily caused by various weak intercalations and geological faults. The mechanical evolution of the faults can obviously affect the leakage quantity, and the variation of the seepage field will change the internal fluid pressures [4–6]. Under the long-term seepage actions, the cracks in the fractured rock mass expand continuously and the seepage resistance decreases correspondingly to formulate low-resistance seepage channels. They can cause the stress redistribution, weaken the mechanical strength, and damage the stability. For example, the primary cause of the Malpasset dam failure was the stress concentration in local zones in the foundation during water storage [7]. Moreover, in the Meishan reservoir, the lateral seepage pressure on the weak intercalations increased gradually during high-water level operation and caused the arch stacks to move sideways abnormally in 1962 [8].In practical engineering projects, the seepage in the fractured rock mass involves complex fluid-structure coupling behavior [9, 10] and can usually be generalized to the seepage in fractured networks [11–13]. The coupling behavior between the seepage field and the stress field was simulated by Zhang et al. [14] using extended finite element method. Shi and Yang [15] established a coupled flow model to simulate the seepage property and the structural evolution of the fractured rock mass. The hydromechanical effect was studied by Figueiredo et al. [16] to characterize the permeability change at 1000 m deep rock mass. Gan and Elsworth [17] presented a continuum coupled model on the stress and the fluid flow in the fractured rock mass. Liu et al. [18] developed a numerical approach to estimate the shear effects on the fluid flow in the fractured networks. Pappalardo [19] acknowledged that the cooling behavior could affect the hydraulic conductivity by the infrared thermography. Furthermore, various numerical modelling methods of the seepage flow have been proposed and implemented to reveal the complex mechanism. Ren et al. [20] proposed an equivalent permeability method for modelling the fluid flow numerically. Sundararajan and Sankaran [21] developed a groundwater flow and transport model to estimate the seepage loss of Musi basin through simulation. Zhang et al. [22] produced physical models to simulate the seepage flow behavior in rough fractures with different fractal dimensions. The aforementioned studies have promoted further development of the seepage models on the fractured rock mass.Owing to various disturbances in dam foundations, the propagation of microcracks and joint fissures in the fractured rock mass produce strong heterogeneity and nonlinearity of flow characteristics [23–26]. Relevant experimental studies are generally subject to technical means of obtaining effective and implicit test information [27, 28], including structural properties, seepage deformations, and seepage patterns. Considering the correlation verification of multidimension test methods, it is urgent to study the seepage evolution characteristics of the fractured water and the structural control effect. In addition, the seepage of the fractured rock mass owns obvious characteristics of dynamic evolution due to the transformation of the seepage pressures [29, 30]. Few achievements have been made on the seepage damage effect of the foundation rock mass in controlling its permeability performance.In this study, the seepage characteristics of the fractured rock mass and their correlations with the structural permeable mediums are firstly extracted by combining with seepage tests inSection 2. Secondly, a nonlinear viscoplastic body is newly introduced to improve the Nishihara model to establish novel creep damage constitutive models in Section 3. Thirdly, an innovative seepage evolution model considering the creep damage is proposed in Section 4. The proposed methods are finally validated and applied to simulate the seepage evolution of the practical dam foundation in Section 5.


## 2. Seepage Characteristics of the Fractured Rock Mass

According to triaxial compression tests of the fractured rock, the progressive failure is often accompanied by the compression, crack, expansion, and mutual intersection of internal fractures [31, 32]. The permeability changes obviously with the formation of the fractured networks, owing to the expansion and the penetration of primary cracks. They can formulate seepage channels in the conglomerate and decrease the seepage resistance, resulting in the enhancement of the permeability.


### 2.1. Permeability of the Soft and Hard Rocks

Under certain constraint pressures and pore pressures, the permeability change can be obtained by rock electrohydraulic servocontrolled permeability tests (Figure1(a)) through applying axial pressures gradually [33]. The permeability evolutions of the hard and soft rocks in the whole stress-strain process are shown in Figure 1 [34, 35]. Figure 1 reflects the close relationship between the permeability and the displacement, indicating the permeability change in different stages. In the elastic stage (primary microcrack closure stage), the permeability coefficient decreases slightly with the increase of the axial stress. After entering the nonlinear strain hardening stage, the permeability coefficient gradually increases and its maximum lags behind the peak strength, due to the propagation of microcracks and the new formation of macrocracks. Subsequently, with the increase of the axial strain, the permeability coefficient may continue to increase and decrease gradually or sharply, depending on rock types, fabric, and void development. The soft rock shows strong plastic rheology after failure and the permeability attenuation range is generally large, while the permeability of the hard rock is relatively stable after failure, with a small change range. They are both caused by the internal correlations between the rock permeability and the crack propagation.Figure 1
Permeability variations of the hard and soft rocks. (a) Electrohydraulic servocontrolled testing system; (b) hard rock: fine sandstone; (c) soft rock: bauxite.
(a)(b)(c)


### 2.2. Permeability of the Fractured Rock Mass

The permeability measurement of the fractured rock mass can be obtained by laboratory permeability tests or in situ water pumping tests. The seepage pressure and the fractured permeability channel shape are two critical factors affecting the seepage evolution of the fractured rock mass.(1)
According to the laboratory test results (Figure2) [34], the seepage-guiding state of the penetrated fracture channel is generally characterized by the high seepage resistance, and the seepage of the fractured network is turbulent. Under certain seepage pressures, the seepage velocity decreases with the increase of the seepage channel width, while the seepage pressure along the fractured network gradually decreases with the increase of the seepage distance. They indicate that the decrease of the seepage pressure and the increase of the seepage resistance are correlated and can be transformed.(2)
In in situ water pumping tests [34, 35], the selected sections include the normal fault zone, the reverse fault zone, and the unloading disturbance zone. The dual-porosity seepage method is taken to test the correlations of the exerted water pressures and the flow duration change under high seepage pressures. The whole process of the rock mass reaching the conduction seepage state is shown in Figure 3. The permeability primarily depends on structural planes, that is, the permeability and mechanical characteristics of the fractures change synchronously. Meanwhile, the mechanical, physical, and chemical effects of the internal water in the fractures, microcracks, and rock pores can decrease the structural strength. Moreover, the additional seepage volume force will also change the rock mass displacement and adjust the internal stress states.Figure 2
Influence of the seepage channel width and the particle size on the seepage velocity. (a) The seepage testing model using the acrylic glass and the battery carbon rods; (b) crack width: 2 mm; (c) crack width: 5 mm.
(a)(b)(c)Figure 3
In situ water pumping tests on the fractured rock mass. (a) Normal fault; (b) reverse fault; (c) unloading disturbance zone.
(a)(b)(c)The in situ test results can indicate the permeability of the rock masses in practical geological conditions. As shown in Figure3, the seepage evolution and the structural control effect can be regarded as the essential basis for evaluating the permeability and antiseepage failure strength of the fractured rock masses.


## 2.1. Permeability of the Soft and Hard Rocks

Under certain constraint pressures and pore pressures, the permeability change can be obtained by rock electrohydraulic servocontrolled permeability tests (Figure1(a)) through applying axial pressures gradually [33]. The permeability evolutions of the hard and soft rocks in the whole stress-strain process are shown in Figure 1 [34, 35]. Figure 1 reflects the close relationship between the permeability and the displacement, indicating the permeability change in different stages. In the elastic stage (primary microcrack closure stage), the permeability coefficient decreases slightly with the increase of the axial stress. After entering the nonlinear strain hardening stage, the permeability coefficient gradually increases and its maximum lags behind the peak strength, due to the propagation of microcracks and the new formation of macrocracks. Subsequently, with the increase of the axial strain, the permeability coefficient may continue to increase and decrease gradually or sharply, depending on rock types, fabric, and void development. The soft rock shows strong plastic rheology after failure and the permeability attenuation range is generally large, while the permeability of the hard rock is relatively stable after failure, with a small change range. They are both caused by the internal correlations between the rock permeability and the crack propagation.Figure 1
Permeability variations of the hard and soft rocks. (a) Electrohydraulic servocontrolled testing system; (b) hard rock: fine sandstone; (c) soft rock: bauxite.
(a)(b)(c)


## 2.2. Permeability of the Fractured Rock Mass

The permeability measurement of the fractured rock mass can be obtained by laboratory permeability tests or in situ water pumping tests. The seepage pressure and the fractured permeability channel shape are two critical factors affecting the seepage evolution of the fractured rock mass.(1)
According to the laboratory test results (Figure2) [34], the seepage-guiding state of the penetrated fracture channel is generally characterized by the high seepage resistance, and the seepage of the fractured network is turbulent. Under certain seepage pressures, the seepage velocity decreases with the increase of the seepage channel width, while the seepage pressure along the fractured network gradually decreases with the increase of the seepage distance. They indicate that the decrease of the seepage pressure and the increase of the seepage resistance are correlated and can be transformed.(2)
In in situ water pumping tests [34, 35], the selected sections include the normal fault zone, the reverse fault zone, and the unloading disturbance zone. The dual-porosity seepage method is taken to test the correlations of the exerted water pressures and the flow duration change under high seepage pressures. The whole process of the rock mass reaching the conduction seepage state is shown in Figure 3. The permeability primarily depends on structural planes, that is, the permeability and mechanical characteristics of the fractures change synchronously. Meanwhile, the mechanical, physical, and chemical effects of the internal water in the fractures, microcracks, and rock pores can decrease the structural strength. Moreover, the additional seepage volume force will also change the rock mass displacement and adjust the internal stress states.Figure 2
Influence of the seepage channel width and the particle size on the seepage velocity. (a) The seepage testing model using the acrylic glass and the battery carbon rods; (b) crack width: 2 mm; (c) crack width: 5 mm.
(a)(b)(c)Figure 3
In situ water pumping tests on the fractured rock mass. (a) Normal fault; (b) reverse fault; (c) unloading disturbance zone.
(a)(b)(c)The in situ test results can indicate the permeability of the rock masses in practical geological conditions. As shown in Figure3, the seepage evolution and the structural control effect can be regarded as the essential basis for evaluating the permeability and antiseepage failure strength of the fractured rock masses.


## 3. Creep Damage Model of the Fractured Rock Mass

In the dam foundations, the joints, fissures, and gaps are not only the migration channels of the groundwater but also the primary mediums of the foundation deformation. The displacement increases with time under continuous loads [36]. Hence, the interior fluid flow is mainly controlled by the structural properties and the connectivity of the fractured networks. The test results of the granite foundation of the Huangdeng Dam are shown in Figure 4 [37], including the initial strain, as well as the steady and accelerated creep strain under different confining pressures. In Figure 4(a), the irreversible plastic deformation became clearer with the increase of the stress levels and the accumulation creep effect, indicating that the internal structure adjustment is more severe under the high axial stress levels. When the stress levels reached 100 MPa and 120 MPa, respectively, the creep failure occurred due to the damage accumulation. In Figure 4(b), under the confining stresses of 2 MPa and 6 MPa, the seepage flow rates decreased remarkably, due to the variations of the fractures and pores under the initial loading, and they became stable gradually after 0.21 h and 0.44 h, respectively.Figure 4
Test results of the granite foundation under different confining pressures. (a) Creep properties for volcanic breccia; (b) seepage flow rate at the instantaneous strain stage.
(a)(b)Under the long-term infiltration action of high hydraulic pressures, the foundation creep leads to certain changes in microcracks of the rock mass and geometric shapes of the structural planes in the faults. They can affect the hydraulic conductivity and change the seepage laws. The creep effect can be simulated numerically combining with the bionic algorithms [38]. On the contrary, the change of the seepage field can cause the variation of the stress field and the creep of the dam foundation accordingly. As shown in Figure 4, the creep displacement finally tends to be stable when the shear load is small. It can be assumed that the fractured rock mass will not be damaged in the elastic stage and the damage only occurs after entering the viscoplastic stage. Karalar and Çavuşli [39, 40] simulated the creep settlements of the Ilısu concrete faced rockfill dam and the nonlinear seismic behaviors of the Boyabat concrete gravity dam with special viscoplastic models using the FLAC3D software based on the finite difference method. Therefore, the creep damage model should contain the stress-strain constitutive equation and the damage evolution equation [41, 42]. Because the creep damage of the fractured rock mass can lead to stiffness deteriorating and strength weakening, the establishment of the damage function and the yield function should consider the damage evolution.


### 3.1. Creep Damage Function

The increase of the viscoplastic creep displacement is accompanied by the damage accumulation in Figure4(a). The evolution equation of damage variables should be firstly established to indicate the creep damage. According to the second law of thermodynamics, it is a function of the damage dissipation potential, defined as follows:(1)D=I−exp−BY−Y01/n,where I is the unit matrix, B and n are material parameters, Y0 is the initial dissipated energy, and Y is the total dissipated energy after the creep damage.Due to the functional relationship between the damage dissipation energy and the tangential plastic displacement, the elastic-plastic tangential damageDτ and the normal damage Dv are(2)Dτ=1−exp−κζDR,Dv=1−exp−κsζDRs,where κ, R, κs, and Rs are material parameters and ζD is the plastic shear displacement.Assuming that the shear stress in the complete damage isτc and the normal displacement is vc, the following equation can be obtained:(3)τa=Dττc+1−Dττi,va=ve+Dvvc,where τa and va are the actual shear stress and the normal displacement, respectively; τi is the shear stress without considering the damage.Equation (3) is consistent with the elastic-plastic damage function. Because there is a functional relationship between the damage dissipation energy and its plastic work accumulation, the functional relationship between the creep damage dissipation energy and the viscoplastic displacement can be constructed, regardless of the normal displacement. Namely,(4)Y=a′ωvpb′,where ωvp=∫du′vpdu′vp+dvvpdv′vp1/2; du′vp and dv′vp are the viscoplastic displacement increments in the local axes of x′ and y′, respectively; a′ and b′ are the material parameters.Assuming that the tangential damage is isotropic, without original damage, the evolution function of the tangential creep damage and the normal creep damage in the foundation faults can be defined as follows:(5)Dτ=1−exp−aωvpb,Dv=1−exp−cωvpd,where a, b, c, and d are material parameters.According to equation (5), the tangential and normal damage degrees are both equal to 0 when the viscoplastic displacement is 0, while they are nearly 1.0 when the fractured rock mass is completely damaged. The friction coefficient of the fractured rock mass changes slightly after the shear damage and is 0.95–1-fold of the original. When the cohesion almost disappears completely, only 0.1-0.2-fold of the original is preserved. The shear creep damage yield function is(6)f=τ+μσn−αc,where μ is the friction coefficient; σn is the normal stress (the positive value represents the tensile stress and the compressive stress otherwise); c is the cohesion; α is the damage function: α=1−Dτ. For the undamaged rock mass, Dτ=0 and α=1; consequently, f=τ+μσn−c, the traditional Mohr–Coulomb yield criterion; with the increase of the creep displacement, Dτ increases, whereas α decreases until the rock mass is completely damaged; that is, α=0 and f=τ+μσn.


### 3.2. Creep Damage Constitutive Model

For the fractured rock mass under external loads, the tangential displacement is different from the normal one. The tangential model and the normal model of the creep damage should be established, respectively. The viscoelastic-plastic constitutive equation without considering the damage is(7)σ˙ij=Dijkl0ε˙kl−ε˙klvp,where Dijkl0 is the stiffness tensor, σ˙ij is the stress tensor rate, ε˙kl is the total strain rate, and ε˙klvp is the viscous strain rate.According to the equivalent strain hypothesis, the viscoelastic-plastic damage constitutive equation can be deduced from equation (7), replacing the apparent stress by the effective stress.(8)σ˙ij∗=Dijkl0ε˙kl−ε˙klvp.Assuming that the tangential elastic stiffness coefficient iskx′i=ky′i=ksi under the positive stress σn (x′ and y′ are mutually perpendicular axes parallel to the structural planes) and that ksc remains in the complete damage condition, the tangential elastic stiffness coefficient in the evolution is defined as follows:(9)ks=Dτksc+1−Dτksi.According to equation (8), the elastic strain rate can be deduced.(10)ε˙ije=ε˙ij−ε˙ijv=Cijkl0σ˙kl∗,where Cijkl0 is the flexibility tensor.WhenCijkl0 is replaced by the flexibility tensor Cijkl0−d of the damaged medium and the effective stress is replaced by the apparent stress σkl, the following equation can be obtained:(11)ε˙ije=Cijkl0−dσ˙kl+C˙ijkl0−dσkl.The elastic strain rate for the shear creep damage is(12)γ˙je=1ksτ˙j−1ks2k˙sτj=1ksτ˙j+1ks2ksi−kscD˙ττj,where j can be x′ or y′.When the shear stress is small in the fractured rock mass, the shear creep deformation tends to be stable, but the shear deformation after unloading is plastic and does not decrease with time, due to the relative slip and dislocation during the creep development. The tangential creep damage equation can be deduced by the model shown in Figure5, through improving the Nishihara model with a nonlinear viscoplastic body. The elastic body simulates the instantaneous elastic deformation, the Kelvin body simulates the decaying creep deformation, the plastic body simulates the plastic yield, and the viscous body simulates the irreversible creep. The viscoelastic strain rate can be got according to the equivalent strain assumption.(13)γ˙jve=τ˜jηkexp−kkηkt,where τ˜j is the effective shear stress: τ˜j=ks/ksiτj=1−Dτ+ks/ksiDττj.Figure 5
Shear creep model of the fractured rock mass.The viscoelastic strain can be got as follows:(14)γj,t+Δtve=γj,tveexp−kkηkΔt+τ˜jkk1−exp−kkηkΔt.When the yielding phenomenon occurs, the plastic potential functionQ should be introduced:(15)Q=Qσ,ε˙vp.The viscoplastic flow rule is(16)γ˙jvp=γ<f>∂Q∂τj,where γ is the flow coefficient; <f>=f,f>0,0,f≤0.For the creep model shown in Figure5, γ=1/ηs, where ηs is the viscosity coefficient after the creep yield. The associated flow rule is adopted.(17)γ˙jvp=γ<f>∂Q∂τj=1ηs<f>∂Q∂τj=1ηs<f>τjτx′2+τy′2.Whenτx′ and τy′ are constants and f=τx′2+τy′21/2+μσz′−1−Dτc>0, σz′ is constant and z′ is the axis perpendicular to the structural plane. With the accumulation of the creep damage, f and γ˙vp will increase, indicating the creep acceleration.On this basis, the tangential creep damage constitutive equation can be deduced.(18)γ˙j=γ˙je+γ˙jv+γ˙jve+γ˙jvp=1ksτ˙j+1ks2ksi−kscD˙ττj+1−Dτ+ks/ksiDτηMτj+1−Dτ+ks/ksiDττjηk⋅exp−kkηkt+1ηs<f>τjτx′2+τy′2.Influenced by the pores in the fractured rock mass, the significant compaction can be produced under the normal force, that is, the instantaneous elastic deformation occurs first, and then the creep rate decreases, and finally the creep deformation tends to be stable. Most of the normal compression deformation cannot be recovered after unloading. The attenuation creep characteristics along the normal direction can be simulated by the model shown in Figure6. The elastic body reflects the instantaneous deformation, the Kelvin body reflects the decay creep and the hysteresis recovery deformation after unloading, and the unrecoverable creep can be indicated by the viscous body. The normal displacement can be defined as(19)εz′=εz′e+εz′ve+Dvεz′c,where the viscous displacement is(20)εz′ve=σz′k2,n1−e−k2,n/ηn+∫0tσz′ηMtdt.Figure 6
Normal compression creep model of the fractured rock mass.In Figure6, in addition to the instantaneous recovery during unloading, the hysteresis recovery time is actually determined by the two viscous bodies. The normal compression creep damage constitutive equation is(21)εz′=σz′kn+σz′k2,n1−e−k2,n/ηnt+∫t=0tσz′ηMtdt+Dvεz′0exp−kσz′.


## 3.1. Creep Damage Function

The increase of the viscoplastic creep displacement is accompanied by the damage accumulation in Figure4(a). The evolution equation of damage variables should be firstly established to indicate the creep damage. According to the second law of thermodynamics, it is a function of the damage dissipation potential, defined as follows:(1)D=I−exp−BY−Y01/n,where I is the unit matrix, B and n are material parameters, Y0 is the initial dissipated energy, and Y is the total dissipated energy after the creep damage.Due to the functional relationship between the damage dissipation energy and the tangential plastic displacement, the elastic-plastic tangential damageDτ and the normal damage Dv are(2)Dτ=1−exp−κζDR,Dv=1−exp−κsζDRs,where κ, R, κs, and Rs are material parameters and ζD is the plastic shear displacement.Assuming that the shear stress in the complete damage isτc and the normal displacement is vc, the following equation can be obtained:(3)τa=Dττc+1−Dττi,va=ve+Dvvc,where τa and va are the actual shear stress and the normal displacement, respectively; τi is the shear stress without considering the damage.Equation (3) is consistent with the elastic-plastic damage function. Because there is a functional relationship between the damage dissipation energy and its plastic work accumulation, the functional relationship between the creep damage dissipation energy and the viscoplastic displacement can be constructed, regardless of the normal displacement. Namely,(4)Y=a′ωvpb′,where ωvp=∫du′vpdu′vp+dvvpdv′vp1/2; du′vp and dv′vp are the viscoplastic displacement increments in the local axes of x′ and y′, respectively; a′ and b′ are the material parameters.Assuming that the tangential damage is isotropic, without original damage, the evolution function of the tangential creep damage and the normal creep damage in the foundation faults can be defined as follows:(5)Dτ=1−exp−aωvpb,Dv=1−exp−cωvpd,where a, b, c, and d are material parameters.According to equation (5), the tangential and normal damage degrees are both equal to 0 when the viscoplastic displacement is 0, while they are nearly 1.0 when the fractured rock mass is completely damaged. The friction coefficient of the fractured rock mass changes slightly after the shear damage and is 0.95–1-fold of the original. When the cohesion almost disappears completely, only 0.1-0.2-fold of the original is preserved. The shear creep damage yield function is(6)f=τ+μσn−αc,where μ is the friction coefficient; σn is the normal stress (the positive value represents the tensile stress and the compressive stress otherwise); c is the cohesion; α is the damage function: α=1−Dτ. For the undamaged rock mass, Dτ=0 and α=1; consequently, f=τ+μσn−c, the traditional Mohr–Coulomb yield criterion; with the increase of the creep displacement, Dτ increases, whereas α decreases until the rock mass is completely damaged; that is, α=0 and f=τ+μσn.


## 3.2. Creep Damage Constitutive Model

For the fractured rock mass under external loads, the tangential displacement is different from the normal one. The tangential model and the normal model of the creep damage should be established, respectively. The viscoelastic-plastic constitutive equation without considering the damage is(7)σ˙ij=Dijkl0ε˙kl−ε˙klvp,where Dijkl0 is the stiffness tensor, σ˙ij is the stress tensor rate, ε˙kl is the total strain rate, and ε˙klvp is the viscous strain rate.According to the equivalent strain hypothesis, the viscoelastic-plastic damage constitutive equation can be deduced from equation (7), replacing the apparent stress by the effective stress.(8)σ˙ij∗=Dijkl0ε˙kl−ε˙klvp.Assuming that the tangential elastic stiffness coefficient iskx′i=ky′i=ksi under the positive stress σn (x′ and y′ are mutually perpendicular axes parallel to the structural planes) and that ksc remains in the complete damage condition, the tangential elastic stiffness coefficient in the evolution is defined as follows:(9)ks=Dτksc+1−Dτksi.According to equation (8), the elastic strain rate can be deduced.(10)ε˙ije=ε˙ij−ε˙ijv=Cijkl0σ˙kl∗,where Cijkl0 is the flexibility tensor.WhenCijkl0 is replaced by the flexibility tensor Cijkl0−d of the damaged medium and the effective stress is replaced by the apparent stress σkl, the following equation can be obtained:(11)ε˙ije=Cijkl0−dσ˙kl+C˙ijkl0−dσkl.The elastic strain rate for the shear creep damage is(12)γ˙je=1ksτ˙j−1ks2k˙sτj=1ksτ˙j+1ks2ksi−kscD˙ττj,where j can be x′ or y′.When the shear stress is small in the fractured rock mass, the shear creep deformation tends to be stable, but the shear deformation after unloading is plastic and does not decrease with time, due to the relative slip and dislocation during the creep development. The tangential creep damage equation can be deduced by the model shown in Figure5, through improving the Nishihara model with a nonlinear viscoplastic body. The elastic body simulates the instantaneous elastic deformation, the Kelvin body simulates the decaying creep deformation, the plastic body simulates the plastic yield, and the viscous body simulates the irreversible creep. The viscoelastic strain rate can be got according to the equivalent strain assumption.(13)γ˙jve=τ˜jηkexp−kkηkt,where τ˜j is the effective shear stress: τ˜j=ks/ksiτj=1−Dτ+ks/ksiDττj.Figure 5
Shear creep model of the fractured rock mass.The viscoelastic strain can be got as follows:(14)γj,t+Δtve=γj,tveexp−kkηkΔt+τ˜jkk1−exp−kkηkΔt.When the yielding phenomenon occurs, the plastic potential functionQ should be introduced:(15)Q=Qσ,ε˙vp.The viscoplastic flow rule is(16)γ˙jvp=γ<f>∂Q∂τj,where γ is the flow coefficient; <f>=f,f>0,0,f≤0.For the creep model shown in Figure5, γ=1/ηs, where ηs is the viscosity coefficient after the creep yield. The associated flow rule is adopted.(17)γ˙jvp=γ<f>∂Q∂τj=1ηs<f>∂Q∂τj=1ηs<f>τjτx′2+τy′2.Whenτx′ and τy′ are constants and f=τx′2+τy′21/2+μσz′−1−Dτc>0, σz′ is constant and z′ is the axis perpendicular to the structural plane. With the accumulation of the creep damage, f and γ˙vp will increase, indicating the creep acceleration.On this basis, the tangential creep damage constitutive equation can be deduced.(18)γ˙j=γ˙je+γ˙jv+γ˙jve+γ˙jvp=1ksτ˙j+1ks2ksi−kscD˙ττj+1−Dτ+ks/ksiDτηMτj+1−Dτ+ks/ksiDττjηk⋅exp−kkηkt+1ηs<f>τjτx′2+τy′2.Influenced by the pores in the fractured rock mass, the significant compaction can be produced under the normal force, that is, the instantaneous elastic deformation occurs first, and then the creep rate decreases, and finally the creep deformation tends to be stable. Most of the normal compression deformation cannot be recovered after unloading. The attenuation creep characteristics along the normal direction can be simulated by the model shown in Figure6. The elastic body reflects the instantaneous deformation, the Kelvin body reflects the decay creep and the hysteresis recovery deformation after unloading, and the unrecoverable creep can be indicated by the viscous body. The normal displacement can be defined as(19)εz′=εz′e+εz′ve+Dvεz′c,where the viscous displacement is(20)εz′ve=σz′k2,n1−e−k2,n/ηn+∫0tσz′ηMtdt.Figure 6
Normal compression creep model of the fractured rock mass.In Figure6, in addition to the instantaneous recovery during unloading, the hysteresis recovery time is actually determined by the two viscous bodies. The normal compression creep damage constitutive equation is(21)εz′=σz′kn+σz′k2,n1−e−k2,n/ηnt+∫t=0tσz′ηMtdt+Dvεz′0exp−kσz′.


## 4. Seepage Evolution Model considering the Creep Damage

Generally speaking, the seepage in the fractured rock mass conforms to Darcy’s law, and the permeability coefficient is a critical index. The relationship between the porosity, the permeability, and damage variables is shown in Figure7 [43], where β is related to the moisture content and the temperature. The permeability gradually increases with the damage accumulation and tends to be stable when the dominant penetrating pores are formed. Most previous porous media models assume that the porosity and the permeability coefficients are material constants independent of time. However, the gradual damage of the foundation rock mass (including the grouting curtain rock mass) will greatly affect the seepage evolution. Here, the influence of the creep damage on the permeability is used to characterize the increase of the hydraulic conductivity, which is primarily caused by the spatial change of the water flow inside the filling. The damage of the grouting curtain rock mass is used to characterize the degradation of its antiseepage function.Figure 7
Relationship between porosity, permeability, and damage variables in the fractured rock mass.


### 4.1. Evolution Equation of the Permeability Coefficient

As shown in Figure8, the lengths of three sides of a damaged unit are dx1, dx2, and dx3, and the anisotropic damage characteristics can be indicated by the damage principal value. Damage variables are defined as(22)Di=Ai−Ai∗Ai=dsi∗dxjdxk,where Di (i = 1, 2, 3) is the primary damage variable; dsi∗ is the reduction area of the section surface dxidxk, caused by the porosity variation.Figure 8
Diagram of a damaged unit.The relationship between the porosityn and the damage variable D is(23)n=D1D2D31/2.For isotropic materials,D1=D2=D3=D and n=D3/2, the relationship between the change rate of the porosity and the damage developing rate is(24)n˙=32D1/2D˙.According to equation (23), the porosity evolution can be described in view of the continuum damage mechanics due to the quantitative relationship between n and D. Assuming that the permeability coefficient depends on the infilling properties, regarded as the porous rock, the isotropous permeability coefficient can be defined as(25)k=gμd2180n31−n2,where μ is the viscosity coefficient, d is the average size of solid particles, and n is the porosity of the infillings.Due to the variation of the porosityn along the depth z, where n=nz, the permeability coefficient can be(26)kz=gμd2180nz31−nz2.Therefore, the relation equations of the permeability coefficient, the damage variable, and their change rates are(27)k=gμd2180D9/21−D3/22=gμd2180D9/21−2D3/2+D3,k˙=gμd21803n21−n2n˙+2n31−nn˙1−n4=gμd21809D7/21−2D3/2+D3D˙+6D51−D3/2D˙21−D3/24=gμd21809D7/2−2D5+D13/2D˙+6D5−D13/2D˙21−D3/24.In view of the permeability evolution of the fractured rock mass, with respect to the timet, the above variables can be transformed as(28)nz,t=Dz,t3/2,(29)n˙z,t=32Dz,t1/2D˙z,t,(30)kz,t=gμd2180Dz,t9/21−Dz,t3/22=gμd2180Dz,t9/21−2Dz,t3/2+Dz,t3,(31)k˙z,t=gμd21803nz,t21−nz,t2n˙z,t+2nz,t31−nz,tn˙z,t1−nz,t4=gμd21809Dz,t7/2−2Dz,t5+Dz,t13/2D˙z,t+6Dz,t5−Dz,t13/2D˙z,t21−Dz,t3/24.Equations (30) and (31) are evolution equations of the permeability coefficient represented by D and D˙. When D = 0 and k(z, t) = 0, the rock mass is impervious and the porosity is n(t, z) = 0. The initial damage is unavoidable in practical bedrocks with certain permeability. The grouting curtain will lose the impervious effectiveness when its permeability coefficient increases to that of the bedrocks, D < 1 for this condition.


### 4.2. Numerical Simulation of the Seepage Evolution

The anisotropy equivalent permeability coefficient tensor should be known in advance in the numerical simulation. For infinitely extended and regularly arranged structural planes of the foundation faults, the equivalent permeability tensor is(32)kij=∑l=1mgbl312μSlδij−nilnjl,where m is the grouping number of cracks; bl is the equivalent hydraulic gap width of the lth group of fractures; Sl is the spacing of cracks; δij is the Kronecker symbol; nil (i = 1, 2, 3) is the cosine along the normal direction; μ is the kinematic viscosity coefficient.If the viscous displacement is equivalent to an instantaneous stress increment, the equivalent stress increment of the viscous strainεklv will be(33)Δσijv=dijklεklv.LetKx′, Ky′, and Kz′ be principal permeability coefficients of the permeability coefficient tensor. The effective principal stresses σx′s, σy′s, and σz′s along the penetration principal axes can be obtained from the total stress σij+Δσijv. According to the relationship between the permeability coefficient tensor and the stress K=K0eλσ, the permeability coefficient under continuous loading is(34)K′=Kx′λσx′sKy′λσy′sKz′λσz′s,where λ is an influence coefficient, determined by test results.The rotation matrix between the infiltration principal axes and the global coordinates is recorded asT. The permeability coefficient tensor is(35)K=TTK′T.The seepage motion equation for the fractured rock mass is(36)Kij,th,jc,i=0,where Kij,t is the permeability tensor after the creep damage, deduced by equations (34) and (35); hc is the water head in continuous rock mass.For 8-node isoparametric element, the governing equation is(37)Kcehce=Qce,where the component of Kce is(38)Ka′b′c=∫ΩeNa,iKij,tNb,jdΩ,a′,b′=1∼8;i,j=1∼3.The seepage influence of the fractured rock mass on its stress is considered as(39)F⟶=F⟶s+F⟶d=γwgradz−γwgradh,where γw is the bulk density of water; F⟶s is the hydrostatic pressure, and its effect can be equivalent by an effective stress; F⟶d is the hydrodynamic pressure, considered as an external load in the calculation.The equilibrium equation of the microelement considering the seepage is(40)∂σx∂x+∂τxy∂y+∂τxz∂z−∂h∂x=0,∂τxy∂x+∂σy∂y+∂τyz∂z−∂h∂y=0,∂τxz∂x+∂τyz∂y+∂σz∂z−∂h∂z+γw−γr=0,where γr is the saturated bulk density of rock; h is the total water head.Combining the stress-strain constitutive equation and the geometric equation, the equilibrium equations expressed by the displacement component and the water head are(41)G∇2δx+G1−2υ∂εv∂x−∂h∂x+X0=0,G∇2δy+G1−2υ∂εv∂y−∂h∂y+Y0=0,G∇2δz+G1−2υ∂εv∂z−∂h∂z+γw−γr=0,where G=E/21+υ is the shear modulus; ∇2=∂2/∂x2+∂2/∂y2+∂2/∂z2 is the Laplace operator; εv=εx+εy+εz=∂δx/∂x+∂δy/∂y+∂δz/∂z is the volumetric strain; X0, Y0, and Z0 are equivalent volume forces caused by the initial strain ε0.The equilibrium equation of the whole unit matrix derived from equation (41) is(42)Keδe=F1e+F0e+K′ehe,where F1e is the joint force caused by external loads; F0e is the joint force caused by the initial strain; K′e is the coupling matrix of the seepage and stress.


## 4.1. Evolution Equation of the Permeability Coefficient

As shown in Figure8, the lengths of three sides of a damaged unit are dx1, dx2, and dx3, and the anisotropic damage characteristics can be indicated by the damage principal value. Damage variables are defined as(22)Di=Ai−Ai∗Ai=dsi∗dxjdxk,where Di (i = 1, 2, 3) is the primary damage variable; dsi∗ is the reduction area of the section surface dxidxk, caused by the porosity variation.Figure 8
Diagram of a damaged unit.The relationship between the porosityn and the damage variable D is(23)n=D1D2D31/2.For isotropic materials,D1=D2=D3=D and n=D3/2, the relationship between the change rate of the porosity and the damage developing rate is(24)n˙=32D1/2D˙.According to equation (23), the porosity evolution can be described in view of the continuum damage mechanics due to the quantitative relationship between n and D. Assuming that the permeability coefficient depends on the infilling properties, regarded as the porous rock, the isotropous permeability coefficient can be defined as(25)k=gμd2180n31−n2,where μ is the viscosity coefficient, d is the average size of solid particles, and n is the porosity of the infillings.Due to the variation of the porosityn along the depth z, where n=nz, the permeability coefficient can be(26)kz=gμd2180nz31−nz2.Therefore, the relation equations of the permeability coefficient, the damage variable, and their change rates are(27)k=gμd2180D9/21−D3/22=gμd2180D9/21−2D3/2+D3,k˙=gμd21803n21−n2n˙+2n31−nn˙1−n4=gμd21809D7/21−2D3/2+D3D˙+6D51−D3/2D˙21−D3/24=gμd21809D7/2−2D5+D13/2D˙+6D5−D13/2D˙21−D3/24.In view of the permeability evolution of the fractured rock mass, with respect to the timet, the above variables can be transformed as(28)nz,t=Dz,t3/2,(29)n˙z,t=32Dz,t1/2D˙z,t,(30)kz,t=gμd2180Dz,t9/21−Dz,t3/22=gμd2180Dz,t9/21−2Dz,t3/2+Dz,t3,(31)k˙z,t=gμd21803nz,t21−nz,t2n˙z,t+2nz,t31−nz,tn˙z,t1−nz,t4=gμd21809Dz,t7/2−2Dz,t5+Dz,t13/2D˙z,t+6Dz,t5−Dz,t13/2D˙z,t21−Dz,t3/24.Equations (30) and (31) are evolution equations of the permeability coefficient represented by D and D˙. When D = 0 and k(z, t) = 0, the rock mass is impervious and the porosity is n(t, z) = 0. The initial damage is unavoidable in practical bedrocks with certain permeability. The grouting curtain will lose the impervious effectiveness when its permeability coefficient increases to that of the bedrocks, D < 1 for this condition.


## 4.2. Numerical Simulation of the Seepage Evolution

The anisotropy equivalent permeability coefficient tensor should be known in advance in the numerical simulation. For infinitely extended and regularly arranged structural planes of the foundation faults, the equivalent permeability tensor is(32)kij=∑l=1mgbl312μSlδij−nilnjl,where m is the grouping number of cracks; bl is the equivalent hydraulic gap width of the lth group of fractures; Sl is the spacing of cracks; δij is the Kronecker symbol; nil (i = 1, 2, 3) is the cosine along the normal direction; μ is the kinematic viscosity coefficient.If the viscous displacement is equivalent to an instantaneous stress increment, the equivalent stress increment of the viscous strainεklv will be(33)Δσijv=dijklεklv.LetKx′, Ky′, and Kz′ be principal permeability coefficients of the permeability coefficient tensor. The effective principal stresses σx′s, σy′s, and σz′s along the penetration principal axes can be obtained from the total stress σij+Δσijv. According to the relationship between the permeability coefficient tensor and the stress K=K0eλσ, the permeability coefficient under continuous loading is(34)K′=Kx′λσx′sKy′λσy′sKz′λσz′s,where λ is an influence coefficient, determined by test results.The rotation matrix between the infiltration principal axes and the global coordinates is recorded asT. The permeability coefficient tensor is(35)K=TTK′T.The seepage motion equation for the fractured rock mass is(36)Kij,th,jc,i=0,where Kij,t is the permeability tensor after the creep damage, deduced by equations (34) and (35); hc is the water head in continuous rock mass.For 8-node isoparametric element, the governing equation is(37)Kcehce=Qce,where the component of Kce is(38)Ka′b′c=∫ΩeNa,iKij,tNb,jdΩ,a′,b′=1∼8;i,j=1∼3.The seepage influence of the fractured rock mass on its stress is considered as(39)F⟶=F⟶s+F⟶d=γwgradz−γwgradh,where γw is the bulk density of water; F⟶s is the hydrostatic pressure, and its effect can be equivalent by an effective stress; F⟶d is the hydrodynamic pressure, considered as an external load in the calculation.The equilibrium equation of the microelement considering the seepage is(40)∂σx∂x+∂τxy∂y+∂τxz∂z−∂h∂x=0,∂τxy∂x+∂σy∂y+∂τyz∂z−∂h∂y=0,∂τxz∂x+∂τyz∂y+∂σz∂z−∂h∂z+γw−γr=0,where γr is the saturated bulk density of rock; h is the total water head.Combining the stress-strain constitutive equation and the geometric equation, the equilibrium equations expressed by the displacement component and the water head are(41)G∇2δx+G1−2υ∂εv∂x−∂h∂x+X0=0,G∇2δy+G1−2υ∂εv∂y−∂h∂y+Y0=0,G∇2δz+G1−2υ∂εv∂z−∂h∂z+γw−γr=0,where G=E/21+υ is the shear modulus; ∇2=∂2/∂x2+∂2/∂y2+∂2/∂z2 is the Laplace operator; εv=εx+εy+εz=∂δx/∂x+∂δy/∂y+∂δz/∂z is the volumetric strain; X0, Y0, and Z0 are equivalent volume forces caused by the initial strain ε0.The equilibrium equation of the whole unit matrix derived from equation (41) is(42)Keδe=F1e+F0e+K′ehe,where F1e is the joint force caused by external loads; F0e is the joint force caused by the initial strain; K′e is the coupling matrix of the seepage and stress.


## 5. Validation and Application

The Longyangxia hydropower station is the first large-scale reservoir in the upper reaches of the Yellow River, with area of 383 km2 and storage capacity of 24.7 billion m3. It is composed of main dam, auxiliary dam, and discharge structures (shown in Figure 9). Dam concrete began to be poured in June 1982, and the diversion tunnel was put down to store water on October 15, 1986. The foundation elevation is 2432 m, the dam crest elevation is 2610 m, the maximum dam height is 178 m, the dam length is 1227 m (including the main dam length of 396 m), the dam crest width is 23.5 m, and the maximum bottom width is 80 m.Figure 9
The Longyangxia hydropower station.The dam foundation is granodiorite with hard lithology. The dam abutment faults develop well through complex geotectonic movement and the main faults are shown in Figure10. Most of the NE-trending faults are filled with wide quartzite veins and form altered rock belts, among which the larger ones are the F120 fault and the A2 quartzite, located at the right bank. The F120 fault locates in the deep part near the dam abutment, with width of 2 m∼6 m and poor behavior, and its strike is nearly orthogonal to the thrust direction of the dam. The filled mylonite and breccia are not very continuous, but there is a fully and strongly weathered altered rock zone with an average thickness of 1.2 cm. The A2 quartzite has a pulse width of 5 m∼10 m. Its NE-trending vertical fractures and nearly horizontal fractures develop well, with strong water permeability. The F120 + A2 affected zone has several to tens of centimeters of intercalated mud containing breccia and detritus, which are dense and impervious to the leakage water. However, the fractured rock mass in the F120 + A2 zone constitutes a network-like seepage channel with strong permeability, and the seepage damage has always been in progress under the high hydraulic pressures.Figure 10
Distribution of main faults and layout of the groundwater monitoring holes.The reservoir water level change is shown in Figure11. Since 1986, the water level had never reached the design water level of 2600 m during 27 years of operation, but it rose continuously in the flood season of 2005 from 2558.70 m on April 26th to 2597.62 m on November 19th. As shown in Figure 11, the leakage at the F120 + A2 fault increased obviously and the monitoring leakage at 2463 m elevation increased by 252.7 mL/s during the water rising period.Figure 11
Measured time series of the leakage at 2463 m elevation of theF120 + A2 fault.


### 5.1. Validation by the 2D Model Simulation

An impervious curtain was formed in the fractured rock mass of the dam foundation after grouting, which could prolong the seepage path and reduce the seepage pressure. However, the curtain rock mass was gradually damaged under the long-term action of the stress field, the temperature field, and the seepage field. To simulate the situation, a 2-dimensional finite element model is firstly established, shown in Figure12, where the F120 + A2 fault passes through the gravity pier on the right bank (the foundation elevation is 2580 m). The model range is as follows: the length near the upstream gravity pier is 10 times of the curtain depth, the length near the downstream gravity pier is 3 times, and it is 10 times along the depth. Displacement-pore pressure quadrilateral coupling element is used for modelling. Convergence tests have been implemented to determine the number of elements and nodes: 57131 coupling elements and 57624 nodes in total.Figure 12
2D model and the zoning results of the dam foundation. (a) 2D finite element model; (b) divisions of theF120 + A2 fault.
(a)(b)Through geological prospecting, the permeability coefficient of the layered fault decreases along the depth. According to the GB 50487-2008 (code for engineering geological investigation of water resources and hydropower), the permeability coefficient range is10−4≤k<10−2 m/s. In this study, the F120 + A2 fault is divided into ten layers of rock 1– rock 10 (shown in Figure 12(b)) and the main parameters are listed in Table 1. The permeability coefficients are distributed as an exponential function with depth, as shown in Figure 13.Table 1
Material mechanical and seepage parameters.
ZoneElastic modulus (MPa)Poisson’s ratioIntensity (kg/m3)Permeability coefficient (m/s)Gravity pier2.8 × 1040.1672400k12 = 1.0 × 10−7Downstream rock mountain2.0 × 1040.222600k13 = 1.0 × 10−6Curtain2.8 × 1040.1672600k2 = 1.0 × 10−7The 1st rock layer5.4 × 1030.22600k1 = 2.64 × 10−4The 2nd rock layer5.4 × 1030.22600k3 = 2.32 × 10−4The 3rd rock layer5.4 × 1030.22600k4 = 2.04 × 10−4The 4th rock layer5.4 × 1030.22600k5 = 1.79 × 10−4The 5th rock layer5.4 × 1030.22600k6 = 1.58 × 10−4The 6th rock layer5.4 × 1030.22600k7 = 1.39 × 10−4The 7th rock layer5.4 × 1030.22600k8 = 1.22 × 10−4The 8th rock layer5.4 × 1030.22600k9 = 1.07 × 10−4The 9th rock layer5.4 × 1030.22600k10 = 9.42 × 10−5The 10th rock layer5.4 × 1030.22600k11 = 8.28 × 10−5Figure 13
Distribution of permeability coefficient of theF120 + A2 fault along the depth.The proposed methods are used to study the variation laws of the seepage field of the dam foundation. When the damage variables of theF120 + A2 fault at different depths are 0.10–0.55, the corresponding increment distributions of the permeability coefficient and the rock porosity are shown in Figure 14. They all increase gradually when the damage continues to develop, and they decrease along the rock depth, due to the rock mass confining pressures.Figure 14
Evolution of the permeability coefficients of the 1st–11th rock layer.Take 16 groups of calculating conditions for illustration, where the impervious curtain rock mass damage is 0.1–0.84, respectively, with the upstream normal water level of 2600 m and the downstream water level of 2450 m. According to equations (28) and (30), when different damage variables are taken for the curtain rock mass, corresponding curtain porosities and relative coefficients (the ratio of the curtain rock mass permeability coefficient to the first rock stratum permeability coefficient) are shown in Figure 15. The piezometer head values H2 at a typical point behind the curtain and corresponding osmotic pressure coefficient α2 both increase gradually. When the damage degree reaches 0.84, the seepage control of the curtain rock mass can be considered completely lost, since the relative coefficient equals 1.Figure 15
Damage evolution of the permeability characteristics of the curtain rock mass.The variations of the pore water pressure in the exposed local area of the downstream slope under different damage degrees are shown in Figure16. With the increase of the damage degree, the saturation line of the exposed local area gradually rises, indicating that the impervious curtain gradually weakens and also verifies the effectiveness of the proposed methods.Figure 16
Distributions of the pore water pressure of the fractured rock mass (unit: Pa). (a)D = 0.1, (b) D = 0.45, (c) D = 0.5, (d) D = 0.6, (e) D = 0.75, and (f) D = 0.84.
(a)(b)(c)(d)(e)(f)


### 5.2. Application on the 3D Model Simulation

According to the distribution of the main faults in Figure10, the 3D finite element model is established, as shown in Figure 17, including the main boundary range that affects the seepage behavior. The model range is 1 time the dam height on the left and right banks, 761 m in total, and the upstream and downstream length is 613 m. The model simulates the arch dam, the foundation rock mass, main faults, and antiseepage curtain, especially the F120 + A2 fault and the sandwiched mud. After convergence tests, there are 27,291 nodes and 24,631 elements.Figure 17
The 3D model of the Longyangxia dam.The simulation of the seepage field primarily considers the following boundary conditions. (1) Interception boundary of the regional basement: the foundation depth of the 3D model is 314.8 m, which is 2 times the curtain depth and it is treated as an impervious boundary condition. (2) The upstream and downstream interception boundary of the dam located area is treated as an impervious boundary. (3) Known boundary of the reservoir water level: the highest water level reached in 2005 is 2597.62 m upstream and 2453.52 m downstream at the corresponding time. The riverbed below the upstream and downstream water levels is regarded as the known boundary of constant water head. (4) Possible overflow boundaries are regarded as possible seepage boundaries, including the surfaces above the riverbed water level downstream and the intercepted boundary below the water level at the two banks. (5) Interception boundary of the two banks: the distribution of the groundwater is simulated with the measured data of groundwater level holes (Figure10) to obtain the calculation boundary. The water head on November 19, 2005, is selected for the simulation. The permeability coefficient values of each region are listed in Table 2.Table 2
Permeability coefficients of each region in the 3D model.
Material zoneShallow rock massDeep rock massThe curtainTheG4 faultTheF18 faultTheF73 faultPermeability coefficient (m/s)5.0 × 10−62.0 × 10−61.0 × 10−71.0 × 10−61.0 × 10−61.0 × 10−6In order to obtain the leakage through the F120 + A2 fault, the flow velocities of the elements corresponding to 2463 m, 2497 m, and 2530 m below the saturation lines after the curtain are extracted, and the element areas corresponding to the fault from 2463 m to 2497 m, from 2497 m to 2530 m, and from 2530 m to the saturation line are calculated. Figure 18(a) indicates the seepage extraction cross section through the F120 + A2 fault. The saturation lines at the highest water level of 2597.62 m are shown in Figure 18(b). The leakages passing the 2463 m elevation, the 2497 m elevation, and the 2530 m elevation of the fault are calculated. The results and corresponding measured data are listed in Table 3.Figure 18
TheF120 + A2 fault model and the calculated saturation lines. (a) Seepage extraction cross section; (b) the saturation line with water level of 2597.62 m.
(a)(b)Table 3
Calculation results of the leakages through theF120 + A2 fault (mL/s).
Storage date (water level)2463 m elevation2497 m elevation2530 m elevationCalculatedMeasuredCalculatedMeasuredCalculatedMeasured2005-6-4 (2575.45 m)322.49318.65359.35265.81224.79219.562005-8-18 (2580.85 m)359.28359.45377.98385.12238.82227.722005-10-5 (2591.54 m)431.51422.89567.36527.17290.61292.352005-11-16 (2597.62 m)470.22462.8584.83580.79315.96332.9The calculated results at typical elevations significantly increase during the rise of the reservoir water level and are in good agreement with the measured values. With the increase of the permeability coefficient, theF120 + A2 fault becomes the main leakage channel, and the grouted curtain has significant effect to reduce the seepage water head. The results also indicate that the water inflow rate is primarily affected by both the hydraulic conductivity and induced microcracks of the damage fault during the storage period. Once the curtain fails, the water level at each measured pore will become higher and the corresponding leakage amount will also increase obviously.According to the simulated and measured results, the seepage of the dam foundation on the right bank has a remarkable correlation with the change of the reservoir, but there has not been a persistent increase phenomenon. Namely, the seepage behavior has not deteriorated, and it is not necessary to repair the curtain temporarily. Because the high seepage pressures will gradually increase the damage of fractured rock mass, it is also suggested that the monitoring of the fault leakage should be paid special attention during the future operation, especially the high reservoir water period. If necessary, the drainage facilities of the dam foundation should be newly added to reduce the seepage pressure and ensure the safe operation of the dam.


## 5.1. Validation by the 2D Model Simulation

An impervious curtain was formed in the fractured rock mass of the dam foundation after grouting, which could prolong the seepage path and reduce the seepage pressure. However, the curtain rock mass was gradually damaged under the long-term action of the stress field, the temperature field, and the seepage field. To simulate the situation, a 2-dimensional finite element model is firstly established, shown in Figure12, where the F120 + A2 fault passes through the gravity pier on the right bank (the foundation elevation is 2580 m). The model range is as follows: the length near the upstream gravity pier is 10 times of the curtain depth, the length near the downstream gravity pier is 3 times, and it is 10 times along the depth. Displacement-pore pressure quadrilateral coupling element is used for modelling. Convergence tests have been implemented to determine the number of elements and nodes: 57131 coupling elements and 57624 nodes in total.Figure 12
2D model and the zoning results of the dam foundation. (a) 2D finite element model; (b) divisions of theF120 + A2 fault.
(a)(b)Through geological prospecting, the permeability coefficient of the layered fault decreases along the depth. According to the GB 50487-2008 (code for engineering geological investigation of water resources and hydropower), the permeability coefficient range is10−4≤k<10−2 m/s. In this study, the F120 + A2 fault is divided into ten layers of rock 1– rock 10 (shown in Figure 12(b)) and the main parameters are listed in Table 1. The permeability coefficients are distributed as an exponential function with depth, as shown in Figure 13.Table 1
Material mechanical and seepage parameters.
ZoneElastic modulus (MPa)Poisson’s ratioIntensity (kg/m3)Permeability coefficient (m/s)Gravity pier2.8 × 1040.1672400k12 = 1.0 × 10−7Downstream rock mountain2.0 × 1040.222600k13 = 1.0 × 10−6Curtain2.8 × 1040.1672600k2 = 1.0 × 10−7The 1st rock layer5.4 × 1030.22600k1 = 2.64 × 10−4The 2nd rock layer5.4 × 1030.22600k3 = 2.32 × 10−4The 3rd rock layer5.4 × 1030.22600k4 = 2.04 × 10−4The 4th rock layer5.4 × 1030.22600k5 = 1.79 × 10−4The 5th rock layer5.4 × 1030.22600k6 = 1.58 × 10−4The 6th rock layer5.4 × 1030.22600k7 = 1.39 × 10−4The 7th rock layer5.4 × 1030.22600k8 = 1.22 × 10−4The 8th rock layer5.4 × 1030.22600k9 = 1.07 × 10−4The 9th rock layer5.4 × 1030.22600k10 = 9.42 × 10−5The 10th rock layer5.4 × 1030.22600k11 = 8.28 × 10−5Figure 13
Distribution of permeability coefficient of theF120 + A2 fault along the depth.The proposed methods are used to study the variation laws of the seepage field of the dam foundation. When the damage variables of theF120 + A2 fault at different depths are 0.10–0.55, the corresponding increment distributions of the permeability coefficient and the rock porosity are shown in Figure 14. They all increase gradually when the damage continues to develop, and they decrease along the rock depth, due to the rock mass confining pressures.Figure 14
Evolution of the permeability coefficients of the 1st–11th rock layer.Take 16 groups of calculating conditions for illustration, where the impervious curtain rock mass damage is 0.1–0.84, respectively, with the upstream normal water level of 2600 m and the downstream water level of 2450 m. According to equations (28) and (30), when different damage variables are taken for the curtain rock mass, corresponding curtain porosities and relative coefficients (the ratio of the curtain rock mass permeability coefficient to the first rock stratum permeability coefficient) are shown in Figure 15. The piezometer head values H2 at a typical point behind the curtain and corresponding osmotic pressure coefficient α2 both increase gradually. When the damage degree reaches 0.84, the seepage control of the curtain rock mass can be considered completely lost, since the relative coefficient equals 1.Figure 15
Damage evolution of the permeability characteristics of the curtain rock mass.The variations of the pore water pressure in the exposed local area of the downstream slope under different damage degrees are shown in Figure16. With the increase of the damage degree, the saturation line of the exposed local area gradually rises, indicating that the impervious curtain gradually weakens and also verifies the effectiveness of the proposed methods.Figure 16
Distributions of the pore water pressure of the fractured rock mass (unit: Pa). (a)D = 0.1, (b) D = 0.45, (c) D = 0.5, (d) D = 0.6, (e) D = 0.75, and (f) D = 0.84.
(a)(b)(c)(d)(e)(f)


## 5.2. Application on the 3D Model Simulation

According to the distribution of the main faults in Figure10, the 3D finite element model is established, as shown in Figure 17, including the main boundary range that affects the seepage behavior. The model range is 1 time the dam height on the left and right banks, 761 m in total, and the upstream and downstream length is 613 m. The model simulates the arch dam, the foundation rock mass, main faults, and antiseepage curtain, especially the F120 + A2 fault and the sandwiched mud. After convergence tests, there are 27,291 nodes and 24,631 elements.Figure 17
The 3D model of the Longyangxia dam.The simulation of the seepage field primarily considers the following boundary conditions. (1) Interception boundary of the regional basement: the foundation depth of the 3D model is 314.8 m, which is 2 times the curtain depth and it is treated as an impervious boundary condition. (2) The upstream and downstream interception boundary of the dam located area is treated as an impervious boundary. (3) Known boundary of the reservoir water level: the highest water level reached in 2005 is 2597.62 m upstream and 2453.52 m downstream at the corresponding time. The riverbed below the upstream and downstream water levels is regarded as the known boundary of constant water head. (4) Possible overflow boundaries are regarded as possible seepage boundaries, including the surfaces above the riverbed water level downstream and the intercepted boundary below the water level at the two banks. (5) Interception boundary of the two banks: the distribution of the groundwater is simulated with the measured data of groundwater level holes (Figure10) to obtain the calculation boundary. The water head on November 19, 2005, is selected for the simulation. The permeability coefficient values of each region are listed in Table 2.Table 2
Permeability coefficients of each region in the 3D model.
Material zoneShallow rock massDeep rock massThe curtainTheG4 faultTheF18 faultTheF73 faultPermeability coefficient (m/s)5.0 × 10−62.0 × 10−61.0 × 10−71.0 × 10−61.0 × 10−61.0 × 10−6In order to obtain the leakage through the F120 + A2 fault, the flow velocities of the elements corresponding to 2463 m, 2497 m, and 2530 m below the saturation lines after the curtain are extracted, and the element areas corresponding to the fault from 2463 m to 2497 m, from 2497 m to 2530 m, and from 2530 m to the saturation line are calculated. Figure 18(a) indicates the seepage extraction cross section through the F120 + A2 fault. The saturation lines at the highest water level of 2597.62 m are shown in Figure 18(b). The leakages passing the 2463 m elevation, the 2497 m elevation, and the 2530 m elevation of the fault are calculated. The results and corresponding measured data are listed in Table 3.Figure 18
TheF120 + A2 fault model and the calculated saturation lines. (a) Seepage extraction cross section; (b) the saturation line with water level of 2597.62 m.
(a)(b)Table 3
Calculation results of the leakages through theF120 + A2 fault (mL/s).
Storage date (water level)2463 m elevation2497 m elevation2530 m elevationCalculatedMeasuredCalculatedMeasuredCalculatedMeasured2005-6-4 (2575.45 m)322.49318.65359.35265.81224.79219.562005-8-18 (2580.85 m)359.28359.45377.98385.12238.82227.722005-10-5 (2591.54 m)431.51422.89567.36527.17290.61292.352005-11-16 (2597.62 m)470.22462.8584.83580.79315.96332.9The calculated results at typical elevations significantly increase during the rise of the reservoir water level and are in good agreement with the measured values. With the increase of the permeability coefficient, theF120 + A2 fault becomes the main leakage channel, and the grouted curtain has significant effect to reduce the seepage water head. The results also indicate that the water inflow rate is primarily affected by both the hydraulic conductivity and induced microcracks of the damage fault during the storage period. Once the curtain fails, the water level at each measured pore will become higher and the corresponding leakage amount will also increase obviously.According to the simulated and measured results, the seepage of the dam foundation on the right bank has a remarkable correlation with the change of the reservoir, but there has not been a persistent increase phenomenon. Namely, the seepage behavior has not deteriorated, and it is not necessary to repair the curtain temporarily. Because the high seepage pressures will gradually increase the damage of fractured rock mass, it is also suggested that the monitoring of the fault leakage should be paid special attention during the future operation, especially the high reservoir water period. If necessary, the drainage facilities of the dam foundation should be newly added to reduce the seepage pressure and ensure the safe operation of the dam.


## 6. Conclusions

The mechanical characteristics of the fractured rock mass in seepage fields and the interaction mechanism of water-rock are critical for evaluation of the instability of dam foundations and even large-scale geological disasters. In this work, the creep damage and its effect on the permeability evolution are investigated combining with multiple test evidences. The permeability of the fractured rock mass can significantly increase with the growth and coalescence of microcracks under certain seepage pressures. Hence, the creep damage constitutive models are established to indicate the strength weakening under continuous loads. The evolution equation of the permeability coefficient is furtherly constructed with the crack tensor to characterize the mechanical behavior under the seepage-compression and shear action.The proposed approaches are applied to the Longyangxia dam and the simulated results agree well with the measured values, indicating that the presented model can describe the permeability variation of theF120 + A2 fault. It was further confirmed that the seepage pressures can intensify the crack expansion and penetration in the rock masses and may lead to the gradual failure of the foundation. As the leakage behavior of dam foundations is constantly changing due to the transformation of the seepage pressures, compared to preexisting models, the presented microcrack-based creep damage and seepage coupled method is capable of appraising the fracturing evolution and the leakage variation in the fractured rock mass of practical dam foundations. The research results can also provide a reference for engineering repair and supervision through controlling the permeability performance. For the future work, the criteria for judging the change of the seepage behavior of dam foundations and the early warning indicators need to be further studied for detecting the potential dangers in time.



---

*Source: 8832774-2021-02-13.xml*