Heightening of an Existing Embankment Dam: Results from Numerical Simulations

Abstract

The old dam of the Zhushou Reservoir is a clay core rock-debris dam with a maximum height of 63.4 m. After heightening, the new dam is a concrete-faced rockfill dam with a maximum height of 98.1 m. In the initial design stage, a rigid connection is proposed between the cutoff wall and toe slab. After the concrete cutoff wall is built at the axis of the old dam, a complete cutoff system is composed of cutoff wall, toe slab, and face slab. In this paper, based on the static and dynamic tests of dam materials, the Shen Zhujiang double-yield surface elastic-plastic model is used as the static constitutive model, and the contact friction model is used as the contact surface model. The three-dimensional finite element method is used to simulate the construction filling and water storage process during operation. The simulation results show that the maximum horizontal displacement occurs in the dam body of the old dam and the maximum settlement occurs at the interface between the old and new dams. During the storage period, the cutoff wall will not be damaged, and the tensile stress of the local area at the junction of toe slab and bank slope has exceeded the allowable value for C30 plain concrete, so the reinforcement should be strengthened at this location.

1. Introduction

With the rapid development of the economy and the large-scale development of water energy, the construction of reservoir dams has become an important engineering initiative to meet the needs of social and economic development. Over time, the sediment in the reservoir continues to accumulate, and the storage capacity for prosperity and flood control continues to decrease or even loses capacity completely. In addition, due to a lack of water level data and drainage area data in the original design or a lack of labor, equipment, funds, or other resources during construction, the construction of small storage capacity reservoirs cannot meet the current demand for water resources. Therefore, the construction of new water conservancy facilities or the heightening of the old dams has become an urgent problem to consider. Compared with the construction of a new dam, raising an original dam body does not require the consideration of the location of a new dam, and it can obtain a larger storage capacity at a lower economic cost. Therefore, the dam heightening scheme has attached increasing attentions from engineers [1].

Addressing the technical problems that rise during the process of heightening is becoming a top priority due to the large amount of work and the complexity of construction technology. There are different key problems in dam heightening engineering due to the dam type and heightening method. Earth-rock dams are a widely used type of dam. Due to the permeability of earth-rock materials, it is urgent to study the impact of seepage on the earth-rock dam during the heightening process [2]. For slope-type heightened and thickened concrete dams, the key issues related to dam heightening are the stress concentrations and deformation of the dam body during construction and operation, stress analysis and structural form of the interface between new and old concrete, and design of drainage and water stop [3]. Periodic changes in the temperature and changes in the temperature of the old dam after new concrete is placed will cause problems such as deterioration of the dam heel stress, cracks in the joint surface, and cracks in the downstream dam surface [4].

There are many engineering precedents for dam elevation, such as the Goscheneralp Dam and Grande Dixence Dam in Switzerland, Steenbras Dam in South Africa, Roseires Dam in Sudan, and Danjiangkou Dam and Songyue Dam in China [3, 5, 6, 7, 8, 9]. Due to the rapid increases in the urban population of Cape Town, raising the Steenbras Dam offered an effective solution to the problem of a serious water shortage. During the course of anchoring the dam, engineers considered that post-stressing would have advantages in terms of cost and expedition. Essentially the process is one of placing vertical cables through the wall of a mass concrete dam from the crest into the foundation and stressing the cables to produce stabilizing compressive forces on the upstream face. Similar to the Steenbras Dam, the Songyue Dam also raised the dam to meet the water supply needs of Helong City. The Songyue Dam is located in a severely cold area, with an average annual temperature of 4.8°C, and the temperature changes greatly during the year. Therefore, the heightened structure needed to adapt to the characteristics of the severely cold area. The calculation research on the Songyue Dam heightening scheme shows that setting a sliding joint in the middle of the joint surface can absorb the shrinkage and deformation of a part of the newly poured concrete, which has a significant effect on improving the tensile stress of the upstream and downstream dam surfaces.

The Zhushou Reservoir is located in Sichuan Province, China, which is located in a seismically active area. The dam of Zhushou Reservoir is a clay core rock-debris dam. To meet the production and domestic water demand of nearby cities, it is necessary to expand the capacity of the Zhushou Reservoir. Under the action of gravity loads, water loads, and earthquake loads, effectively coordinating the deformation of the rockfill of the new and old dams to allow the stress and deformation of the seepage control system to be within the allowable range of the materials is a major technical difficulty to be solved. Therefore, based on the experience of previous engineering technologies, the necessary theoretical research is carried out to accurately predict the stress and deformation of the dam, especially the coordination between the deformation of the old and new dams, to improve the rationality of engineering design and to improve future engineering operations.

2. Project description

2.1 General situation of the Zhushou Reservoir project

The Zhushou Reservoir pivotal project is located in Liangshan Prefecture, Sichuan Province, and is a medium-sized reservoir. The dam is made of a clay core rock-debris dam. Its top elevation is 2416.10 m, the dam length is 190.00 m, the top elevation of the wave-proof wall is 2417.10 m, the dam height is 63.4 m, and the width of the dam top is 6.0 m. Both the upper and lower dams are provided with rockfilled prisms. The upstream slope is protected by dry block stone, while the downstream slope is protected by a dry block stone arch ring and turf in the circle. The thickness of the dry block stone is 40 cm. The top width of the gravel soil core wall is 6.0 m, the top elevation is 2415.3 m, and the upper and lower slopes are 1:0.4.

2.2 Dam heightening scheme

According to the water supply project planning of the Baihetan hydropower station resettlement area, to meet the production and domestic water demand of the resettlement area, the Zhushou Reservoir should be expanded and matched to the corresponding water diversion project. The dam should be increased from 63.4 m to 98.1 m. At the same time, when the dam is heightened, the impervious body of the original dam should be strengthened [10].

The objective of dam heightening is to make use of the water-retaining capacity of the original core wall dam to produce rockfill heightening on the top and downstream slope of the old dam so that the original dam body becomes a part of the heightened dam. At the same time, a core wall and foundation anti-seepage system of the original dam is strengthened, a concrete cutoff wall is added, and the foundation anti-seepage curtain grouting is strengthened. The anti-seepage type of the heightening dam body adopts the upstream reinforced concrete-faced slab, the slope ratio of the upstream dam is 1:1.4, and the comprehensive slope of the downstream rockfill body is 1:1.6 [11]. Figures 1 and 2 show general view of the Zhushou Reservoir dam.

Figure 1.

Plane figure of heightening of the Zhushou Reservoir dam.

Figure 2.

Standard profile of heightening of the Zhushou Reservoir dam.

2.3 Dam heightening construction procedure

To avoid excessive deformation and cracking of the lower core wall caused by the compression of the upper high rockfill, the cutoff wall is constructed after the upper rockfill body is filled and settled for 3 months. The concrete connecting plate between the cutoff wall and the toe slab shall be constructed after the toe slab and the face plate are completed.

The overall construction procedure is as follows: old dam filling → new dam filling to 2447.90 m → core wall reinforcement and cutoff wall construction → toe slab construction → panel construction → connecting plate construction → new dam filling to 2451 m. The water level remains at 2395.0 m during the construction period. The construction period of dam heightening is 31 months, which are as follows:

From September of the first year to February of the second year, the construction period of the old dam filling is 6 months.

From March of the second year to November of the second year, the construction period of the new dam filling to an elevation of 2447.9 m is 9 months.

From December of the second year to May of the third year, the construction period of core wall reinforcement and cutoff wall construction is 6 months.

During June of the third year, the construction period of toe slab is 1 month.

From July of the third year to August of the third year, the construction period of concrete panel and wave wall construction is 2 months.

From September of the third year to November of the third year, the construction period of connecting plate construction is 3 months.

From December of the third year to July of the fourth year, the construction period of new dam filling to 2451 m is 4 months.

3. Finite element modeling

3.1 Finite element meshing

Figure 3 shows a finite element mesh diagram of a typical riverbed section, Figure 4 shows a three-dimensional finite element mesh diagram, and Figure 5 shows an anti-seepage system (core wall, connecting plate, toe slab and panel) meshing diagram, where the X forward direction is defined as from the left bank to the right bank, the Y forward direction is defined as upstream to downstream, and the Z forward direction is defined as the opposite direction of gravity. The three-dimensional solid element adopts an 8-node hexahedral isoparametric element and a 4-node tetrahedral isoparametric element, and the latter is treated as a degenerated hexahedral element. There are 29,905 generating units and 33,482 nodes in total.

Figure 3.

Finite element mesh diagram of a typical riverbed section.

Figure 4.

Three-dimensional finite element mesh diagram.

Figure 5.

Anti-seepage system meshing diagram.

3.2 Construction sequence and simulation of the water storage process

According to the above construction and water storage process, the order of the filling and storage simulation in the finite element calculation is as follows: old dam filling → new dam filling to 2447.90 m (the water level remained at 2395.0 m) → cutoff wall construction → toe slab construction → panel construction → connecting plate construction → new dam filling to 2451 m → upstream water storage to a normal water level elevation of 2444 m. There are 70 stages for simulation, including 42 stages for dam filling and 38 stages for water storage. Figure 6 shows the simulation diagram of the Zhushou Reservoir construction and water storage process. Figure 7 shows the water level-time curve of the Zhushou Reservoir during the construction and water storage process.

Figure 6.

The Zhushou reservoir construction and water storage process simulation diagram.

Figure 7.

The Zhushou reservoir water level-time curve during the construction and water storage process.

4. Constitutive model of the dam material

4.1 Constitutive model of the soil and rockfill

As the main body of the concrete-faced rockfill dam, reasonable simulation of its stress–strain relationship is very important to improve the rationality of the calculation results of the stress and deformation of the concrete-faced rockfill dam. In this project, the constitutive model of rockfill material is based on the Shen Zhujiang double-yield surface elastic-plastic model proposed by Shen Zhujiang. Compared with the nonlinear elastic model, the model can consider the dilatancy and shear-shrinkage characteristics of rockfill bodies and can more accurately reflect the stress-strain characteristics of dam bodies than other models.

In the Shen Zhujiang double-yield surface elastic-plastic model, the two-yield surfaces are only regarded as the boundary of elastic region and are no longer related to hardening parameters. The double-yield surface is used to establish the unloading criterion, make the elastic-plastic matrix symmetrical, and specify the direction of plastic strain. As shown in Figure 8, due to the double-yield surface, not only the loading direction B will produce plastic strain, but also the loading directions A and C will produce plastic strain.

Figure 8.

Double hardening model.

The two-yield surfaces of the Shen Zhujiang double-yield surface elastic-plastic model are

$$\left.\begin{array}{l}{\mathit{F}}_1={\mathit{p}}^2+{\mathit{r}}^2{\mathit{q}}^2-{\mathit{f}}_1=0\\ {\mathit{F}}_2={\mathit{q}}^{\mathit{s}}/\mathit{p}-{\mathit{f}}_2=0\end{array}\right\}$$E1

where $\mathit{p}=\frac{1}{3}\left({\sigma }_1+{\sigma }_2+{\sigma }_3\right)$, $\mathit{q}=\frac{1}{\sqrt{2}}{\left[{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2\right]}^{1/2}$, and $\mathit{r}$ and $\mathit{s}$ are model parameters and can be taken as 2 for rockfill materials. The expression of the strain increment of the double-yield surface model is as follows:

$$\left\{\mathrm{\Delta }\epsilon \right\}={\left[\mathit{D}\right]}^{-1}\left\{\mathrm{\Delta }\sigma \right\}+{\mathit{A}}_1\left\{{\mathit{n}}_1\right\}{\left\{\frac{\partial {\mathit{f}}_1}{\partial \sigma }\right\}}^{\mathit{T}}\left\{\mathrm{\Delta }\sigma \right\}+{\mathit{A}}_2\left\{{\mathit{n}}_2\right\}{\left\{\frac{\partial {\mathit{f}}_2}{\partial \sigma }\right\}}^{\mathit{T}}\left\{\mathrm{\Delta }\sigma \right\}$$E2

where $\left[\mathit{D}\right]$ is the elastic matrix, $\left\{{\mathit{n}}_1\right\}$ and $\left\{{\mathit{n}}_2\right\}$ are the normal direction cosines of the yield surface, and ${\mathit{A}}_1$ and ${\mathit{A}}_2$ are plasticity coefficients. $\Delta {\mathit{f}}_1$ and $\Delta {\mathit{f}}_2$ can be written as.

$$\left.\begin{array}{l}\Delta {\mathit{f}}_1=2\mathit{p}\Delta \mathit{p}+2{\mathit{r}}^2\mathit{q}\Delta \mathit{q}\\ \Delta {\mathit{f}}_2={\mathit{q}}^{\mathit{s}}/\mathit{p}\left(-\Delta \mathit{p}/\mathit{p}+\mathit{s}\Delta \mathit{q}/\mathit{q}\right)\end{array}\right\}$$E3

The model adopts the normal flow rule, so the plastic potential surface is orthogonal to the direction of the plastic strain increase and ${\mathit{Q}}_1={\mathit{F}}_1$ and ${\mathit{Q}}_2={\mathit{F}}_2$. According to the normal flow rule, $\partial {\mathit{Q}}_1/\partial {\sigma }_1$, $\partial {\mathit{Q}}_1/\partial {\sigma }_3$, $\partial {\mathit{Q}}_3/\partial {\sigma }_1$, and $\partial {\mathit{Q}}_3/\partial {\sigma }_3$ can be calculated. Considering that $\mathit{P}=\frac{1}{3}\left({\sigma }_1+2{\sigma }_3\right)$, $\mathit{q}={\sigma }_1-{\sigma }_3$, and $\Delta {\epsilon }_{\mathit{v}}=\Delta {\epsilon }_1+2\Delta {\epsilon }_3$ under triaxial conditions, Eq. (2) can be expressed as:

$$\left.\begin{array}{l}\frac{\Delta {\epsilon }_1}{\Delta {\sigma }_1}=\frac{1}{\mathit{E}}+\frac{4}{9}{\left(\mathit{p}+3{\mathit{r}}^2\mathit{q}\right)}^2{\mathit{A}}_1+\frac{1}{9}{\left(\frac{1}{\mathit{p}}-\frac{3\mathit{s}}{\mathit{q}}\right)}^2\frac{{\mathit{q}}^{2\mathit{s}}}{{\mathit{p}}^2}{\mathit{A}}_2\\ \frac{\Delta {\epsilon }_{\mathit{v}}}{\Delta {\sigma }_1}=\frac{1-2\nu }{\mathit{E}}+\frac{4}{3}\mathit{p}\left(\mathit{p}+3{\mathit{r}}^2\mathit{q}\right){\mathit{A}}_1+\frac{1}{3}\left(\frac{1}{\mathit{p}}-\frac{3\mathit{s}}{\mathit{q}}\right)\frac{{\mathit{q}}^{2\mathit{s}}}{{\mathit{p}}^2}{\mathit{A}}_2\end{array}\right\}$$E4

Tangent Young’s modulus is defined as ${\mathit{E}}_{\mathit{t}}=\frac{\Delta {\sigma }_1}{\Delta {\epsilon }_1}$, and tangent volume ratio is defined as ${\mu }_{\mathit{t}}=\frac{\Delta {\epsilon }_{\mathit{v}}}{\Delta {\epsilon }_1}$. According to ${\mathit{E}}_{\mathit{t}}$ and ${\mu }_{\mathit{t}}$, ${\mathit{A}}_1$ and ${\mathit{A}}_2$ can be expressed as

$$\left.\begin{array}{l}{\mathit{A}}_1=\frac{1}{4{\mathit{p}}^2}\frac{\eta \left(\frac{9}{{\mathit{E}}_{\mathit{t}}}-\frac{3{\mu }_{\mathit{t}}}{{\mathit{E}}_{\mathit{t}}}-\frac{3}{{\mathit{G}}_{\mathit{e}}}\right)+2\mathit{s}\left(\frac{3{\mu }_{\mathit{t}}}{{\mathit{E}}_{\mathit{t}}}-\frac{1}{{\mathit{B}}_{\mathit{e}}}\right)}{2\left(1+3{\mathit{r}}^2\eta \right)\left(\mathit{s}+{\mathit{r}}^2{\eta }^2\right)}\\ {\mathit{A}}_2=\frac{{\mathit{p}}^2{\mathit{q}}^2}{{\mathit{q}}^{2\mathit{s}}}\frac{\left(\frac{9}{{\mathit{E}}_{\mathit{t}}}-\frac{3{\mu }_{\mathit{t}}}{{\mathit{E}}_{\mathit{t}}}-\frac{3}{{\mathit{G}}_{\mathit{e}}}\right)-2{\mathit{r}}^2\eta \left(\frac{3{\mu }_{\mathit{t}}}{{\mathit{E}}_{\mathit{t}}}-\frac{1}{{\mathit{B}}_{\mathit{e}}}\right)}{2\left(3\mathit{s}-\eta \right)\left(\mathit{s}+{\mathit{r}}^2{\eta }^2\right)}\end{array}\right\}$$E5

where $\eta =\mathit{q}/\mathit{p}$ and ${\mathit{G}}_{\mathit{e}}$ and ${\mathit{B}}_{\mathit{e}}$ are the elastic shear modulus and bulk modulus, respectively:

$${\mathit{G}}_{\mathit{e}}={\mathit{E}}_{\mathit{ur}}/2\left(1+\nu \right)$$E6

$${\mathit{B}}_{\mathit{e}}={\mathit{E}}_{\mathit{ur}}/3\left(1-2\nu \right)$$E7

In the formula, the elastic Poisson’s ratio $\nu$ is 0.3, and ${\mathit{E}}_{\mathit{ur}}$ is the modulus of the unloading resilience. The tangent Young’s modulus ${\mathit{E}}_{\mathit{t}}$ and tangent volume ratio ${\mu }_{\mathit{t}}$ in Eq. (5) are two basic variables of the model, which are expressed as follows:

$${\mathit{E}}_{\mathit{t}}={\mathit{KP}}_{\mathit{a}}{\left(\frac{{\sigma }_3}{{\mathit{P}}_{\mathit{a}}}\right)}^{\mathit{n}}{\left(1-{\mathit{R}}_{\mathit{f}}{\mathit{S}}_{\mathit{l}}\right)}^2$$E8

$${\mu }_{\mathit{t}}=2{\mathit{c}}_{\mathit{d}}{\left(\frac{{\sigma }_3}{{\mathit{P}}_{\mathit{a}}}\right)}^{{\mathit{n}}_{\mathit{d}}}\frac{{\mathit{E}}_{\mathit{i}}{\mathit{R}}_{\mathit{s}}}{{\sigma }_1-{\sigma }_3}\frac{1-{\mathit{R}}_{\mathit{d}}}{{\mathit{R}}_{\mathit{d}}}\left(1-\frac{{\mathit{R}}_{\mathit{s}}}{1-{\mathit{R}}_{\mathit{s}}}\cdot \frac{1-{\mathit{R}}_{\mathit{d}}}{{\mathit{R}}_{\mathit{d}}}\right)$$E9

$${\mathit{S}}_{\mathit{l}}=\frac{\left({\sigma }_1-{\sigma }_3\right)\left(1-sin\varphi \right)}{2\mathit{c}cos\varphi +2{\sigma }_{3}sin\varphi }$$E10

where ${\mathit{P}}_{\mathit{a}}$ is the atmospheric pressure, $\mathit{K}$ is the Young’s modulus coefficient, and $\mathit{n}$ is the power of the tangent Young’s modulus ${\mathit{E}}_{\mathit{t}}$, which increases with the increase in the confining pressure ${\sigma }_3$. ${\mathit{R}}_{\mathit{f}}$ is the failure ratio, ${\mathit{S}}_{\mathit{l}}$ is the stress level, and $\mathit{c}$ and $\varphi$ are the shear strength indexes; ${\mathit{R}}_{\mathit{s}}={\mathit{R}}_{\mathit{f}}{\mathit{S}}_{\mathit{l}}$. ${\mathit{R}}_{\mathit{d}}$, ${\mathit{C}}_{\mathit{d}}$, and ${\mathit{n}}_{\mathit{d}}$ are calculation parameters; ${\mathit{C}}_{\mathit{d}}$ corresponds to the maximum shrinkage volume strain when ${\sigma }_3$ equals the unit atmospheric pressure; ${\mathit{n}}_{\mathit{d}}$ is the power of the shrinkage volume strain which increases with the increase in ${\sigma }_3$, and ${\mathit{R}}_{\mathit{d}}$ is the ratio of ${\left({\sigma }_1-{\sigma }_3\right)}_{\mathit{d}}$ to the asymptotic value of the deviating stress ${\left({\sigma }_1-{\sigma }_3\right)}_{\mathit{ult}}$ when the maximum shrinkage occurs. The elastoplastic matrix of the double-yield surface model can be obtained from the inverse of Eq. (2):

$$\left\{\mathit{\Delta \sigma }\right\}={\left[\mathit{D}\right]}_{\mathit{ep}}\left\{\mathit{\Delta \epsilon }\right\}$$E11

However, the expression of ${\left[\mathit{D}\right]}_{\mathit{ep}}$ is quite complex. In order to simplify the ${\left[\mathit{D}\right]}_{\mathit{ep}}$, Prandtl-Reuss flow rule is adopted in $\pi$ plane; ${\left[\mathit{D}\right]}_{\mathit{ep}}$ can be expressed as

$${\left[\mathit{D}\right]}_{\mathit{ep}}=\left[\begin{array}{cccccc}{\mathit{d}}_{11}& & & & & \\ {\mathit{d}}_{21}& {\mathit{d}}_{22}& & & \text{symmetric}& \\ {\mathit{d}}_{31}& {\mathit{d}}_{32}& {\mathit{d}}_{33}& & & \\ {\mathit{d}}_{41}& {\mathit{d}}_{42}& {\mathit{d}}_{43}& {\mathit{d}}_{44}& & \\ {\mathit{d}}_{51}& {\mathit{d}}_{52}& {\mathit{d}}_{53}& {\mathit{d}}_{54}& {\mathit{d}}_{55}& \\ {\mathit{d}}_{61}& {\mathit{d}}_{62}& {\mathit{d}}_{63}& {\mathit{d}}_{64}& {\mathit{d}}_{65}& {\mathit{d}}_{66}\end{array}\right]$$E12

where ${\mathit{d}}_{11}={\mathit{M}}_1-\mathit{P}\left({\mathit{S}}_{\mathit{x}}+{\mathit{S}}_{\mathit{x}}\right)/\mathit{q}-{\mathit{QS}}_{\mathit{x}}^2/{\mathit{q}}^2$, ${\mathit{d}}_{22}={\mathit{M}}_1-\mathit{P}\left({\mathit{S}}_{\mathit{y}}+{\mathit{S}}_{\mathit{y}}\right)/\mathit{q}-{\mathit{QS}}_{\mathit{y}}^2/{\mathit{q}}^2$, ${\mathit{d}}_{33}={\mathit{M}}_1-\mathit{P}\left({\mathit{S}}_{\mathit{z}}+{\mathit{S}}_{\mathit{z}}\right)/\mathit{q}-{\mathit{QS}}_{\mathit{z}}^2/{\mathit{q}}^2$, ${\mathit{d}}_{44}={\mathit{G}}_{\mathit{e}}-\mathit{Q}{\tau }_{\mathit{xy}}^2/{\mathit{q}}^2$, ${\mathit{d}}_{55}={\mathit{G}}_{\mathit{e}}-\mathit{Q}{\tau }_{\mathit{yz}}^2/{\mathit{q}}^2$, ${\mathit{d}}_{66}={\mathit{G}}_{\mathit{e}}-\mathit{Q}{\tau }_{\mathit{zx}}^2/{\mathit{q}}^2$, ${\mathit{d}}_{12}={\mathit{M}}_2-\mathit{P}\left({\mathit{S}}_{\mathit{x}}+{\mathit{S}}_{\mathit{y}}\right)/\mathit{q}-{\mathit{QS}}_{\mathit{x}}{\mathit{S}}_{\mathit{y}}/{\mathit{q}}^2$, ${\mathit{d}}_{13}={\mathit{M}}_2-\mathit{P}\left({\mathit{S}}_{\mathit{x}}+{\mathit{S}}_{\mathit{z}}\right)/\mathit{q}-{\mathit{QS}}_{\mathit{x}}{\mathit{S}}_{\mathit{z}}/{\mathit{q}}^2$, ${\mathit{d}}_{14}=-\mathit{P}{\tau }_{\mathit{xy}}/\mathit{q}-{\mathit{QS}}_{\mathit{x}}{\tau }_{\mathit{xy}}/{\mathit{q}}^2$, ${\mathit{d}}_{15}=-\mathit{P}{\tau }_{\mathit{yz}}/\mathit{q}-{\mathit{QS}}_{\mathit{x}}{\tau }_{\mathit{yz}}/{\mathit{q}}^2$, ${\mathit{d}}_{16}=-\mathit{P}{\tau }_{\mathit{zx}}/\mathit{q}-{\mathit{QS}}_{\mathit{x}}{\tau }_{\mathit{zx}}/{\mathit{q}}^2$, ${\mathit{d}}_{23}={\mathit{M}}_2-\mathit{P}\left({\mathit{S}}_{\mathit{y}}+{\mathit{S}}_{\mathit{z}}\right)/\mathit{q}-{\mathit{QS}}_{\mathit{y}}{\mathit{S}}_{\mathit{z}}/{\mathit{q}}^2$, ${\mathit{d}}_{24}=-\mathit{P}{\tau }_{\mathit{xy}}/\mathit{q}-{\mathit{QS}}_{\mathit{y}}{\tau }_{\mathit{xy}}/{\mathit{q}}^2$, ${\mathit{d}}_{25}=-\mathit{P}{\tau }_{\mathit{yz}}/\mathit{q}-{\mathit{QS}}_{\mathit{y}}{\tau }_{\mathit{yz}}/{\mathit{q}}^2$, ${\mathit{d}}_{26}=-\mathit{P}{\tau }_{\mathit{zx}}/\mathit{q}-{\mathit{QS}}_{\mathit{y}}{\tau }_{\mathit{zx}}/{\mathit{q}}^2$, ${\mathit{d}}_{34}=-\mathit{P}{\tau }_{\mathit{xy}}/\mathit{q}-{\mathit{QS}}_{\mathit{z}}{\tau }_{\mathit{xy}}/{\mathit{q}}^2$, ${\mathit{d}}_{35}=-\mathit{P}{\tau }_{\mathit{yz}}/\mathit{q}-{\mathit{QS}}_{\mathit{z}}{\tau }_{\mathit{yz}}/{\mathit{q}}^2$, ${\mathit{d}}_{36}=-\mathit{P}{\tau }_{\mathit{zx}}/\mathit{q}-{\mathit{QS}}_{\mathit{z}}{\tau }_{\mathit{zx}}/{\mathit{q}}^2$, ${\mathit{d}}_{45}=-\mathit{P}{\tau }_{\mathit{yz}}/\mathit{q}-\mathit{Q}{\tau }_{\mathit{xy}}{\tau }_{\mathit{yz}}/{\mathit{q}}^2$, ${\mathit{d}}_{46}=-\mathit{P}{\tau }_{\mathit{zx}}/\mathit{q}-\mathit{Q}{\tau }_{\mathit{xy}}{\tau }_{\mathit{zx}}/{\mathit{q}}^2$, and ${\mathit{d}}_{56}=-\mathit{P}{\tau }_{\mathit{zx}}/\mathit{q}-\mathit{Q}{\tau }_{\mathit{yz}}{\tau }_{\mathit{zx}}/{\mathit{q}}^2$. In the formula, ${\mathit{S}}_{\mathit{x}}={\sigma }_{\mathit{x}}-\mathit{p}$, ${\mathit{S}}_{\mathit{y}}={\sigma }_{\mathit{y}}-\mathit{p}$, ${\mathit{S}}_{\mathit{z}}={\sigma }_{\mathit{z}}-\mathit{p}$, ${\mathit{M}}_1={\mathit{K}}_{\mathit{p}}+4{\mathit{G}}_{\mathit{e}}/3$, ${\mathit{M}}_2={\mathit{K}}_{\mathit{p}}-2{\mathit{G}}_{\mathit{e}}/3$, $\mathit{P}={\mathit{B}}_{\mathit{e}}{\mathit{G}}_{\mathit{e}}\gamma /\left(1+{\mathit{B}}_{\mathit{e}}\alpha +{\mathit{G}}_{\mathit{e}}\delta \right)$, and $\mathit{Q}={\mathit{G}}_{\mathit{e}}^2\delta /\left(1+{\mathit{B}}_{\mathit{e}}\alpha +{\mathit{G}}_{\mathit{e}}\delta \right)$ are as follows:

$${\mathit{K}}_{\mathit{p}}=\frac{{\mathit{B}}_{\mathit{e}}}{1+{\mathit{B}}_{\mathit{e}}\alpha }\left(1+\frac{{\mathit{B}}_{\mathit{e}}{\mathit{G}}_{\mathit{e}}{\gamma }^2}{1+{\mathit{B}}_{\mathit{e}}\alpha +{\mathit{G}}_{\mathit{e}}\delta }\right)$$E13

$$\alpha =\frac{{\mathit{A}}_1}{{\mathit{r}}^2}+{\eta }^2{\mathit{A}}_2$$E14

$$\beta ={\mathit{r}}^2{\eta }^2{\mathit{A}}_1+{\mathit{s}}^2{\mathit{A}}_2$$E15

$$\gamma =\eta \left({\mathit{A}}_1-{\mathit{sA}}_2\right)$$E16

$$\delta =\beta +{\mathit{B}}_{\mathit{e}}\left(\mathit{\alpha \beta }-{\gamma }^2\right)$$E17

The Shen Zhujiang elastic-plastic model has eight model parameters, which are $\mathit{K}$, $\mathit{n}$, ${\mathit{R}}_{\mathit{f}}$, $\mathit{c}$, $\varphi$, ${\mathit{R}}_{\mathit{d}}$, ${\mathit{C}}_{\mathit{d}}$, and ${\mathit{n}}_{\mathit{d}}$. The parameters can also be determined from the results of conventional triaxial tests. Compared with the parameters of the Duncan $\mathit{E}-\nu$ model, only the latter three parameters ${\mathit{R}}_{\mathit{d}}$, ${\mathit{C}}_{\mathit{d}}$, and ${\mathit{n}}_{\mathit{d}}$ of the Shen Zhujiang elastic-plastic model are different from those of the Duncan model.

The Shen Zhujiang elastic-plastic model can also be calculated by the parameters of the model Duncan $\mathit{E}-\nu$ model. Its tangent volume ratio ${\mu }_{\mathit{t}}$ can be obtained from the tangent Poisson’s ratio ${\nu }_{\mathit{t}}$:

$${\mu }_{\mathit{t}}=1-2{\nu }_{\mathit{t}}$$E18

The tangent Poisson’s ratio ${\nu }_{\mathit{t}}$ in the Duncan $\mathit{E}-\nu$ model is as follows:

$${\nu }_{\mathit{t}}=\frac{{\mathit{G}}_{\mathit{e}}-\mathit{F}\cdot lg\left({\sigma }_3/\mathit{Pa}\right)}{{\left[1-\mathit{D}\left({\sigma }_1-{\sigma }_3\right)/{\mathit{E}}_{\mathit{i}}\left(1-{\mathit{R}}_{\mathit{f}}{\mathit{S}}_{\mathit{l}}\right)\right]}^2}$$E19

For unloading, the modulus of resilience is calculated as follows:

$${\mathit{E}}_{\mathit{ur}}={\mathit{K}}_{\mathit{ur}}{\mathit{P}}_{\mathit{a}}{\left(\frac{{\sigma }_3}{{\mathit{P}}_{\mathit{a}}}\right)}^{\mathit{n}}$$E20

where ${\mathit{K}}_{\mathit{ur}}$ is the modulus of resilience.

The loading and unloading criteria of the Shen Zhujiang elastic-plastic model are as follows:

1. If ${\mathit{F}}_1>{\mathit{F}}_{1max}$ and ${\mathit{F}}_2>{\mathit{F}}_{2max}$, it is fully loaded. $\mathit{A}{}_1>0$ and $\mathit{A}{}_2>0$.

2. If ${\mathit{F}}_1>{\mathit{F}}_{1max}$ and ${\mathit{F}}_2\le {\mathit{F}}_{2max}$ and ${\mathit{F}}_1\le {\mathit{F}}_{1max}$and ${\mathit{F}}_2>{\mathit{F}}_{2max}$, it is partially unloaded. $\mathit{A}{}_1=0$ or $\mathit{A}{}_2=0$.

3. If $F_1\le F_{1max}$ and $F_2\le F_{2max}$, it is fully unloaded. $\mathit{A}{}_1={\mathit{A}}_2=0$.

For coarse-grained materials, $\mathit{c}=0$ and $\varphi$ is calculated by the following formula:

$$\varphi ={\varphi }_{\mathit{o}}-\mathrm{\Delta }\varphi lg\left(\frac{{\sigma }_3}{{\mathit{P}}_{\mathit{a}}}\right)$$E21

where ${\varphi }_0$ and $\mathrm{\Delta }\varphi$ are material parameters determined by triaxial test results.

4.2 Constitutive model of concrete

The linear elastic model is used for concrete materials, and the stress–strain relationship conforms to the following generalized Hooke’s law:

$$\left\{\sigma \right\}=\left[\mathit{D}\right]\left\{\epsilon \right\}$$E22

where $\left[\mathit{D}\right]$ is an elastic matrix.

4.3 Interface model

At present, the Goodman thickness-free elements and Desai thin-layer elements are commonly used. Because the interface is a kind of interface without a thickness, it is more suitable to use the Goodman element without a thickness to theoretically simulate the interface. However, in fact, a Goodman element without a thickness must obtain a large normal stiffness to avoid overlap. In addition, shear dislocation does not necessarily occur on the interface and may penetrate into the soil at a certain distance. Desai thin-layer elements reflect normal deformation to a certain extent, but the choice of the thickness of thin-layer elements has a great influence on the calculation results. A large element thickness will introduce errors in the physics, and a small element thickness will introduce errors in the mathematics. Desai suggests that the ratio of the thickness $\mathit{t}$ to the length $\mathit{B}$ is as follows:

$$\frac{\mathit{t}}{\mathit{B}}=\frac{1}{10}\sim \frac{1}{100}$$E23

For the constitutive model of the contact surface, the hyperbolic model and ideal elastic-plastic model of the relationship between the shear stress and relative displacement proposed by Clough and Duncan are most commonly used. The results show that the shear stress on the interface between the soil and structure is not uniform, the shear deformation is actually a rigid-plastic deformation, and the contact friction model can be well simulated.

Before the shear stress $\tau$ on the contact surface reaches the destructive shear stress $\tau {}_{\mathit{f}}$, the dislocation deformation is very small, and the displacement of the contact surface is mainly caused by shear deformation. When the shear stress $\tau$ reaches the destructive shear stress $\tau {}_{\mathit{f}}$, the displacement of the contact surface is mainly caused by the dislocation deformation and can develop indefinitely.

The deformation on the contact surface can be divided into two parts: basic deformation and failure deformation. The basic deformation is similar to the deformation calculation model of other soils, expressed as $\left\{\mathrm{\Delta }{\epsilon }^{\prime }\right\}$. The failure deformation includes sliding failure and tension cracking failure, which only exists when the shear stress of the element reaches the shear strength and sliding failure along the contact surface or tension cracking failure occurs on the contact surface. In $\left\{\mathrm{\Delta }{\epsilon }^{″}\right\}$, the total deformation of the contact surface is as follows:

$$\left\{\mathrm{\Delta }\epsilon \right\}=\left\{\mathrm{\Delta }{\epsilon }^{\prime }\right\}+\left\{\mathrm{\Delta }{\epsilon }^{″}\right\}=\left[{\mathit{C}}^{\prime }\right]\left\{\mathrm{\Delta }\sigma \right\}+\left[{\mathit{C}}^{″}\right]\left\{\mathrm{\Delta }\sigma \right\}$$E24

There are two forms of failure and deformation of the elements: tension cracking and slip. The rigid-plastic model is used to calculate the relative shear deformation of the element. There is no relative slip on the contact surface before failure, and after failure, the relative slip will continue to develop.

For the three-dimensional thin-layer contact surface element, the Y direction is the normal direction of the contact surface:

$$\left\{\begin{array}{c}\mathrm{\Delta }{\epsilon }_{\mathit{x}}\\ \mathrm{\Delta }{\epsilon }_{\mathit{y}}\\ \mathrm{\Delta }{\epsilon }_{\mathit{z}}\\ \mathrm{\Delta }{\gamma }_{\mathit{xy}}\\ \mathrm{\Delta }{\gamma }_{\mathit{yz}}\\ \mathrm{\Delta }{\gamma }_{\mathit{zx}}\end{array}\right\}=\left[\begin{array}{cccccc}{\mathit{C}}_{11}& {\mathit{C}}_{12}& {\mathit{C}}_{13}& 0& 0& 0\\ {\mathit{C}}_{21}& {\mathit{C}}_{22}+\frac{1}{{\mathit{E}}^{\prime }}& {\mathit{C}}_{23}& 0& 0& 0\\ {\mathit{C}}_{31}& {\mathit{C}}_{32}& {\mathit{C}}_{33}& 0& 0& 0\\ 0& 0& 0& {\mathit{C}}_{44}+\frac{1}{{{\mathit{G}}_{\mathit{xy}}}^{\prime }}& 0& 0\\ 0& 0& 0& 0& {\mathit{C}}_{55}+\frac{1}{{{\mathit{G}}_{\mathit{yz}}}^{\prime }}& 0\\ 0& 0& 0& 0& 0& {\mathit{C}}_{66}\end{array}\right]\left\{\begin{array}{c}\mathrm{\Delta }{\sigma }_{\mathit{x}}\\ \mathrm{\Delta }{\sigma }_{\mathit{y}}\\ \mathrm{\Delta }{\sigma }_{\mathit{z}}\\ \mathrm{\Delta }{\tau }_{\mathit{xy}}\\ \mathrm{\Delta }{\tau }_{\mathit{yz}}\\ \mathrm{\Delta }{\tau }_{\mathit{zx}}\end{array}\right\}$$E25

If the contact surface is under tension, ${\mathit{E}}^{\prime }$ can be set to a very small value, such as ${\mathit{E}}^{\prime }$ = 5.0 kPa. If the contact surface is under compression, a larger value should be taken, or $\frac{1}{{\mathit{E}}^{\prime }}$ = 0.0.

The value of ${\mathit{G}}^{\prime }$ is determined by the stress level of the contact surface element: when the stress level is S < 0.99, the value of ${\mathit{G}}^{\prime }$ is larger, or $\frac{1}{{\mathit{G}}^{\prime }}$ = 0.0; when the stress level is S > 0.99, the shear failure of the contact surface element occurs, and the shear modulus corresponding to the residual stiffness of the element, or ${\mathit{G}}^{\prime }$ = 5.0 kPa, is obtained.

The flexibility matrix $\left[\mathit{C}\right]$ of the contact surface is directional. After forming the stiffness matrix ${\left[\mathit{C}\right]}^{‐1}$ in the local coordinate system, the element stiffness matrix in the global coordinate system needs to be obtained by a coordinate transformation.

4.4 Material parameters

4.4.1 Concrete

The dam uses C25 concrete and C30 concrete. The unit weight, elastic modulus and Poisson’s ratio of C25 and C30 concrete are 2.5 t/m³, 2.8 × 10⁴ N/mm², and 0.167 and 2.5 t/m³, 3.0 × 10⁴ N/mm², and 0.167, respectively.

4.4.2 Dam material

The lithology of the newly filled rockfill material in the Zhushou Reservoir is the Ordovician Hongshiya Formation (O1h) quartz sandstone, fine sandstone with silty mudstone, and Qiaojia Formation (O2q) gray thin-to-medium thick sandstone, dolomite, and limestone. According to design filling standards and field testing data, the triaxial CD test had been carried out on rockfill materials of the heightening dam and the filling materials of the old dam body. The calculated parameters determined by the test are shown in Table 1.

Material name		ρd (g/cm^3)	c (kPa)	φо (°)	Δφ (°)	k	n	Rf	D	F	G
Old dam	Gravel clay core wall	1.84	56.9	29.3	0	164.4	0.46	0.69	3.68	0.1	0.3
	Stone slag in the upper part of the dam hell (elevation above 2390 m)	2.04	0	41.8	9.1	318.8	0.46	0.79	2.78	0.04	0.35
	Stone slag in the lower part of the dam shell (elevation above 2390 m)	2.12	0	44.9	9.1	431.3	0.38	0.72	3.3	0.09	0.37
	Rockfill	2.14	0	47.7	10.1	811	0.31	0.54	10.4	0.12	0.4
New dam	Cushion zone	2.2	0	58.8	10.9	1245.6	0.35	0.60	10.4	0.12	0.4
	Transition region	2.17	0	59.5	13.3	1405.4	0.29	0.65	10.1	0.15	0.39
	Main rockfill area	2.16	0	59.4	13.6	1301.5	0.27	0.60	9.3	0.15	0.39
	Secondary rockfill area	2.14	0	56.1	11.9	954.1	0.37	0.63	9.5	0.13	0.36

Table 1.

Calculated parameters of the dam material.

5. Results and discussion

5.1 Stress and deformation of the dam body

Considering the stress and deformation of the new dam after filling and storage period and influence of the stress and deformation of the new dam on the old dam, Table 2 lists the characteristic values of the stress and deformation of the dam body.

Statistical items		Dam body
		Completion period	Storage period
Displacement along the river (cm)	Downstream	−18.2	−10.2
	Upstream	6.9	9.25
Settlement (cm)		47.5	48.8
Major principal stress (MPa)		2.13	2.14
Minor principal stress (MPa)		1.21	1.23

Table 2.

Characteristic values of stress and deformation of the dam body.

Figures 9 and 10 show the contour of the displacements of the dam body during the completion period and the storage period. The simulation results show that the maximum horizontal displacement occurs in the dam body of the old dam and the maximum settlement occurs at the interface between the old and new dams. During the completion period, the maximum settlement of the dam is 47.5 cm, and the horizontal displacement to the upstream and downstream is 18.2 cm and 6.90 cm, respectively. After the water storage, the maximum deformation of the dam body under upstream water load was reduced to 10.2 cm, while the horizontal displacement towards the downstream was increased to 9.25 cm, and the maximum settlement was increased to 48.8 cm.

Figure 9.

Contour of the displacements of the dam body during the completion period (cm). (a) Displacement along the river and (b) settlement.

Figure 10.

Contour of the displacements of the dam body during the storage period (cm). (a) Displacement along the river and (b) settlement.

The results of principal stress calculation show that due to the large modulus of cutoff wall and pile foundation, significant stress concentration has occurred in the dam.

5.2 Stress and deformation of the cutoff wall

Table 3 lists the characteristic values of the stress and deformation of the cutoff wall.

Statistical items		Cutoff wall
		Completion period	Storage period
Dam axial displacement (cm)	Left side bank	/	−0.11
	Right side bank	/	0.12
Displacement along the river (cm)	Downstream	/	10.6
Settlement (cm)		/	0.48
Dam axial stress (MPa)	Tensile stress	−0.21	−2.53
	Compressive stress	1.18	3.21
Major principal compressive stress (MPa)		2.25	12.0
Minor principal tensile stress (MPa)		−0.23	−1.74

Table 3.

Characteristic values of stress and deformation of the cutoff wall.

Since the cutoff wall is constructed after the new dam is filled to 2447.9 m, the deformation of the cutoff wall will not occur during the completion period, so only the deformation distribution during the storage period is given. Figure 11 shows contour of the displacement of the cutoff wall during the storage period. The axial displacement of the dam is represented by the compression from both sides towards the riverbed, and the deformation in the direction of the right bank and the left bank is 0.12 cm and 0.11 cm, respectively. The axial displacement of the dam is generally small. For the displacement along the river, the water load shows a deformation towards the downstream, and the maximum value is 10.6 cm. Because the upper part of the impervious wall is filled with rockfill and supported by the connecting plate, the deformation along the river of the impervious wall increases first and then decreases slightly from the bottom to the top. For the vertical displacement, the maximum value is 0.48 cm, which increases gradually from the bottom to top under the action of the upper water load.

Figure 11.

Contour of the displacements of the cutoff wall during the storage period (cm). (a) Dam axial direction, (b) displacement along the river and (c) settlement.

Figure 12 shows the contour of the dam axial stresses on the downstream and upstream sides of the cutoff wall during the completion period. Figure 13 shows the contour of the dam axial stresses on the downstream and upstream sides of the cutoff wall during the storage period. Because the cutoff wall will be built after the new dam is basically completed, the stress difference between the upstream and downstream faces of the completion period is small, the stress of the cutoff wall is mainly caused by the self-weight, and the tensile and compressive stresses are small. During the storage period, the axial stress of the dam corresponds to the deformation direction. After storage, the upstream face is in tension at both ends of the middle compression zone, while the downstream face is basically in compression, but the pressure stress at both sides is significantly greater than that at the riverbed. The maximum value of the tensile and compressive stress is −2.53 MPa and 3.21 MPa, respectively. For the major principal stress, the downstream stress is greater than the upstream stress because the deformation is oriented downstream during the storage period. At the same time, due to the relatively small height of the wall near the bank slope and the influence of the boundary constraints, the local stress near the bank slope is concentrated, so the stress at the bank slope on both banks is large, and the maximum pressure stress is 12.0 MPa. For the minor principal stress, the upstream and downstream faces are all in compression at the middle part and tension at both sides. The maximum tensile stress is −1.74 MPa.

Figure 12.

Contour of the dam axial stresses on the downstream and upstream sides of the cutoff wall during the completion period (MPa). (a) Downstream side and (b) upstream side.

Considering the ultimate compressive strain of 700 με and ultimate tensile strain of 100 με, the allowable compressive strength and tensile strength of C25 concrete are 19.6 MPa and −2.8 MPa, respectively. From the above calculation results, the tensile and compressive stresses of the cutoff wall are all within the allowable range for C25 plain concrete (Figure 13).

Figure 13.

Contour of the dam axial stresses on the downstream and upstream surface of the cutoff wall during the storage period (MPa). (a) Downstream side and (b) upstream side.

5.3 Stress and deformation of the connecting plate and toe slab

Table 4 lists the characteristic values of the stress and deformation of the connecting plate and toe slab during the storage period.

Statistical items			Storage period
Connecting plate and toe slab	Dam axial displacement (cm)	Left side bank	−0.71
		Right side bank	0.89
	Displacement along the river (cm)	Upstream	/
		Downstream	5.36
	Settlement (cm)		5.63
Connecting plate	Dam axial stress (MPa)	Tensile stress	−1.81
		Compressive stress	0.56
	Major principal stress (MPa)	Compressive stress	0.86
	Minor principal stress (MPa)	Tensile stress	−1.82
		Compressive stress	0.32
Toe slab	Dam axial stress (MPa)	Tensile stress	−4.78
		Compressive stress	1.53
	Major principal stress (MPa)	Compressive stress	6.33
	Minor principal stress (MPa)	Tensile stress	−4.80
		Compressive stress	0.90

Table 4.

The characteristic values of the stress and deformation of the connecting plate and toe slab during the storage period.

Figure 14 shows the contour of the deformation of the connecting plate and toe slab during the storage period. For the axial displacement of the dam, the water displacement is represented by the compression from both sides of the riverbed. The axial displacement of the dam is generally small. The maximum displacements of the left and right banks after water storage are 0.71 cm and 0.89 cm, respectively, which occur in the 0 + 209 and 0 + 65 sections. The displacement of the river is characterized by a downward-directed deformation under the water load during the storage period, with a maximum value of 5.36 cm, which occurs in the 0 + 125 section of the riverbed. For the vertical displacement, the maximum value is 5.63 cm during the storage period, which also occurs at the 0 + 125 section of the riverbed. It can also be seen from Figure 14 that due to the deformation joint between the connecting plate and the toe slab, the connection between the toe slab and the connecting plate is staggered, but the magnitude is small, and the setting of the toe slab length is appropriate.

Figure 14.

Contour of the deformation of the connecting plate and toe slab during the storage period (cm). (a) Dam axial displacement, (b) displacement along the river and (c) settlement.

Figure 15 shows the contour of the dam axial stresses of the connecting plate and toe slab during the storage period. Under the action of water loading, the deformation of the connecting plate is constrained by the toe slab, and the deformation of the toe slab is constrained by the face slab, so the stress of the toe slab is greater than that of the connecting plate. The dam axial stress, corresponding to the deformation direction, is mainly manifested as tension at both ends and compression in the middle, and the downstream compressive stress is greater than the upstream compressive stress. After the storage period, the maximum tensile compressive stress is −4.78 MPa and 1.53 MPa, respectively, which occurs at the right end of the toe slab and in the 0 + 95 section of the riverbed.

Figure 15.

Contour of the dam axial stresses of the connecting plate and toe slab during the storage period (MPa).

Considering the ultimate compressive strain of 700 με and ultimate tensile strain of 100 με for C30 concrete, the allowable compressive strength and tensile strength are 27.3 MPa and −3.9 MPa, respectively. It can be seen from the above calculation results that the compressive stress and tensile stress of the connecting plate and toe slab are within the allowable range for C30 plain concrete, but the maximum tensile stress of the toe plate exceeds the allowable value of C30 plain concrete, and the exceeding area is mainly located in the local area at the junction of the toe slab and the bank slope, which could be resolved by adding reinforcement.

6. Conclusions

The Zhushou Reservoir was transformed from a clay core rock-debris dam to a concrete-faced rockfill dam, with the maximum dam height increasing from 63.4 m to 98.1 m. The three-dimensional finite element method was used to simulate the operation process of construction filling and the storage period, and the conclusions are discussed as follows:

1. The simulation results show that the maximum horizontal displacement occurs in the dam body of the old dam, and the maximum settlement occurs at the interface between the old and new dams. Due to the large modulus of cutoff wall and pile foundation, significant stress concentration has occurred in the dam.

2. During the storage period, the maximum axial tensile and compressive stresses of the cutoff wall are −2.53 MPa and 3.21 MPa, respectively, and the maximum major and minor principal stresses are 12.0 MPa and −1.74 MPa, respectively. The tensile and compressive stresses are all within the allowable range for C25 plain concrete, and the cutoff wall will not be damaged under static conditions.

3. During the storage period, the maximum axial tensile and compressive stresses of the toe slab (connecting plate) dam are −4.78 MPa and 1.53 MPa, respectively, and the maximum major and minor principal stresses are 6.33 MPa and −4.80 MPa, respectively. The compressive stress of toe slab and connecting plate and the tensile stress of connecting plate are all within the allowable range for C30 plain concrete, but the tensile stress of the local area at the junction of toe slab and bank slope has exceeded the allowable value for C30 plain concrete, so the reinforcement should be strengthened at this location.