# Equivalent Continuum Coupling-Based Slope Stability Analysis of Zhouning Pumped Storage Power Station

Authors: Jinhong Lin; Yongsheng Zhu; Jing Hou; Bing Pan; Weijiang Chu; Wenhong Huang
Journal: Geofluids
 (2021)

---

## Abstract

A 3D equivalent continuum coupling analysis of the slope of the Zhouning pumped storage power station was proposed in this study. Firstly, the hydraulic properties of rock mass, groundwater, and deformation mechanism under high-frequency fluctuations of the water level are analyzed. Secondly, the three-dimensional Digital Terrain Model (DTM) with surface image information was established by digital photogrammetry technology for cataloging field fracture information. Third, stress-seepage coupling simulation using a hybrid grid approach and the FLAC3D equivalent continuum model was implemented. The results show that the permeability of rock mass is stress-dependent and anisotropic, and the direction of maximum principal permeability of both bank slopes is approximately parallel to the river. Slope deformations induced by impoundment during the operation period are the superposition results of two mechanisms, referred to as “mattress effect” and “swelling effect.” Attention should be paid to the tensile strength of the rock mass during operation.


---

## Body

## 1. Introduction

Slope rock mass of water conservancy and hydropower projects will experience a complex stress-seepage coupling process due to reservoir water level variation during impoundment and operation period. While reducing the effective stress in the rock mass, water further deteriorates the strength condition of the rock mass by participating in the physical and chemical actions, which brings stability problems to the slope and hydraulic structure [1, 2]. At present, carrying out the experimental test, building a theoretical model, and performing numerical analysis are the main technologies to solve this kind of problem.As for the hydromechanical coupling analysis of fracture rock, there are three solutions, namely, the equivalent continuum seepage model, discrete fracture network seepage model, and continuum-discrete coupling seepage model. The first method simplifies the rock mass to a homogeneous porous media and is most widely used in engineering. The equivalent continuum model based on the Biot consolidation theory was first used to describe the stress-seepage coupling in soil. Since the cubic law for describing flow in a single fracture was proposed, the emergence of a discrete joint network simulation for rock mass structure description has promoted the development of a discontinuum model [3]. Barenblatt et al. [4] proposed a double-medium model to consider processes of non-steady-state infiltration in fissured rocks. In this model, besides considering fractures as infiltration channels, rock blocks can also play a role in water conduction. Due to the discontinuum model accurately describing the structure of rock mass, it is especially suitable for further study of the seepage and stress coupling mechanism from the perspective of the mesoscopic structure of the rock mass [5, 6]. However, the parameters used to describe the structure and permeability of the rock mass in the model are difficult to obtain directly by conventional means, and it is also time-consuming in data management for joint contact structure and detection in a large-scale model [7].Up to now, two numerical algorithms, namely, iterative coupling and full coupling, have been developed for stress-seepage coupling calculation. The former establishes the governing equations for the seepage field and the stress field, while the latter adopts the unified field governing equations. Because the seepage field and the stress field are obtained simultaneously in every numerical iteration, the fully coupling algorithm has higher computational efficiency. Wang et al. [8] proposed a four degrees of freedom fully coupling algorithm, which can also consider the permeability of the rock mass adjusted with stress. However, the fully coupling algorithm treats the rock mass as a linear elastic continuum and simplifies the influence of rock mass deformation on the seepage field. Therefore, the equivalent continuum model based on an iterative coupling algorithm is still the mainstream approach to study the coupling effect of seepage and stress in the rock mass, especially for the large-scale model study.The permeability characteristics of the rock mass mainly depend on the distribution and properties of fractures. It can be seen that the statistical analysis of rock mass fracture distribution and properties are the critical premises for performing stress-seepage coupling analysis. However, the researches of onsite fracture data acquisition and three-dimensional stress-seepage coupling analysis for study deformation and stability of the slope under the condition of water level variation are relatively few.Taking the slope of the lower reservoir of Zhouning Pumping and Reservoir Power Station in Fujian China as an example, this paper introduced 3D digital photogrammetry as an effective manner to collect the outcrop of the slope structural surface, then statistical analysis was implemented by using geological database technology to obtain fracture distribution parameters (orientation, spacing, etc.) to assist in determining rock mass permeability. A FLAC3D equivalent continuum model for stress-seepage analysis was built by using a hybrid grid approach proposed to have a higher computational efficiency. The model simulation was performed to investigate groundwater pressure, deformation behavior, and stability of the slope during the operation period with frequent fluctuation of the reservoir water level.


## 2. Project Overview and Rock Permeability

### 2.1. Project Overview

Fujian Zhouning pumped storage power station is located in Zhouning County, Ningde City, Fujian Province, China. The lower reservoir of the power station is located in Qibuxi valley, and the reservoir area is generally a long strip-shaped mountain canyon. The project is mainly composed of an upper reservoir, lower reservoir, and hydraulic structure including a water conveyance system, underground powerhouse, and switch station (Figure1). The designed capacity of the power station is 1200 MW.Figure 1
Topography and building layout of the Zhouning pumped storage power station.The bedrock strata are mainly the Nanyuan formation of the upper Jurassic of Mesozoic, the late Yanshanian intrusive rocks, and the Quaternary strata on the surface. The lithology is mainly potassium feldspar granite, and the rocks are composed of potassium feldspar, plagioclase, and other minerals with medium-fine grain granitic structure and massive structure. On both banks of the reservoir area, strongly to weakly weathered bedrock is widely exposed, and a completely weathered layer with limited thickness is exposed locally. There are fault F3 with crushed zone and several dikes (diabase dikes (βμ) and a quartz monzonite porphyry dike (ηoπ)) crossing through the lower reservoir of the power station, as shown in Figure 2.Figure 2
Plan view for engineering geology at lower reservoir of the power station.


### 2.2. Permeability Characteristics of Rock

Packer test were carried out in field boreholes, and a total of 287 test results were obtained in the lower reservoir area of the Zhouning pumped storage power station. The permeability characteristics of the rock mass were stratified and statistically analyzed according to different weathering degrees (Figure3). In strongly weathered rock mass, the medium permeable section (10Lu≤q<100Lu) accounts for nearly 20% of the total test section length and the weak permeable section (1Lu≤q<10Lu) accounts for more than 80%. With the decrease of the weathering degree, the length proportion of the medium permeability section gradually decreases, but the length proportion of micropermeability (q<1Lu) increases. In general, the permeability of the rock mass decreases with the decrease of the weathering degree, which is dominated by weak micropermeability.Figure 3
Permeability of granite in the Zhouning pumped storage power station.The permeability characteristics of the rock mass depend on not only the structural conditions of the rock mass but also the in situ stress state. The integrity of the slope rock mass increases with the increase of the buried depth or stress; meanwhile, closure of the fracture decreases resulting in the fracture becoming tighter due to stress compression, which leads to the decrease of permeability characteristics of rock mass [9, 10].A fracture is the main component of the rock mass structure, to which the permeability of the rock mass closely correlates. Geometric characteristics of fractures (mainly shape, size, width, and connectivity), statistical distribution, and physical and mechanical properties of rock walls or fillings can be objective factors affecting fracture property and behavior including permeability characteristics [11]. At the same time, the permeability characteristics of the rock mass can also have complex properties such as anisotropy and size effect due to the group, dominance, and uncertainty of fracture distribution [12].Therefore, the fractures in the slope rock mass play a vital role in the stability of the slope and engineering structures during the operation period. The successful acquisition of fracture distribution in the field and determination of the permeability characteristic parameters are critical concerns in stress-seepage coupling analysis.


## 2.1. Project Overview

Fujian Zhouning pumped storage power station is located in Zhouning County, Ningde City, Fujian Province, China. The lower reservoir of the power station is located in Qibuxi valley, and the reservoir area is generally a long strip-shaped mountain canyon. The project is mainly composed of an upper reservoir, lower reservoir, and hydraulic structure including a water conveyance system, underground powerhouse, and switch station (Figure1). The designed capacity of the power station is 1200 MW.Figure 1
Topography and building layout of the Zhouning pumped storage power station.The bedrock strata are mainly the Nanyuan formation of the upper Jurassic of Mesozoic, the late Yanshanian intrusive rocks, and the Quaternary strata on the surface. The lithology is mainly potassium feldspar granite, and the rocks are composed of potassium feldspar, plagioclase, and other minerals with medium-fine grain granitic structure and massive structure. On both banks of the reservoir area, strongly to weakly weathered bedrock is widely exposed, and a completely weathered layer with limited thickness is exposed locally. There are fault F3 with crushed zone and several dikes (diabase dikes (βμ) and a quartz monzonite porphyry dike (ηoπ)) crossing through the lower reservoir of the power station, as shown in Figure 2.Figure 2
Plan view for engineering geology at lower reservoir of the power station.


## 2.2. Permeability Characteristics of Rock

Packer test were carried out in field boreholes, and a total of 287 test results were obtained in the lower reservoir area of the Zhouning pumped storage power station. The permeability characteristics of the rock mass were stratified and statistically analyzed according to different weathering degrees (Figure3). In strongly weathered rock mass, the medium permeable section (10Lu≤q<100Lu) accounts for nearly 20% of the total test section length and the weak permeable section (1Lu≤q<10Lu) accounts for more than 80%. With the decrease of the weathering degree, the length proportion of the medium permeability section gradually decreases, but the length proportion of micropermeability (q<1Lu) increases. In general, the permeability of the rock mass decreases with the decrease of the weathering degree, which is dominated by weak micropermeability.Figure 3
Permeability of granite in the Zhouning pumped storage power station.The permeability characteristics of the rock mass depend on not only the structural conditions of the rock mass but also the in situ stress state. The integrity of the slope rock mass increases with the increase of the buried depth or stress; meanwhile, closure of the fracture decreases resulting in the fracture becoming tighter due to stress compression, which leads to the decrease of permeability characteristics of rock mass [9, 10].A fracture is the main component of the rock mass structure, to which the permeability of the rock mass closely correlates. Geometric characteristics of fractures (mainly shape, size, width, and connectivity), statistical distribution, and physical and mechanical properties of rock walls or fillings can be objective factors affecting fracture property and behavior including permeability characteristics [11]. At the same time, the permeability characteristics of the rock mass can also have complex properties such as anisotropy and size effect due to the group, dominance, and uncertainty of fracture distribution [12].Therefore, the fractures in the slope rock mass play a vital role in the stability of the slope and engineering structures during the operation period. The successful acquisition of fracture distribution in the field and determination of the permeability characteristic parameters are critical concerns in stress-seepage coupling analysis.


## 3. Parametric Analysis for Equivalent Continuum Stress-Seepage Coupling Model

### 3.1. Biot Theory

In the finite difference program FLAC3D, an equivalent continuum model for studying stress-seepage coupling regards rock or soil as continuous porous media, and the classical Darcy law is used to describe the flow process of the fluid. At the same time, the Biot principle is satisfied; that is, the rock and soil skeleton and mineral particles are compressible. The fluid seepage and its coupling with stress satisfy the following main governing equations:
(1)Flowequation:qi=−kilk^sp−ρwxjgj,l,(2)Equilibriumequation:1M∂p∂t+ns∂s∂t=1s∂ζ∂t−a∂ε∂t,where qi is the flow rate per unit area; kil, k^s are the permeability coefficient tensor and the associated coefficient reflecting the influence of saturation on permeability, respectively; kil is correlated to the tensor of permeability coefficient Kil and satisfies kil=Kil/ρwgi (ρw and gi are fluid density and gravitational acceleration, respectively); p is the pore pressure; ε is the porosity of the medium and the volumetric strain of the solid skeleton; ζ is the imbalance of the fluid volume in porous medium per unit volume; and M and α are the Biot modulus and Biot coefficient (dimensionless and in the range between 0 and 1), respectively.


### 3.2. Hydraulic and Mechanical Properties of Rock Mass

Snow [13], Romm [14], and Oda [15] proposed the expression of the permeability tensor of fractured rock mass successively, and the commonly used description form is as follows:
(3)K=g12vw∑i=1Nbi3Si1−nxi2−nxinyi−nzinxi−nxinyi1−nyi2nyinzi−nzinxinyinzi1−nzi2,where g and vw are the gravitational acceleration and the fluid viscosity coefficient, respectively; N is the number of fracture sets; bi and Si are the averaged closure and spacing in the i-th fracture set, respectively; and nxi, nyi, and nzi are the averaged cosine components of a normal vector in the i-th fracture set.In Equation (1), it is assumed that the fractures are completely through the rock mass and the filling is not considered so that the permeability coefficient of the rock mass is in general overestimated. It is suggested to determine a more reasonable modified principal permeability coefficient tensor Km by further correcting the principal permeability coefficient tensor KT of K using the packer tests in borehole results [16], as represented in the following expression:
(4)Km=mKTmK1000mK2000mK3,where K1, K2, and K3 are the principal permeability coefficients; m is the modified coefficient and equals to m=k1k2⋯knn/K1K1K3, in which k1, k2,…, kn are the permeability data points of packer test results. Note that the orientation of the principal permeability coefficient remains unchanged after modification.In addition to drilling, adit, and other conventional exploration, 3DM digital photogrammetry technology [17] is introduced to collect fracture outcrops on the slope rock mass by each weathered zone, in accordance to which the slope permeability parameters are comprehensively determined. The principle of 3DM digital photogrammetry technology is finding the same objects in multiple images, then connecting the perspective center and target points and estimating the associated coordinates in 3D space by projection, finally reproducing a digital 3D terrain model (DTM) and elevation model (DEM) containing the aerial image. The spatial location, length, and even closure of each fracture can be acquired by geometric logging and analysis tools provided in 3DM software (Figure 4).Figure 4
3D digital model and data extraction of fracture outcrop cataloging in 3DM software.
(a)
Slope model of the left bank
(b)
Slope model of the right bankThe fracture logging data were summarized into the geological database management and analysis system GeoDataManage [18], in which dominant joint set number and distribution parameter (average normal cosine and spacing) for each fracture set were statistically analyzed for the different weathered zones of the slope rock mass at both banks. In particular, the average closure of each fracture set is quantified according to the relationship of geological description (closed, micro open, and open) to fracture closure established by Peng [19]. The modified equivalent principal permeability coefficients Ki (Table 1) are obtained by firstly using these known parameters in Equation (1) and then further corrected by packer test results in boreholes.Table 1
Principal permeability coefficient of slopeKi (10-7 m/s).
Rock mass
Left bank
Right bank
Weakly weathered upper
K
1
=
1.81, K2=9.31, K3=9.98
K
1
=
1.87, K2=5.58, K3=7.40
Weakly weathered lower
K
1
=
0.58, K2=4.59, K3=4.87
K
1
=
0.92, K2=2.64, K3=3.54
Fresh/microweathered
K
1
=
0.28, K2=2.21, K3=2.35
K
1
=
0.45, K2=1.30, K3=1.73It can be seen that the rock permeability characteristics on both bank slopes are varied due to differing fracture distributions, and the permeability characteristics of the rock mass in each weathered zone are different and anisotropic. Note that theK3 is the maximum principal permeability coefficient and the K1 is the minor.It was found that the orientations of the principal permeability coefficient of the slope rock mass on the same bank were almost the same. Therefore, the orientations of permeability coefficients are adopted on both sides identically from the orientations in the upper weakly weathered rock section, as expressed in the following matrices:(5)Vleft=−0.02440.7421−0.6699−0.01630.66970.74240.99960.0290−0.0043,Vright=0.8203−0.1088−0.56160.4412−0.50440.74220.36400.85660.3657.The azimuth conversion of the permeability tensor shows that the plunge ofK3 of the left bank and right bank are nearly horizontal and close to 30°, respectively. The angles between the trend of K3 and the strike direction of the slope are in the range of 10°-20° at both banks. It can be seen that the medium and steep fracture sets along the river play a controlling role in rock mass permeability.Since Bear [12] proposed the concept of representative elementary volume (REV), the REV has been widely applied to evaluate the applicability of the equivalent continuum model to describe the mechanical properties and the behavior of the rock mass from the perspective of size effect. Maini et al. [20] and Winson and Witherspoon [21] illustrated the REV of the rock mass by using the ratio of the averaged distance between fractures and model size, which was supposed to be smaller than 1/20 and 1/50 as a judgment for the applicability of the equivalent continuum model. The spacing of fractures in the Zhouning project is far less than the slope scale, which means the equivalent continuum model can be applied to study the seepage-stress coupling in the rock mass.The constitutive parameters of Mohr-Coulomb (Table2) are estimated based on the results of field exploration, laboratory tests, and engineering classification of the rock mass. It should be stated the maximum width of the fault F3 is up to 2.2 m, the rock mass in F3 is composed of ferromanganese rendering cataclasite with calcareous and a small amount of fault gouge locally, and the thickness of the influence belt of fault F3 is approximately 1.5 m. From a point of view of unfavorable consideration, the mechanical properties of the influence belt were chosen to be consistent with F3. It should also be noted that the softening behavior of the rock mass due to groundwater immersion was not considered in this study.Table 2
Physical and mechanical parameters of rock masses.
Geological units and rock mass classification
Weathered state
Nature weightγ (kN/m3)
Elastic modulus
Poisson’s ratiov
Eff. friction coefficient/cohesion
Porosity
E (GPa)
f
′
C
′ (MPa)
n (%)
K-feldspar granite
III2A
Weakly weathered (upper)
25.6
5.00
0.28
0.95
1.100
7.5
K-feldspar granite
III1A
Weakly weathered (lower)
25.5
12.00
0.26
1.20
1.300
5.0
K-feldspar granite
IIA
Fresh/weakly weathered
25.8
15.00
0.24
1.30
1.500
2.0
Dykes
Completely weathered
19.5
0.04
0.32
0.47
0.018
10.0
Weakly weathered
27.6
5.00
0.28
1.00
1.100
6.0
Weakly weathered
27.5
9.00
0.24
1.20
1.200
3.0
Fault F3 and effecting zones
25.0
1.10
0.40
0.50
0.150
20.0The permeability coefficients in dikes are1.0×10−6m/s for completely the weathered, 4.0×10−7m/s for weakly weathered, and 2.0×10−7m/s for slightly weathered. The permeability coefficient of fault F3 and its influence belt is united at 2.0×10−6m/s. It should be noted that the permeability of these geological structures is also orthotropic, and the permeability coefficient along the normal direction to the structural plane is adopted as 1/10 of the above-mentioned permeability coefficients.


### 3.3. Variation of Reservoir Water Level during Operation

During the operation period of the power station, the water level of the lower reservoir takes one day (24 h) as an operating cycle and will sequentially go through three stages. Specifically, the reservoir water level increases from the dead water level to the storage water level firstly then decreases once again to the dead water level (stable state) and keeps unchanged in the third stage. It is represented in Figure5 that the duration of different stages is 6 h, 9 h, and 9 h.Figure 5
Variation characteristics of reservoir water level in one cycle (1 day/d) during the operating period.


## 3.1. Biot Theory

In the finite difference program FLAC3D, an equivalent continuum model for studying stress-seepage coupling regards rock or soil as continuous porous media, and the classical Darcy law is used to describe the flow process of the fluid. At the same time, the Biot principle is satisfied; that is, the rock and soil skeleton and mineral particles are compressible. The fluid seepage and its coupling with stress satisfy the following main governing equations:
(1)Flowequation:qi=−kilk^sp−ρwxjgj,l,(2)Equilibriumequation:1M∂p∂t+ns∂s∂t=1s∂ζ∂t−a∂ε∂t,where qi is the flow rate per unit area; kil, k^s are the permeability coefficient tensor and the associated coefficient reflecting the influence of saturation on permeability, respectively; kil is correlated to the tensor of permeability coefficient Kil and satisfies kil=Kil/ρwgi (ρw and gi are fluid density and gravitational acceleration, respectively); p is the pore pressure; ε is the porosity of the medium and the volumetric strain of the solid skeleton; ζ is the imbalance of the fluid volume in porous medium per unit volume; and M and α are the Biot modulus and Biot coefficient (dimensionless and in the range between 0 and 1), respectively.


## 3.2. Hydraulic and Mechanical Properties of Rock Mass

Snow [13], Romm [14], and Oda [15] proposed the expression of the permeability tensor of fractured rock mass successively, and the commonly used description form is as follows:
(3)K=g12vw∑i=1Nbi3Si1−nxi2−nxinyi−nzinxi−nxinyi1−nyi2nyinzi−nzinxinyinzi1−nzi2,where g and vw are the gravitational acceleration and the fluid viscosity coefficient, respectively; N is the number of fracture sets; bi and Si are the averaged closure and spacing in the i-th fracture set, respectively; and nxi, nyi, and nzi are the averaged cosine components of a normal vector in the i-th fracture set.In Equation (1), it is assumed that the fractures are completely through the rock mass and the filling is not considered so that the permeability coefficient of the rock mass is in general overestimated. It is suggested to determine a more reasonable modified principal permeability coefficient tensor Km by further correcting the principal permeability coefficient tensor KT of K using the packer tests in borehole results [16], as represented in the following expression:
(4)Km=mKTmK1000mK2000mK3,where K1, K2, and K3 are the principal permeability coefficients; m is the modified coefficient and equals to m=k1k2⋯knn/K1K1K3, in which k1, k2,…, kn are the permeability data points of packer test results. Note that the orientation of the principal permeability coefficient remains unchanged after modification.In addition to drilling, adit, and other conventional exploration, 3DM digital photogrammetry technology [17] is introduced to collect fracture outcrops on the slope rock mass by each weathered zone, in accordance to which the slope permeability parameters are comprehensively determined. The principle of 3DM digital photogrammetry technology is finding the same objects in multiple images, then connecting the perspective center and target points and estimating the associated coordinates in 3D space by projection, finally reproducing a digital 3D terrain model (DTM) and elevation model (DEM) containing the aerial image. The spatial location, length, and even closure of each fracture can be acquired by geometric logging and analysis tools provided in 3DM software (Figure 4).Figure 4
3D digital model and data extraction of fracture outcrop cataloging in 3DM software.
(a)
Slope model of the left bank
(b)
Slope model of the right bankThe fracture logging data were summarized into the geological database management and analysis system GeoDataManage [18], in which dominant joint set number and distribution parameter (average normal cosine and spacing) for each fracture set were statistically analyzed for the different weathered zones of the slope rock mass at both banks. In particular, the average closure of each fracture set is quantified according to the relationship of geological description (closed, micro open, and open) to fracture closure established by Peng [19]. The modified equivalent principal permeability coefficients Ki (Table 1) are obtained by firstly using these known parameters in Equation (1) and then further corrected by packer test results in boreholes.Table 1
Principal permeability coefficient of slopeKi (10-7 m/s).
Rock mass
Left bank
Right bank
Weakly weathered upper
K
1
=
1.81, K2=9.31, K3=9.98
K
1
=
1.87, K2=5.58, K3=7.40
Weakly weathered lower
K
1
=
0.58, K2=4.59, K3=4.87
K
1
=
0.92, K2=2.64, K3=3.54
Fresh/microweathered
K
1
=
0.28, K2=2.21, K3=2.35
K
1
=
0.45, K2=1.30, K3=1.73It can be seen that the rock permeability characteristics on both bank slopes are varied due to differing fracture distributions, and the permeability characteristics of the rock mass in each weathered zone are different and anisotropic. Note that theK3 is the maximum principal permeability coefficient and the K1 is the minor.It was found that the orientations of the principal permeability coefficient of the slope rock mass on the same bank were almost the same. Therefore, the orientations of permeability coefficients are adopted on both sides identically from the orientations in the upper weakly weathered rock section, as expressed in the following matrices:(5)Vleft=−0.02440.7421−0.6699−0.01630.66970.74240.99960.0290−0.0043,Vright=0.8203−0.1088−0.56160.4412−0.50440.74220.36400.85660.3657.The azimuth conversion of the permeability tensor shows that the plunge ofK3 of the left bank and right bank are nearly horizontal and close to 30°, respectively. The angles between the trend of K3 and the strike direction of the slope are in the range of 10°-20° at both banks. It can be seen that the medium and steep fracture sets along the river play a controlling role in rock mass permeability.Since Bear [12] proposed the concept of representative elementary volume (REV), the REV has been widely applied to evaluate the applicability of the equivalent continuum model to describe the mechanical properties and the behavior of the rock mass from the perspective of size effect. Maini et al. [20] and Winson and Witherspoon [21] illustrated the REV of the rock mass by using the ratio of the averaged distance between fractures and model size, which was supposed to be smaller than 1/20 and 1/50 as a judgment for the applicability of the equivalent continuum model. The spacing of fractures in the Zhouning project is far less than the slope scale, which means the equivalent continuum model can be applied to study the seepage-stress coupling in the rock mass.The constitutive parameters of Mohr-Coulomb (Table2) are estimated based on the results of field exploration, laboratory tests, and engineering classification of the rock mass. It should be stated the maximum width of the fault F3 is up to 2.2 m, the rock mass in F3 is composed of ferromanganese rendering cataclasite with calcareous and a small amount of fault gouge locally, and the thickness of the influence belt of fault F3 is approximately 1.5 m. From a point of view of unfavorable consideration, the mechanical properties of the influence belt were chosen to be consistent with F3. It should also be noted that the softening behavior of the rock mass due to groundwater immersion was not considered in this study.Table 2
Physical and mechanical parameters of rock masses.
Geological units and rock mass classification
Weathered state
Nature weightγ (kN/m3)
Elastic modulus
Poisson’s ratiov
Eff. friction coefficient/cohesion
Porosity
E (GPa)
f
′
C
′ (MPa)
n (%)
K-feldspar granite
III2A
Weakly weathered (upper)
25.6
5.00
0.28
0.95
1.100
7.5
K-feldspar granite
III1A
Weakly weathered (lower)
25.5
12.00
0.26
1.20
1.300
5.0
K-feldspar granite
IIA
Fresh/weakly weathered
25.8
15.00
0.24
1.30
1.500
2.0
Dykes
Completely weathered
19.5
0.04
0.32
0.47
0.018
10.0
Weakly weathered
27.6
5.00
0.28
1.00
1.100
6.0
Weakly weathered
27.5
9.00
0.24
1.20
1.200
3.0
Fault F3 and effecting zones
25.0
1.10
0.40
0.50
0.150
20.0The permeability coefficients in dikes are1.0×10−6m/s for completely the weathered, 4.0×10−7m/s for weakly weathered, and 2.0×10−7m/s for slightly weathered. The permeability coefficient of fault F3 and its influence belt is united at 2.0×10−6m/s. It should be noted that the permeability of these geological structures is also orthotropic, and the permeability coefficient along the normal direction to the structural plane is adopted as 1/10 of the above-mentioned permeability coefficients.


## 3.3. Variation of Reservoir Water Level during Operation

During the operation period of the power station, the water level of the lower reservoir takes one day (24 h) as an operating cycle and will sequentially go through three stages. Specifically, the reservoir water level increases from the dead water level to the storage water level firstly then decreases once again to the dead water level (stable state) and keeps unchanged in the third stage. It is represented in Figure5 that the duration of different stages is 6 h, 9 h, and 9 h.Figure 5
Variation characteristics of reservoir water level in one cycle (1 day/d) during the operating period.


## 4. 3D Seepage-Stress Coupling Model

Same as thermal and dynamic problems, when the transient seepage-stress coupling process is iteratively solved in a time domain using the finite difference method, the critical time step of iterationΔtcr is closely related to the mesh quality in the model, and even the computational efficiency is, in general, lower than static analysis remarkably [22]. The Δtcr is defined in the seepage-stress coupling model as expressed
(6)Δtcr=1a⋅Lc2c,where the diffusion coefficient of seepage c can be expressed as k/1/M+α2/K+4G/3; the k and G are the permeability coefficient and shear modulus of rock, respectively; a is the index representing mesh quality, which equals 6 if the model is united and structurally gridded; and Lc is the characteristic length of the model, commonly adopted as the grid length.It can be seen that the computational efficiency of the seepage-stress coupling model is affected by the shape and length of the grid and the grid length has a critical weight of influence. Note thatΔtcr∝1/Lc2.The structured grid, unstructured grid, and octree grid are supported in program FLAC3D (Figure 6). Different from the unstructured grid, every 8 grids (hexahedral grid) share one grid point in a structured grid. The octree grid is one kind of unstructured grid, which can be divided into 8 hexahedral subgrids as needed. The numerical element Attach in FLAC3D with an internal interpolation function needs to be introduced to ensure continuity of mechanical response at grid discontinuities in a large-scale complicated model.Figure 6
Typical grid types in FLAC3D.The above-mentioned various types of grid andAttach element are applied in the seepage-stress coupling model for improving the computational accuracy and efficiency. In detail, the slope is gridded by using a structured grid and octree grid, and the dam is gridded by using an unstructured grid. By assuming that the dam is impermeable in the analysis, its grid type and size do not affect the critical time step Δtcr. However, this combination of grid types will bring discontinuities in the model; the continuity of the mechanical response of the model is further improved by introducing the Attach element at these discontinuities as shown in Figure 7.Figure 7
Model mesh settings for foundation base of the dam.The survey scope according to the geological data is650m×600m. Depending on the sensitivity analysis for the effects of the grid length on the computational efficiency, the time step Δtcr of 16 s and 36 s corresponds to the grid length of 2.5 m and 4.5 m, respectively. The model built with an element length of 4.5 m is shown in Figure 8. It is noted that the model is composed of 300 thousand structured elements, 310 thousand grid points, and 15 thousand Attach elements. The constitutive model and associated boundary conditions are assigned:
(1)
For stress-strain analysis: the Mohr-Coulomb model is chosen for performing mechanical analysis, and normal displacement constraints are set at the bottom and side of the model(2)
For seepage analysis: an anisotropic equivalent continuum seepage model is adopted to the model seepage behavior of the rock mass, and varying hydraulic boundary conditions in the time domain are defined by the built-in development language FISH in FLAC3DFigure 8
Seepage-stress coupling analysis model of slope in the lower reservoir area of the Zhouning hydropower station.The geological characteristics in the slope on both sides:(1)
A large area of weakly weathered rock in a shallow layer of the slope, and a small area of completely weathered rock exposed merely above the storage water level in this layer(2)
Faults such as the fault F3 lying under the dam are usually detrimental to slope stability; however, the compression on the F3 due to the self-weight of the dam is beneficial


## 5. Hydraulic and Mechanical Behavior in the Slope Rock Mass

### 5.1. Slope Groundwater and Deformation Analysis in One Operation Cycle

Simulation of seepage-stress coupling in the slope during the operation period was performed. This section is aimed at illustrating the variation regulation of groundwater pressure and deformations in the slope rock mass due to water level change during one operation cycle (1 day/d).The groundwater pressure distributions in the representative geological section of the upstream slope (Figure9) are demonstrated at the key time points A, B, C, and D (Figure 4). The calculation results indicate that hysteresis of groundwater pressure adjustment in the slope rock mass is the critical characteristic with the variation of the reservoir water level in one operation cycle. The hysteresis characteristic in phases A-B can be observed while increasing the water level of the reservoir, i.e., infiltration rate of the reservoir water into the slope lags behind the lifting speed of the reservoir water due to low permeability of the rock mass, resulting in an upturned parabolic shape of the ground pressure contour in the slope. The observation of hysteresis as a result of outflow speed from the slope lower than the reservoir waterfall speed in phases B-C is also recorded.Figure 9
Distribution and evolution of pore pressure in slope.Besides being the hydraulic boundary, reservoir water also acts as a ponded load on the slope to affect the deformation behavior of the rock mass.It is represented in Figure10 that the deformations of the slope induced by reservoir water level lifting in phases A-B mainly occur in the rock mass below the normal storage water level and with the maximum magnitude of 2 mm. It is also observed that dykes and fault F3 have no obvious influence on slope deformation behavior. Note that the induced deformations pointing into the slope are a kind of compressive behavior, which is particularly different from deformations in the slope composed of loose rock usually pointing to the outside of the slope during water level rising. The compressive deformations are reduced in phases B-C as the result of a rapid decrease in water level and groundwater pressure adjustment in the slope. It should be mentioned that, at the end of one operation cycle (time D in Figure 10), the slope deformations finally point to the outside of the slope and with the maximum magnitude smaller than 1 mm. This deformation behavior usually means the slope stability weakens.Figure 10
Distribution and evolution of deformations in slope.The groundwater pressure and deformation behavior in the slope rock mass during the operation period are dominated by hydraulic-mechanical coupling or referred to as the stress-seepage coupling.Owing to hysteresis of seepage, the excess groundwater pressure appears in the riverbed rock mass in phases A-B for the fluid in rock mass compressed by an increase of ponded water load and will experience dissipation during the decease of the reservoir water level in phases B-C and C-D, as represented in Figure9. Many studies have shown that the hysteresis characteristic of seepage mainly depends on the permeability of the rock and the variation speed of the reservoir water level. The higher permeability of the rock and the lower change speed of the water level cause short-term hysteretic duration. It needs to be mentioned that the distributions of groundwater pressure are almost identical at the time of D and A in Figure 9, which indicates a critical phenomenon; the seepage-stress coupling process would reach the steady-state in one operation period which may be due to a longer time with unchanged water level in phases C-D. In addition, the evolution process of the groundwater pressure in Figure 9 also shows that the depth affected by water level fluctuation varies up to 40 m in the left slope, 45 m in the right slope, and 30 m in the riverbed.


### 5.2. Decoupling Solution for Stress-Seepage Coupling Simulation

The characteristic times of stress propagation and fluid diffusion in equivalent continuum such as united porous media are defined as(7)Stresspropagation:tcm=ρdKu+4/3GLc,Fluiddiffusion:tcf=Lc2c,where Ku is the undrained bulk modulus of the medium.The comparison between the above-mentioned two characteristic times illustrates the coupling degree of stress and seepage, which would help choose the proper approach with high computational efficiency for stress-seepage coupling calculation, i.e., high coupling degree of seepage stress is represented with identical magnitudes of the characteristic times. In general,tcm is less than tcf, while tcm≫tcf and tcm≪tcf, respectively, describe two extreme permeability properties of rock or soil, i.e., completely drained and undrained. Objective to these two extreme situations, the decoupling approach can be introduced to solve the stress-seepage coupling problem, in which the seepage and stress field are solved separately. Although the decoupling approach simplifies the coupling process, it is widely used in reality for obtaining higher computational efficiency.The decoupling approach is also available in this study fortcf more than 500 times to tcm according to the rock properties of the Zhouning project. The decoupling approach is designed to be available to analyze stress-seepage coupling for every typical reservoir water varying phase, i.e., A-B, B-C, and C-D in Figure 5. In one phase of decoupling modeling, the computation needs to be decomposed into two main stages: (1) perform stress calculation with static mode, in which groundwater is not allowed to flow; (2) implement seepage simulation in the real-time domain to have the stress field disturbed by groundwater diffusion.A conceptual plane model as shown in Figure11 was built aimed at illustrating the decoupling approach and further investigating groundwater pressure and deformation behavior in the slope during the impoundment period. In the conceptual model, the rock mass is considered as an elastic material to ignore its nonlinearity for highlighting the stress-seepage coupling mechanism, and the permeability and mechanical parameters are consistent with the upper weakly weathered rock mass of Zhouning slope (Class III). In addition, an idealized steep topography was designed to deform the behavior at the crest of the slope.Figure 11
Deformation mechanism of the slope during impoundment.
(a)
Mattress effect before seepage
(b)
Swelling effect after seepageAccording to Figure11, the decoupling computation results indicate that there are two successive critical effects that exist in the phase of reservoir impoundment, referred to as “mattress effect” and “swelling effect” in this paper. The impoundment was modeled as a sudden water level increase from the dead level to the normal storage level, meanwhile, ignoring the groundwater diffusion for low permeability of the slope rock mass; consequently, the mattress effect occurs, in which compressive deformation occurs under the water level and tensile deformation above the water level. The excess groundwater pressure in the riverbed is also observed as shown in A of Figure 10. After that, seepage simulation is implemented in the real-time domain; consequently, the swelling effect is observed, in which the whole model is uplifted as a result of effective stress reduction due to groundwater diffusion in the slope rock mass. It should be noted that the “swelling effect” is closely related to the time for the groundwater to diffuse; the longer the groundwater diffusion time, the more obvious the “swelling effect” until a state of equilibrium.Objectively, the two mechanisms of “mattress effect” and “swelling effect” run through the whole process of slope stress-seepage coupling dynamically and synchronously until the disturbance reaches an equilibrium state. The slope deformation during the impoundment period is the result of the superposition of the two effects.Decoupling analysis reveals the mechanism of seepage-stress in the rock mass in a manner of simplification and is identical to the analysis of undrained consolidation settlement.The decrease of the water level is the inverse process of impoundment, in which the ponded water load is gradually removed, and the seepage force changes to point to the outside of the slope. At this time, besides the hydraulic behavior, the slope deformation and stability are closely related to the mechanical properties of the rock mass, which is further investigated in a further section.


## 5.1. Slope Groundwater and Deformation Analysis in One Operation Cycle

Simulation of seepage-stress coupling in the slope during the operation period was performed. This section is aimed at illustrating the variation regulation of groundwater pressure and deformations in the slope rock mass due to water level change during one operation cycle (1 day/d).The groundwater pressure distributions in the representative geological section of the upstream slope (Figure9) are demonstrated at the key time points A, B, C, and D (Figure 4). The calculation results indicate that hysteresis of groundwater pressure adjustment in the slope rock mass is the critical characteristic with the variation of the reservoir water level in one operation cycle. The hysteresis characteristic in phases A-B can be observed while increasing the water level of the reservoir, i.e., infiltration rate of the reservoir water into the slope lags behind the lifting speed of the reservoir water due to low permeability of the rock mass, resulting in an upturned parabolic shape of the ground pressure contour in the slope. The observation of hysteresis as a result of outflow speed from the slope lower than the reservoir waterfall speed in phases B-C is also recorded.Figure 9
Distribution and evolution of pore pressure in slope.Besides being the hydraulic boundary, reservoir water also acts as a ponded load on the slope to affect the deformation behavior of the rock mass.It is represented in Figure10 that the deformations of the slope induced by reservoir water level lifting in phases A-B mainly occur in the rock mass below the normal storage water level and with the maximum magnitude of 2 mm. It is also observed that dykes and fault F3 have no obvious influence on slope deformation behavior. Note that the induced deformations pointing into the slope are a kind of compressive behavior, which is particularly different from deformations in the slope composed of loose rock usually pointing to the outside of the slope during water level rising. The compressive deformations are reduced in phases B-C as the result of a rapid decrease in water level and groundwater pressure adjustment in the slope. It should be mentioned that, at the end of one operation cycle (time D in Figure 10), the slope deformations finally point to the outside of the slope and with the maximum magnitude smaller than 1 mm. This deformation behavior usually means the slope stability weakens.Figure 10
Distribution and evolution of deformations in slope.The groundwater pressure and deformation behavior in the slope rock mass during the operation period are dominated by hydraulic-mechanical coupling or referred to as the stress-seepage coupling.Owing to hysteresis of seepage, the excess groundwater pressure appears in the riverbed rock mass in phases A-B for the fluid in rock mass compressed by an increase of ponded water load and will experience dissipation during the decease of the reservoir water level in phases B-C and C-D, as represented in Figure9. Many studies have shown that the hysteresis characteristic of seepage mainly depends on the permeability of the rock and the variation speed of the reservoir water level. The higher permeability of the rock and the lower change speed of the water level cause short-term hysteretic duration. It needs to be mentioned that the distributions of groundwater pressure are almost identical at the time of D and A in Figure 9, which indicates a critical phenomenon; the seepage-stress coupling process would reach the steady-state in one operation period which may be due to a longer time with unchanged water level in phases C-D. In addition, the evolution process of the groundwater pressure in Figure 9 also shows that the depth affected by water level fluctuation varies up to 40 m in the left slope, 45 m in the right slope, and 30 m in the riverbed.


## 5.2. Decoupling Solution for Stress-Seepage Coupling Simulation

The characteristic times of stress propagation and fluid diffusion in equivalent continuum such as united porous media are defined as(7)Stresspropagation:tcm=ρdKu+4/3GLc,Fluiddiffusion:tcf=Lc2c,where Ku is the undrained bulk modulus of the medium.The comparison between the above-mentioned two characteristic times illustrates the coupling degree of stress and seepage, which would help choose the proper approach with high computational efficiency for stress-seepage coupling calculation, i.e., high coupling degree of seepage stress is represented with identical magnitudes of the characteristic times. In general,tcm is less than tcf, while tcm≫tcf and tcm≪tcf, respectively, describe two extreme permeability properties of rock or soil, i.e., completely drained and undrained. Objective to these two extreme situations, the decoupling approach can be introduced to solve the stress-seepage coupling problem, in which the seepage and stress field are solved separately. Although the decoupling approach simplifies the coupling process, it is widely used in reality for obtaining higher computational efficiency.The decoupling approach is also available in this study fortcf more than 500 times to tcm according to the rock properties of the Zhouning project. The decoupling approach is designed to be available to analyze stress-seepage coupling for every typical reservoir water varying phase, i.e., A-B, B-C, and C-D in Figure 5. In one phase of decoupling modeling, the computation needs to be decomposed into two main stages: (1) perform stress calculation with static mode, in which groundwater is not allowed to flow; (2) implement seepage simulation in the real-time domain to have the stress field disturbed by groundwater diffusion.A conceptual plane model as shown in Figure11 was built aimed at illustrating the decoupling approach and further investigating groundwater pressure and deformation behavior in the slope during the impoundment period. In the conceptual model, the rock mass is considered as an elastic material to ignore its nonlinearity for highlighting the stress-seepage coupling mechanism, and the permeability and mechanical parameters are consistent with the upper weakly weathered rock mass of Zhouning slope (Class III). In addition, an idealized steep topography was designed to deform the behavior at the crest of the slope.Figure 11
Deformation mechanism of the slope during impoundment.
(a)
Mattress effect before seepage
(b)
Swelling effect after seepageAccording to Figure11, the decoupling computation results indicate that there are two successive critical effects that exist in the phase of reservoir impoundment, referred to as “mattress effect” and “swelling effect” in this paper. The impoundment was modeled as a sudden water level increase from the dead level to the normal storage level, meanwhile, ignoring the groundwater diffusion for low permeability of the slope rock mass; consequently, the mattress effect occurs, in which compressive deformation occurs under the water level and tensile deformation above the water level. The excess groundwater pressure in the riverbed is also observed as shown in A of Figure 10. After that, seepage simulation is implemented in the real-time domain; consequently, the swelling effect is observed, in which the whole model is uplifted as a result of effective stress reduction due to groundwater diffusion in the slope rock mass. It should be noted that the “swelling effect” is closely related to the time for the groundwater to diffuse; the longer the groundwater diffusion time, the more obvious the “swelling effect” until a state of equilibrium.Objectively, the two mechanisms of “mattress effect” and “swelling effect” run through the whole process of slope stress-seepage coupling dynamically and synchronously until the disturbance reaches an equilibrium state. The slope deformation during the impoundment period is the result of the superposition of the two effects.Decoupling analysis reveals the mechanism of seepage-stress in the rock mass in a manner of simplification and is identical to the analysis of undrained consolidation settlement.The decrease of the water level is the inverse process of impoundment, in which the ponded water load is gradually removed, and the seepage force changes to point to the outside of the slope. At this time, besides the hydraulic behavior, the slope deformation and stability are closely related to the mechanical properties of the rock mass, which is further investigated in a further section.


## 6. Prediction of Slope Stability during the Operation Period

The slope stability relates closely to the hydraulic behavior, i.e., the direction of seepage and the hydraulic gradient. In the process of a sudden drop of reservoir water level, i.e., phases B-C in Figure5, the ponded water load on the slope surface is released rapidly (unloading), and the seepage force in the rock mass changes point to the outside of the slope, as shown in Figure 12. Moreover, because the groundwater adjustment lags behind the drop speed of the reservoir water level, the seepage force formed in the slope may be higher. Subject to these above hydraulic behaviors, the rock mass will bear tensile stress. Therefore, the tensile strength of rock mass σtmass should be emphasized when predicting the slope stability during the operation period with high frequent fluctuation of the reservoir water level.Figure 12
Direction of seepage changes during the operation period.The rock mass is usually characterized by shear failure in engineering practice; the shear behavior is widely concerned, but the tensile strength is not considered as much, under which the rock is easily broken. The tensile strength of rock is determined by some empirical methods instead of laboratory tests in this study.(1)
So-called truncation tensile strength is determined by the maximum value of principal stress (compressive stresses are negative) at the yield surface of Mohr-Coulomb(2)
Proportional tensile strength is determined as 1/5~1/20 of the uniaxial compressive strength obtained in the laboratory(3)
The tensile strength is determined by the Hoek-Brown (HB) strength criterionThe above-mentioned three methods are expressed as the following equations:(8)Method1:σtmass=ctanϕ,Method2:σtmass=λ2cNϕMethod3:σtmass=sσcimb,where Nϕ=1+sinϕ/1−sinϕ, c, and ϕ are the cohesion and internal friction angle of the rock mass; the proportional factor λ is here adopted as 1/10; mb and s are the coefficients related to the lithology and rock mass classification in HB strength criterion; and σci is the uniaxial compressive strength. Figure 13 illustrates that the tensile strength was estimated for various types of rock masses in the Zhouning project. It can be seen that the truncation tensile strength is the largest, which may overestimate slope stability, but the value obtained by the HB method is more conserved. Note that the tensile strength for analysis in Section 4.1 is the truncation tensile strength.Figure 13
Tensile strength of rock masses estimated by three different methods.According to the above three estimated rock mass strengths, the stress seepage coupling process of the slope during the whole operation period (50 years/yr) is simulated by the decoupling approach. Specifically, every operation cycle (1 day/d) is simulated as three sets of decoupling analysis, corresponding to phases A-B, B-C, and C-D.Figure14 illustrates the contours of total and transverse deformation at the end of the operation (50 years) by using three different tensile strengths. It is noted that transverse deformation is approximately perpendicular to the river, used to indicate the direction of deformation. It can be observed that the total deformations estimated under truncation and proportional tensile strengths are under 1 mm at the end of the operation period. However, the maximum total displacement of close to 1 m occurs while using the tensile strength determined by the HB strength criterion, and its transverse deformation points to the outside of the slope, which is unfavorable to stability. The calculation results also show that the influence of fault F3 on the deformation of the slope is not obvious due to the pressure from the dam body on the slope, and so are the dikes.Figure 14
Distribution characteristics of slope deformation at the end of operation period (unit: m).
(a)
Total deformation (unit: m)
(b)
Transverse deformation (unit: m)For further observation on slope deformation behavior and stability, the deformation is traced at the typical location A of the right bank during the whole operation period. For the sake of clear observation, early-stage deformation at location A is represented in Figure15 where the deformation oscillates acutely with frequent fluctuation of the water level during the operation. The deformations determined based on the truncation and proportional tensile strength are identical compressive deformations and close to 4 mm, and the trend of the deformation process converges in a short time, which usually means the slope can keep stable during the operation period. However, the deformation corresponding to the HB tensile strength has no obvious trend of convergence, which means the slope may be in a state of failure.Figure 15
Deformation process of measuring point A.
(a)
Total deformation
(b)
Transverse deformationTo further judge the slope stability condition, Figure16 represents that the deformation at typical location A is extracted at the time interval of 2 years. In the whole operating time, the deformations corresponding to the truncation and proportional tensile strength are few and the instability failure phenomenon is not observed. The deformation rate obtained by HB tensile strength increases suddenly exhibiting divergence behavior at about the 8th year, which indicates that the slope is in the state of instability.Figure 16
Peak deformation process of typical measuring points.It can be concluded that the tensile strengths affect the slope stability, and the computation based on the HB tensile strength does not satisfy stability requirements during the operation. Note that the stability of the left bank slope is similar to that of the right bank slope as discussed above; no further description is detailed.


## 7. Conclusion

Taking the Zhouning pumped storage power station, for example, the stress-seepage coupling process during the operation period was studied with an equivalent continuum model. According to the numerical analysis, the following conclusions are summarized:(1)
The permeability characteristics of rock mass in the Zhouning project relate closely to the field fracture distribution and properties and are stress-dependent and anisotropic. The approach composed of 3D digital photogrammetry technology and a geological database management system is proposed to determine the permeability parameters for large rock slopes(2)
A hybrid grid approach in FLAC3D is proposed to build a large-scale 3D stress-seepage coupling model aimed at improving the computational efficiency, meanwhile ensuring accuracy. A hybrid grid model is the combination of the structured mesh and numerical element Attach(3)
The deformation of the slope during the impoundment period is the result of the superposition of two mechanisms, referred to as the mattress effect and the swelling effect(4)
Slope rock mass will bear tensile stress in the process of a sudden drop of reservoir water level, mainly due to ponded water load releasing rapidly and a much higher seepage force pointing to out of the slope, which means the tensile strength of rock should be considered in predicting rock mass mechanical behavior during the operation period. If the tensile strength is determined by the HB strength criterion, the slope will be in a state of instability while operating and will not meet stability requirementsIt must be pointed out that the Fujian Zhouning pumped storage power station is still in the construction phase, and the numerically simulated results need to be reviewed and even further studied when the power station is in actual operation.



---

*Source: 9926751-2021-06-28.xml*