Coupled Thermo-Hydro-Mechanical Analysis of Valley Narrowing Deformation of High Arch Dam: A Case Study of the Xiluodu Project in China

Abstract

General studies examining reservoir bank deformation during its impoundment primarily consider the coupling effect between the seepage field and the stress field, but thermal field variation in the bedrock and its effect are rarely considered. In this paper, a case study concerning a 285.5 m high arch dam project, where a valley narrowing deformation occurs after the initial impoundment, is implemented. An analysis of in situ measurement is given to interpret the causes of the unique hydro-thermal phenomenon of the project. Possible reasons for the valley narrowing deformation pattern are discussed. A numerical model based on the thermo-hydro-mechanical (THM) coupling theory of porous medium is used to calculate the evolution processes of the thermal, seepage, and stress fields of the area after impoundment of the reservoir. The simulated deformation trend and pattern of the river valley are consistent with the monitoring data. The results demonstrate that water infiltration after impounding cools the bedrock and the temperature decrease makes the bedrock contract, which induces the narrowing deformation of the valley. Factor analysis of the hydrothermal field shows that temperature variation is the main cause of long-term deformation. Thus, it shall be considered as a key factor in terms of structural safety assessment. Furthermore, sensitivity analysis of the hydraulic conductivities of rock strata suggests that future development of the deformation can be eased off if the anti-seepage method is adopted on the bedrock.

1. Introduction

Valley narrowing deformation is a phenomenon that riverbanks deform to their free faces and narrow the valley width. High arch dams, which are mainly located in narrow valleys, are sensitive to this kind of deformation because of their particular shape. Several arch dams were squeezed to crack by the valley narrowing deformation, such as the Beauregard arch dam in northwestern Italy [1] and Zeuzier arch dam in Switzerland [2]. In China, the same kind of deformation has also been reported at the Lijiaxia project [3], Jinping-I project [4], Xiluodu project [5], and several other high arch dams since their initial impoundments. It is a key step to figure out the valley narrowing mechanism to assess the structural safety of a certain dam [6].

Currently, the corresponding mechanisms for valley narrowing deformation of different arch dams are diverse. Classic theory on valley deformation mainly focuses on the slope stability directly affected by stress field changes and seepage field variations [7,8,9]. In the case of the Beauregard arch dam, a deep-seated gravitational slope deformation after impounding was confirmed using long-term monitoring data by integrated techniques [10]. Because of the impoundment, the pore-pressure increased. The increased pore-pressure reduced the effective stresses along the basal shear zone. And the shear strength reduction of the rupture surface induced the plastic deformation of the slope [11]. Research works on the Zeuzier case by Lombardi [2,12] show that the Zeuzier valley was narrowed due to the consolidation settlement of the bedrock, which was induced by tunnel excavation. It was the same case in the Swiss Alps valley. The deformation of several arch dams was observed after the excavation of the Gotthard base tunnel [13,14]. Further studies on horizontal deformation of the Alps valley reveal that cyclic water-table fluctuations in the mountain slopes could induce a cyclic deformation process [15,16,17]. The aforementioned studies examine the mechanism of valley deformation from the perspective of the hydro-mechanical coupling process in the bedrock of arch dams. Both elastic and plastic deformations in the reservoir banks can be induced by the pore-pressure variation because of impoundment, tunnel excavation, or seasonal precipitation changes. In the scope of dam engineering, the valley narrowing deformation is always explained to be a consequence of seepage variation while the effect of thermal field change is often neglected. In fact, the effect of subsurface temperature change on rock deformation has gradually been recognized by studies on tunnels [18], alpine geosystems [19], radioactive waste repositories [20], and many other kinds of fractured rock systems. Gunzburger [21] marks the daily fluctuation of surface temperature as the main cause of the rockfall at Rochers de Valabres slope in France. In the research on a rock slope at Randa (Switzerland), Gischig [22] demonstrates that the thermal effects driven by near-surface temperature cycles can lead to the deformation of deep rock slope. Gruber [23] suggests that a 0.5 °C temperature rise has the potential to induce slope instability.

In Xiluodu, the traditional hydro-mechanical (HM) coupling theory has insufficient explanation because of the project’s unique geological context and deformation pattern. On the one hand, pre-project exploration data spotted the existence of hot spring and abnormal high ground temperature beneath the dam site. On the other hand, the valley narrowing process does not coincide with the reservoir water level or seepage pressure fluctuation. Moreover, the valley narrowing deformation at Xiluodu high arch dam has reached 78 mm in seven years since the initial impoundment, which is near twice as much as the deformation of similar hydropower projects. Thus, new factors shall be taken into account to appropriately explain the deformation mechanisms and evaluate the dam’s long-term safety.

In this paper, basic information and geological features of the reservoir are firstly described and a summary of former analysis [24,25,26] on the geological profile is given to reveal the causes of the project’s geological uniqueness based on observed data (Details in Section 2). After that, seven-year in situ monitoring data of plumb lines and survey lines is closely examined and a unique deformation mode of bedrock is derived (Details in Section 3). Qualitative analysis shows that the valley deformation after the impoundment is consistent along both the river direction and the vertical direction with no sign of interlayer sliding. Based on the analysis above, a three-dimensional coupling thermos-hydro-mechanical (THM) model is reasonably adopted to study the variations of thermal field and seepage field and their effect on the valley narrowing deformation. The numerical results are consistent with displacement and temperature monitoring data recorded at the dam site for over 7 years. The numerical results show that the thermal field variation is the main cause of the long-term deformation after the initial impoundment and it cannot be neglected. The magnitude of valley narrowing is sensitive to the permeability of the bedrock underneath the dam structure. Thus, anti-seepage methods are recommended to reduce further deformation development.

2. Project Overview and Geological Context

Located in the lower reaches of the Jinsha River, Yunnan Province in Southwest China, the Xiluodu hydropower project is China’s second-largest and the world’s third-largest hydropower project, with the total capacity up to 1.267 × 10¹⁰ m³. It is a double curvature arch dam that is 285.5-m high with a 200-m radius. The absolute altitude of the dam’s foundation is 324.5 m. By contrast, the designed normal pool level of the reservoir is 600 m and the corresponding valley width is approximately 500 m. Thirty-one monoliths are planted across the river. The transversal section of the dam is curve-shaped, 14 m wide at the crest, and 60 m wide at the footing. As a typical project for China’s 300 m-high arch dam, it began casting concrete in March 2009. The reservoir began its initial impoundment in March 2012 when the dam was under construction [27]. In September 2014, the reservoir reaches the designed storage water level.

The reservoir area of the dam project is located in the center of a synclinal basin surrounded by limestone outcropped areas as shown in Figure 1. The natural water level of the river varies from 370 m to 390 m above sea level. Geological exploration shows that the river valley is the main outlet of the subsurface runoff.

As shown in Figure 2, the reservoir banks are steep and the valley section is nearly U-shaped. Both sides of the valley slope and the river bedrock at the dam site are mainly composed of Permian basalt strata, which are hard, intact, and suitable for arch dam construction. On top of the basalt strata is Permian sandy shale. Approximately 500 m thick Permian limestone stratum lies under the Permian basalt strata with a shale sediment layer between them.

Geological exploration before the construction of the dam showed signs of hot springs and a relatively high geothermal field at the dam site. Temperature measurements in the boreholes revealed that the ground temperature in limestone could be up to 49.6 °C as shown in Figure 3. The average ground temperature in basalt was around 35 °C while the water temperature in the Jinsha River varied from 8 °C to 22 °C.

Researchers have tried to explain the high geothermal phenomenon in the aspect of rock permeability [25], groundwater circulation [24], and chemical water-rock interaction [26].

According to the previous studies, the river valley is the main outlet of three different strands of groundwater. The sand shale (P_2x) has a low vertical permeability, which means the surface water can only run down along the hillside without infiltrating into the basalt. Thus, the initial groundwater level in basalt is lower than 400 m and mainly controlled by the variation of the river level. The shale sediment layer (P₂β_n) is a fine aquifuge, causing hydraulic head discontinuity between basalt and limestone. The supply region of groundwater in underlying limestone is its stratum outcrops around the basin where the altitude is above 2000 m. The subsurface runoff is heated up by geothermal flux, which causes geothermal springs and relatively high ground temperature at the dam site. The schematic diagram of groundwater circulation is shown in Figure 4.

3. In Situ Measurement and Analysis on Deformation Pattern

3.1. In Situ Measurement Overview

3.1.1. Valley Narrowing Process

Seven survey lines have been set up to monitor the valley deformation at the dam site. The survey line measures the distance between two points on both sides of the valley, i.e., the valley width, by the laser range finder. The information of the survey lines is shown in

Table 1. Information on the survey line to monitor valley deformation.

ID	Elevation/m	Length/m	Distance to the Dam Axis/m	Position
VD01	749.00	791.98	350.00	Upstream
VD02	611.00	593.19	260.00	Upstream
VD03	722.00	786.71	50.00	Upstream
VD04	611.00	568.27	50.00	Upstream
VD05	707.00	632.19	300.00	Downstream
VD06	610.00	524.55	350.00	Downstream
VD07	561.00	415.79	350.00	Downstream

.

Figure 5 shows the location and deformation process of seven survey lines for the valley width with the reservoir level-rising process. The valley width continues decreasing since the beginning of the impoundment. A special trend is observed: All seven survey lines share similar tendencies and deformation amplitudes, regardless of their installation position and elevation. The deformation process has slowed since the first flood season but has not ceased until now. The latest monitoring data shows that the shrinkage deformation is up to 78 mm from 20 December 2012 to 26 July 2019 on average.

3.1.2. Seepage Pressure Process

In the spring of 2019, several boreholes were carried out and seepage measurements were implanted at the bottom of each borehole to obtain the groundwater pressure in limestone. The spatial locations and seepage pressure progress of these boreholes are shown in Figure 6.

3.1.3. Thermal Field Process

Figure 7 shows two typical curves of upstream and downstream bedrock temperature. The temperature has dropped to 17.5 °C as the reservoir impounded, and the descending trend has not fully converged yet. The previous studies [28,29] have shown that the water temperature at the bottom of the reservoir stabilizes at about 14.5 °C. Therefore, the cooling process of bedrock is still in progress.

3.2. Bedrock Deformation Results

As shown in

Table 2. Correlation coefficients between reservoir water level and groundwater level in limestone.

ID	Correlation Coefficients
ZK01	0.9995
ZK02	0.9864
ZK03	0.9893
ZK04	0.9970
ZK05	0.9976
ZK06	0.9920
ZK07	0.9900

As the results turn out, there is a positive correlation between reservoir water level and groundwater pressure in limestone. This indicates that the reservoir water may have been infiltrated into the limestone.

, a strongly positive correlation between reservoir water level and groundwater pressure in limestone is observed while the valley narrowing process has not totally converged yet. Meanwhile, the temperature beneath the dam site decreases to 17.5 °C and is still in progress. Based on the comprehensive consideration of these tendencies, bedrock deformation is derived from the deformation monitoring system and closely examined.

Seven plumb lines and several external deformation supervision points are installed inside the dam to monitor the deformation of the dam. As shown in Figure 8, the red vertical lines represent the plumb lines implanted in the dam, while the blue circles represent external deformation supervision points.

As plumb lines can monitor the relative deformation between anchoring joints and monitoring points, it is possible to derivate the deformation pattern of bedrock by combining the plumb line with external deformation supervision measurements. More details about the calculation method have been discussed by [30].

The result of bedrock deformation is shown in Figure 9. Two deformation characteristics shall be highlighted. Firstly, the transverse deformation along the vertical direction remains constant, which excludes the possibility of interlayer-gliding. Secondly, the deformation of the riverbed is larger than the dam abutment, which indicates that the limestone underneath the riverbed is the main deform zone, rather than the slope.

These are the most unique and important features that have rarely been observed in other similar projects. As many cases can be interpreted as landslides or local surficial movements due to changes in the seepage field. The main deform zone is always the river valley itself. Thus, a new physical model of the valley narrowing deformation shall be proposed to explain the deformation pattern. The effect of thermal field variation has been gradually recognized based on the comprehensive consideration of the abnormal geothermal field before the impoundment and the unique deformation features after the impoundment.

3.3. Qualitative Analysis of the Deformation Mechanism

The novel deformation process of Xiluodu River valley since the initial impoundment has drawn wide attention from researchers and engineers [31,32,33]. It affects the deformation state of the arch dam and threatens the long-term security of the project [34]. Using hydrogeological profiles and analyzing the monitoring deformation process, we list and discuss several potential driven mechanisms for the valley deformation in this study.

3.3.1. Plastic Deformation along the Interlayer Faces of Basalt

For the saturated rupture surface, the increasing pore-pressure reduces the effective normal compressive stress and its shear strength [35]. The reduction in shear strength can induce landslide, e.g., the Beauregard case [1]. For the unsaturated rupture surface, which becomes submerged after the reservoir impoundment, complex physical and chemical reactions may occur between the infiltrated water and the rock mass and they could weaken the rock mass [36,37]. The results in both cases are irreversible plastic deformation of the slope. For the Xiluodu case, the class of unweathered basalt is II and the class of interlayer faces of basalt is III according to the RMR system. The range of the drill core RQD (%) of unweathered basalt is 80–90. The range of drill core RQD (%) of interlayer faces is 50–80. The orientation of the rock strata is 3°–7°. The rock strata at the dam site are intact and nearly horizontal according to pre-project exploration. The monitored deformation pattern during the initial impoundment shows no obvious jump across interlayer weak planes, which should be rupture surfaces (Figure 9). Thus, the possibility of plastic deformation along the rupture surface in the dam abutment is low.

3.3.2. Confined Aquifer Influence

The groundwater level in limestone rises by 100 m after the reservoir impounds, which makes the explanation of confined aquifer influence tempting. However, the groundwater head fluctuates with reservoir water level (Figure 6) while valley deformation does not (Figure 5). Besides, the stratum is horizontal at the dam site. The confined water pressure shall lead to uplift deformation rather than valley shrinkage. Therefore, the possibility of a confined aquifer influence is relatively low.

3.3.3. Thermal Field Variation

A comparison of bedrock temperature before and after reservoir impoundment demonstrates that the reservoir water affects not only the seepage field but also the thermal field. The valley contraction and bedrock temperature drop occur simultaneously. Based on those observed phenomena, we proposed a valley deformation mechanism driven by thermal field variation. The basic principle is shown in Figure 10.

As the water level rises, the reservoir water with low temperature infiltrates into the rock mass, which originally has a higher temperature as described above. The cooling effect is likely to induce rock contraction at the bottom of the valley and consequently narrows the valley.

4. Numerical Simulation and Discussion

The behavior of dam-reservoir-foundation systems, particularly their fluid transport and nonlinear mechanics characteristics, is affected by complex interactions among multiple physical and geochemical processes [38,39]. Dam-reservoir-foundation systems are designed to raise the water level for flood control and hydropower generation. The rising water level in the reservoir causes tremendous changes in pressure and temperature, which induce dynamic redistribution of the in-situ stress. Subsequently, hydraulic properties such as permeability and porosity are affected by the geometrical changes of fractures and pore spaces. Thus, interactions occur between fluid and rock. Additional coupling occurs through pressure- and temperature-dependent material properties.

4.1. Governing Equations and Constitutive Models

The physical system of dam-reservoir-foundation incorporates a non-isothermal fluid flow in a deforming porous medium. The governing equations are derived from the local equilibrium of thermodynamics and macroscopic balance equations. Thus, the liquid phase pressure p, temperature T, and displacement vector U are the corresponding dependent variables.

4.1.1. Seepage Flow

The seepage flow process in a saturated porous medium is described by the following mass balance equation with Darcy’s law,

$$\nabla ·{\rho }_w\left[-\frac{\kappa }{\upsilon }\left(\nabla p+{\rho }_{w}g\nabla D\right)\right]+{\rho }_{w}S\frac{\partial p}{\partial t}+{\rho }_w\frac{\partial {\epsilon }_{vol}}{\partial t}=Q,$$

(1)

where ${\rho }_w$ is the water density, $\kappa$ is the intrinsic permeability of rock mass, $\upsilon$ is the water viscosity, $g$ is gravity, D is the vertical coordinate, S is the storage coefficient, ${\epsilon }_{vol}$ is the volumetric strain of the rock mass, and Q is a fluid source/sink term [38].

The specific storage coefficient S can be interpreted as weighted compressibility of the bulk aquifer material and fluid in the pores. Thus, S is derived as

$$S=\left(\alpha -\varphi \right)/K_s+\varphi /K_w,$$

(2)

where α is Biot–Willis coefficient, $\varphi$ is the porosity, Ks is the solid bulk modulus (the inverse of the fluid compressibility), and Kw is the water bulk modulus [40].

The fluid viscosity significantly affects the fluid flow and strongly depends on the temperature [41]. To consider the effects of temperature variation on seepage flow, a fitting function used to describe the dependency relation is expressed as [42]

$$\nu =2.414\times {10}^{-5+247.8\left(T-140\right)},$$

(3)

where T denotes the temperature in Kelvin.

4.1.2. Heat Transfer

The heat transfer process in a saturated porous medium consists of convective and diffusive fluxes. The energy balance equation for heat transfer is

$${\left(\rho C\right)}_{eff}\frac{\partial T}{\partial t}+{\rho }_w{C}_{w}u·\nabla T+\nabla ·\left(-{\kappa }_{eff}\nabla T\right)=0,$$

(4)

where $\varphi$ is the porosity. ${\left(\rho C\right)}_{eff}=\varphi {\rho }_w{C}_w+\left(1-\varphi \right){\rho }_s{C}_s$ is the equivalent volume heat capacity, $C_w$ is the specific heat capacity of water, $C_s$ is the specific heat capacity of solid, and ${\rho }_s$ is the solid density [43]. The fluid velocity u is included in the convection term. ${\kappa }_{eff}=\varphi {\kappa }_w+\left(1-\varphi \right){\kappa }_s$ denotes the equivalent heat conductivity of a porous medium, ${\kappa }_w$ is heat conductivity of water, and ${\kappa }_s$ is the heat conductivity of solid.

4.1.3. Equilibrium Equation

The deformation process of a saturated porous medium is assumed to be isotropic linear elastic. Considering the contribution of the pore-pressure with the principle of effective stress and thermal stress, the local equilibrium equation in terms of Lame constant can be written as

$$\left(\lambda +\mu \right)U_{j,ii}+\mu U_{i,jj}+F_i+{\left(\alpha p\right)}_{,i}-\left(3\lambda +2\mu \right)\beta T_{,i}=0,$$

(5)

where λ and μ are Lame constants, F is the body force, α is Biot–Willis coefficient, β is the thermal expansion coefficient, and U_j is the displacement component [40].

Equation (5) describes an equilibrium state which neglects the inertial effects. Its application to the case of time-dependent flow and heat transfer model is valid because the time scale of the structural response is generally many orders of magnitude faster than the time scale of the flow and heat transfer. Thus, even the change in stress and strain with time in equation appears stationary, and the structure-to-fluid coupling term, which involves the rate of strain, is nonzero [40].

4.2. The Definition of Model Parameters

4.2.1. Spatial Discretization

Based on the hydrogeological profiles of the reservoir area, a finite element model is developed considering the rock strata distribution and main hydraulic structures, which include the arch dam, underground powerhouses, and their grouted curtain, as shown in Figure 11. The defined model area is 5 km wide perpendicular to the river at the dam site and 4 km long along the river. The bottom of the model is set to −1 km above sea level. Due to the irregular terrain, the model geometry is discretized with 155,951 tetrahedrons and 5508 triangular prisms.

4.2.2. Boundary Conditions

Prior to the numerical calculation, boundary conditions of the reservoir area should be defined according to the pre-project hydrogeological explorations in Section 2. For the mechanical boundary conditions, normal constraints are assigned to the truncated boundaries of the simulated area. A 400 m-head is provided to the vertical boundaries around the model based on linear interpolation of the measured seepage pressure. An impermeable boundary is also assigned to the bottom of the model. The hydraulic head in the river valley is set at 370 m before impounding to simulate a stationary state. After impoundment, the hydraulic head of the upstream valley in the reservoir is set to be the dynamic reservoir level ($h_{rl}$) based on the daily water level measurement which has been shown in Figure 5 and Figure 7. The downstream valley remains unchanged. A heat-insulated boundary condition is assigned to the vertical boundaries. A constant temperature of 45 °C is set for the bottom of the model, according to a former study on the geothermal field [24,25]. The top surface of the model is set to the mean annual temperature, except the river valley is set to the water temperature (14.5 °C) according to the former study [28,29]. The schematic boundary conditions are shown in Figure 12.

4.2.3. Material Parameters

According to field test data, the basic material parameters of the numerical model are shown in

Table 3. Material parameters of the finite element method (FEM) model.

Parameter	Dam	Grouted Curtain	Shale	Limestone	Basalt	Weathered Zone
Density (kg/m3)	2.40 × 103	2.20 × 103	2.60 × 103	2.60 × 103	2.65 × 103	2.60 × 103
Thermal conductivity (W/(m·K))	1.80	2.10	2.72	2.90	4.65	4.65
Specific heat capacity (kJ/(kg·K))	0.88	0.85	0.86	0.86	0.86	0.86
Thermal expansion (1/K)	7.5 × 10−6	8.0 × 10−6	8.0 × 10−6	9.5 × 10−6	8.0 × 10−6	8.5 × 10−6
Elastic modulus (GPa)	43.0	15.0	10.0	20.0	25.0	15.0
Possion’s ratio	0.17	0.20	0.22	0.30	0.20	0.20
Porosity	0.040	0.10	0.25	0.20	0.050	0.20
Biot’s coefficient	-	0.80	0.85	0.95	0.80	0.95
Hydraulic conductivity (m/d)	-	1.0 × 10−4	2.0 × 10−2	5.0 × 10−2	2.0 × 10−2	Table 4.

. Specifically, the hydraulic conductivity of the weathered zone near the river valley varies dramatically with respect to the buried depth, as shown in

Table 4. Hydraulic conductivity of the weathered zone in the river valley.

Weathered Basalt			Weathered Limestone
Buried Depth (m)	Hydraulic Conductivity (m/d)		Buried Depth (m)	Hydraulic Conductivity (m/d)
	Horizontal	Vertical		Horizontal	Vertical
0–40	2.80–1.40	0.84–0.35	0–30	30.0–10.0	20.0–15.0
40–90	1.62–0.60	0.55–0.15	30–80	20.0–5.0	10.0–2.0
90–190	0.54–0.10	0.17–0.02	80–120	10.0–1.0	5.0–0.51
>190	0.06–0.02	0.015–0.003	120–300	3.6–0.5	2.0–0.2

.

4.3. Numerical Results

The multi-field response of the study area was obtained by using the transient finite element solver in COMSOL Multiphysics^®. The calculation is divided into two stages: A stationary state before impoundment and a transient state after impoundment. The pre-project exploration data are used to calibrate the model in the stationary state. Next, the transient state of the study area is calculated with the initial multi-field distribution, which is obtained from the calibrated stationary state results. Two specific sections shown in Figure 13 are used to represent the simulation result in the following chapters.

4.3.1. Thermal and Seepage Field before Impoundment

The simulated thermal and seepage fields are shown in Figure 14. The calculation result fits the previous study on the hydrogeological context which has been discussed in Section 2. As it has been pointed out, the upward water brings heat from limestone and causing abnormally high ground temperature near the dam site. Simulated temperature distribution near the dam site is extracted and compared to the measured values. As shown in Figure 15, the calculated multi-field distribution of the stationary state is consistent with the hydrogeological exploration results in Figure 3. It has to be noted that the underlying limestone stratum with higher permeability is the key, which plays as a seepage channel and makes the upward seepage flow happen.

4.3.2. Thermal and Seepage Variation after Impoundment

The reservoir began the initial impoundment at the end of 2012, and the reservoir level rises from 370 m to 600 m. Because of the impoundment, the hydraulic head in the rock mass near the valley increases drastically to match the reservoir level. As shown in Figure 16, the previous seepage direction at the bottom of the valley was reversed. Though the waterproof curtain can reduce the horizontal seepage flow. It is not able to prevent infiltration along the vertical direction. Thus, the water in the reservoir will infiltrate downward to the limestone and then bypass the waterproof curtain. The infiltration of cold water will form a cooling zone underneath the river bed which leads to the valley shrinkage.

Comparisons between simulated and measured data are as follow. The simulation result of thermal variation extracted from point 1 is consistent with the thermal field measurement represented in Figure 17. The numerical result shows that the temperature will drop by 1 °C in the next two years before it reaches an equilibrium state.

4.3.3. Valley Narrowing Process after Impoundment

As has been previously discussed in Figure 10, the infiltration of cold water will form a cooling zone underneath the riverbed. And the cooling zone will lead to a valley shrinkage which is consistent along both the vertical direction and the river flow direction. Figure 18 shows the distribution of temperature decrease and the displacement vector of the rock mass near the river valley. The maximum temperature decrease was approximately 16 °C at 80–160 m below the river bed. The rock mass cooled and contracted, induced a decrease in valley width. The induced displacement pattern is consistent with the monitored horizontal deformation of the abutment which has been discussed in Figure 9.

Comparison between observed data and simulated data is shown in Figure 19. The simulated results are consistent with the observed data from survey lines which has been shown in Figure 5. According to the simulated results, an average deformation of 15 mm will take place in the next 5 years. The predicted convergence values of survey lines are between 79 mm and 98 mm, which is relatively large but still within the design criteria. After 5 years, the annual mean deformation rate will be less than 0.005 mm/day and meet the national standard of deformation convergence.

4.4. Discussion

4.4.1. Factor Analysis

Comparative studies are also implemented in parallel which considers the hydro-mechanical (HM) interaction process and neglects the influence of thermal field variation. The deformation processes at VD4 and VD7 are displayed in Figure 20.

A comparison between THM and HM results represents the influence of the thermal field. According to the results, thermal field variation shall be responsible for more than 60 percent of the total deformation; it is the main factor that causes the deformation after the initial impoundment. While long term deformation is mainly caused by the temperature changes, the deformation in the first two years is mainly affected by seepage variation.

4.4.2. Sensitivity Analysis of Hydraulic Conductivity

The hydraulic conductivity of the rock mass is often difficult to determine because of the complexity of its intrinsic fractured network. Practically, an estimated range of conductivity is often provided for a specific rock stratum instead of an exact value, as shown in Table 4. A sensitivity analysis of the valley deformation is performed considering the hydraulic conductivities of the weathered basalt and limestone strata near the valley.

The adopted conductivity values are the upper limit, lower limit, half of the lower limits of the suggested ranges in Table 4. Thus, 9 scenarios with different conductivity combinations were implemented. According to the result in Figure 21, a nearly 60% reduction for valley deformation was observed when the conductivities of the weathered rock strata changed from their upper limits to lower limits.

5. Conclusions

In this paper, a thermo-hydro-mechanical coupled approach is used to study the valley narrowing deformation in the Xiluodu project. Similar valley narrowing deformation caused by thermal contraction of bedrock has rarely been reported before. Yet the particular combination of similar hydrogeology conditions and the super-high arch dam project is not common. It should be noted that, because of the extreme complexity of fractured rock mass, the actual cause for the valley narrowing deformation of the Xiluodu arch dam remains an open question, and the research work in this paper is supposed to provide a new perspective. The main results are as follows:

• The proposed THM coupling model explains the valley narrowing deformation mechanism of the Xiluodu Project. According to the analysis of in situ measurement and numerical results, high ground temperature in the river valley before the impoundment is mainly induced by the deep groundwater circulation, which continuously brings heat to the bedrock at the bottom of the valley. After impoundment, the previous seepage direction at the bottom of the valley is reversed. The reservoir water with lower temperature cools the rock mass at the bottom of the valley and thermal contraction occurs. As a result, the “U” shaped valley narrows.
• Seepage field variation is the main factor of the deformation in the first two years, while thermal field variation causes long-term deformation after the initial impoundment. It shall be responsible for more than 60 percent of the total deformation. The numerical results also predict that 15 mm deformation will occur in the next five years. By 2024, the annual deformation rate will be less than 0.005 mm/day.
• Based on the result of sensitivity analysis on hydraulic conductivity, the conductivity of the bedrock is one of the key factors that influence the magnitude of valley narrowing deformation. Hence, anti-seepage methods such as high-pressure grouting are recommended to decrease the heat exchange underneath the riverbed and reduce further development of valley deformation in the feature.

Though the converged value of deformation is relatively large, it is still controlled within the design criteria. Dam safety is guaranteed. As there are many high dam projects still under construction or in planning in China, the analysis method proposed in this article provides reference meaning for future studies on the multi-field evolution process of dam-reservoir-foundation systems, particularly thermal field variation. For future studies on the Xiluodu project and similar cases, hydro-chemical experiments and THMC coupling simulation can be very helpful to gain a more comprehensive understanding of the mechanisms. A recent study [44] on the model has taken fracture into account due to its influence on the multi-field coupling. Fracture scanning through boreholes and THM experiment on rock samples are in progress to reveal the impact of THM process on the propagation of natural fractures. The application of a new model will be made possible by obtaining these experimental data. Further research is also necessary, as a full understanding of the valley narrowing mechanism will be an important reference for dam construction and safety assessment.