Flume tests on the hydraulic/mechanical behavior of the Tangjiashan landslide dam considering the inuences of the different debris materials in heterogeneous strata

Landslide dams (LDs) usually form from natural debris materials and exhibit heterogeneous strata along both the depth and run-out directions. In addition, an LD usually has a weaker structure than that of undisturbed ground and is more vulnerable to seepage loading. Considering that the surface layer of naturally packed LD materials is generally in an unsaturated state, it is undoubtedly important to investigate the stability of the unsaturated debris materials in the heterogeneous strata of LDs. In this paper, a systematic ume test program was rst conducted, in which the Tangjiashan LD was carefully referenced for model design. Three water level rising rates and two stratal arrangements were considered in the ume tests. Then, soil-water-air coupled nite element analyses were conducted to simulate the ume tests, and all the material parameters of the LD materials were carefully determined based on the results of the element tests. A comparison of the test and calculated results shows the possibility of using the proposed numerical method to estimate the occurrence of dam breaching and the risk of LD failure. Moreover, the hydraulic/mechanical behaviors of the LD materials and the heterogeneous strata of the LD were very important to the stability of the Tangjiashan LD. Finally, from an engineering viewpoint, the possibility of utilizing a naturally formed LD and thus not destroying it when it forms is also discussed, e.g., dam breaching risk can be reduced by excavation of a drainage tunnel, and the dam stability can be carefully estimated based on accurate geological data.


Introduction
A landslide dam (LD) usually forms from natural debris materials from landslides, avalanches, debris ows, etc., that are often triggered by earthquakes, rainfall, snowmelt, etc. (Ermini and Casagli, 2003;Peng and Zhang, 2012a). Once an LD is formed, water impoundment will result in a landslide-dammed lake (LD lake). An increase in water level can lead to LD failure within a short time after its formation, which threatens the safety of people and property downstream (Zhu et al., 2003). The failure modes of LDs generally include overtopping, piping and sliding failure (Costa and Schuster, 1988), which are similar to those of arti cial soil and rock ll dams. However, unlike arti cial dams, LD debris materials have a wide grain size distribution from millimeters to meters, even including gap grading with high erodibility (Ermini and Casagli, 2003). In addition, due to the rapid sliding and packing of the debris materials, the strata of LDs are often only partially destroyed and basically are heterogeneous along both the depth and run-out directions (Chang and Zhang, 2010). Therefore, it is undoubtedly important to investigate the stability of LDs considering the in uence of the heterogeneous strata of LDs.
Model testing has become an essential method with which to study the stability of LDs, as in the works by Gregoretti et al. (2010), Awal et al. (2011), Chen et al. (2015), Okeke and Wang (2016), Jiang et al. (2018Jiang et al. ( , 2019 and Dhungana and Wang (2019). In these studies, however, the grain size distribution of the LD material was not considered, and the model LD could only re ect part of a real LD with limitations.
Nevertheless, only a few attempts have been made to make the results of LD model tests more credible. Wang et al. (2018) conducted a large-scale model test to identify the premonitory factors of LD failure induced by seepage using material sourced from the Mihata LD. Xiong et al. (2018) selected three typical LD materials to investigate the failure mechanism of unsaturated LDs under seepage loading. Peng et al. (2019) carried out a series of wave ume tests on the erosion failure modes of LDs, and the typical gradation curve of the Donghekou LD was selected to prepare the dam material. In particular, based on eld investigations of an LD in the Xiaoqinling gold mining area, ume tests were conducted with mine waste material, from which particles with diameters >50 mm were removed .
In the abovementioned studies, however, all LD models were uniformly prepared, and the presence of heterogeneous strata in LD was ignored. Therefore, in this research, a systematic ume test program was conducted, and a speci c LD, named the Tangjiashan LD, was referenced for model design. Particular attention was paid to the in uence of the model ground made of different debris materials and heterogeneous strata from the Tangjiashan LD.
Notably, the Tangjiashan LD did not naturally fail, it was partially triggered due to the erosion of overtopping water after digging a spillway to reduce the water level of the LD lake. Therefore, many studies have mainly focused on the breaching process, as in the works by Liu et al. (2009), Chang and Zhang (2010), Fan et al. (2012), Peng and Zhang (2012b) and Shi et al. (2015b). Basically, an LD has a weaker structure than that of the undisturbed rock mass before landsliding and is more vulnerable to seepage loading. Because it is very complicated to estimate the stability of an LD under seepage loading (Xiong et al., 2018), only a few studies have focused on the seepage stability of the Tangjiashan LD. Shi et al. (2018) proposed a coupled computational uid dynamics and discrete element method (CFD-DEM) model to simulate seepage tests on the material of the Tangjiashan LD, in an attempt to reveal the seepage failure mechanism from macro and micro perspectives. Based on accurate geological data from the Tangjiashan LD, Hu et al. (2012) simulated the seepage ow of the Tangjiashan LD under four water level conditions and discussed its seepage stability. However, the surface layer materials of LDs are usually in an unsaturated state, and their strengths are dependent on the degree of saturation, which has not been properly considered in past studies related to the seepage stability analyses of the Tangjiashan LD.
Therefore, this research aimed to study the stability of an unsaturated LD under seepage loading conditions while considering the in uence of different debris materials in the heterogeneous strata of the LD. A systematic ume test program was conducted rst, in which the Tangjiashan LD was referenced for model design. The development of the phreatic line in the model LD was carefully measured using sensors, including piezometers, water content sensors and tensiometers. Then, in the corresponding simulation, a saturated/unsaturated soil constitutive model and a deformation-dependent water retention curve (WRC) model (Xiong, 2020) were selected. The parameters of the debris materials of the LD were carefully determined using the results of element tests Xiong et al., 2018;Ma et al., 2020), based on which the ume tests on the hydraulic/mechanical behavior of the Tangjiashan LD were simulated with soil-water-air coupled nite element analyses.

Geological Background Of The Tangjiashan Landslide Dam
On May 12, 2008, a 8.0 M s earthquake occurred in Wenchuan County, Sichuan Province, China, and triggered at least 257 LDs (Cui et al., 2009). The Tangjiashan LD was the largest LD located on the Tongkou River 6 km upstream of Beichuan County, as shown in Photograph 1(a). The Tangjiashan LD was 803.4 m in longitudinal length, 611.8 m in maximum transversal width and 82.7-124.4 m in height.
The upstream slope of the Tangjiashan LD was 18°-22°, and its average downstream slope was 38° (Liu et al., 2009). The estimated volumes of the accumulated debris mass and LD lake were 2.04×10 7 m 3 and 3.16×10 8 m 3 , respectively, posing a great threat to Beichuan County.
According to site investigation results, the Tangjiashan LD maintained a part of the original stratum structure (Fig. 1). It mainly consists of four materials: (1) gravely soil, distributed on the top with a thickness of 5-15 m (Photograph 2); (2) strongly weathered cataclasite, a kind of cohesive granular fault rock, distributed in the middle with a thickness of 10-30 m (Photograph 3); (3) weakly weathered cataclasite, distributed on the bottom with a thickness of 50-67 m (Photograph 4); and (4) silty sand, mainly distributed on the upstream slope of the dam (Photograph 5).
After the formation of the LD, the water level of the Tangjiashan LD lake rst rapidly increased and then gradually increased, as shown in Fig. 2. During the rise of the water level, concentrated seepage points with stable ow were observed at elevations of 669 m and 700 m on the downstream slope of the dam. The water from all these seepage points was clear, indicating that piping did not occur in the dam body.
Considering the continuous increase in water level, a spillway with a depth of 12 m was constructed across the Tangjiashan LD (Shi et al., 2015b). Dam breaching nally occurred on June 7, 2008 due to overtopping. The majority of the dam mass, especially the weakly weathered cataclasite layer, however, remains stable, as shown in Photograph 1(b).
In brief, three important characteristics of the Tangjiashan LD were observed: (a) the debris materials of the Tangjiashan LD are obviously zoned, which can satisfy the purpose of this research; (b) before breaching, seepage did develop in the dam body, showing the necessity to clarify its in uence on the Tangjiashan LD; (c) the debris materials of the Tangjiashan LD were not completely breached. These three characteristics were carefully considered in designing the model LD to investigate the in uence of different debris materials in heterogeneous strata on LD stability using ume tests and corresponding numerical simulations.

Ground materials for ume tests
As mentioned before, the debris materials of LDs exhibit a wide range of grain size distributions. The proper selection of debris materials in ume tests is a very important issue. Due to the lack of information about the grain size distribution curves of debris materials in Tangjiashan LD, three typical debris materials of LDs were utilized as the test materials in the ume tests and numerical calculations of this research. Fig. 3 shows the grain size distribution curves of these debris materials. Element tests have been conducted on these materials, including triaxial tests, water retention tests, and permeameter tests, and satisfactory results have been obtained Xiong et al., 2018;Ma et al., 2020).
To assure the operability and repeatability of the ume tests and their corresponding numerical calculations, three different manmade LD debris materials were prepared by mixing different silica sands with grading curves, as shown in Fig. 3. Without losing the general properties of the debris materials, such as permeability and strength, particles greater than 40 mm were removed and not used in the model tests.
These three grading curves represent (1) sand, (2) a well-graded mixture and (3) a gap-graded mixture. For simplicity, hereafter, the manmade LD debris materials are called LD materials.
By comparing the hydraulic/mechanical behaviors of the three typical LD materials Xiong et al., 2018;Ma et al., 2020) with those of the Tangjiashan LD materials, the sand, well-graded mixture and gap-graded mixture were selected to represent the silty sand, the strongly weathered cataclasite, and the gravely soil of the Tangjiashan LD, respectively. Meanwhile, bricks with high strength and high erosion resistance properties were used to represent the weakly weathered cataclasite in this research.
Permeameter tests on these LD materials under a saturated condition were conducted by Shi et al. (2015b), who found that the permeability of sand, well-graded and gap-graded LD materials in a saturated state were k sw = 8.72×10 -4 , 3.33×10 -4 and 1.29×10 -4 m/s, respectively, with an average void ratio of e 0 =0.48.

Setup of ume tests and LD model
A ume with a rectangular section was utilized in this paper; the ume was 50 m in length, 0.8 m in width, 1.25 m in height and had a 0° slope angle. More details about the ume can be found in Xiong et al. (2018).
The longitudinal section of the Tangjiashan LD along the river, as shown in Fig. 1, was simpli ed with layered strata, as shown in Fig. 4. The height-to-length ratio of the model LD was small, only 2:9, which is an important characteristic of LDs and distinctly different from the geometry of arti cial dams (Costa and Schuster, 1988). Based on the maximum dry density obtained from the compaction tests, the initial dry density of the ground materials in the model LD was determined and is presented in Table 1. The designed compaction degree of the sand and gap-graded mixture on the surface of the dam model was 75%, and the designed compaction degree of the well-graded mixture in the interior was 80%. According to the Froude similarity law, when the geometric scale is 1:300, the time scale is , the velocity scale is , and the density scale is 1 in a 1g gravitational eld.
As shown in Fig. 5, the measuring system consists of three piezometers, ve water content sensors, ve tensiometers, three video cameras, two displacement transducers and three data acquisition instruments, which were used to monitor the deformation and seepage of the model LD. The water content sensors were EC-5 sensors, made by METER Group, Inc. USA, and can measure volumetric water contents from 0% to 100%. The tensiometers were ML-2100AM6 tensiometers, made by Mol, Inc., Japan, and can measure pore pressures from -100 kPa to +100 kPa. To reduce the in uence of the tensiometer as an alien component of the model ground on the seepage in the dam body, an L-shaped plastic plate was xed in the ume, as shown in Fig. 6. Through small holes on the vertical plastic plate, tensiometers can be arranged more reliably.

Test methodology
To investigate the stability of the LD under seepage loading conditions, four ume tests were conducted considering the different debris materials in heterogeneous strata of the LD. Based on the simpli ed model LD shown in Fig. 4, the layers of the gap-graded mixture, well-graded mixture and brick had different thicknesses in different cases. As shown in Table 2, the layer thicknesses of the different materials in Case A-1 and Case A-2 are the same as the corresponding layer thicknesses of the Tangjiashan LD after scaling down the observations 300 times. Meanwhile, in Case B-1 and Case B-2, the layer thicknesses of the gap-graded mixture and well-graded mixture were increased to 120 mm to investigate their in uence on the development of the phreatic line.
Moreover, the in uence of the water level rising rate (WLRR) on the hydraulic/mechanical behavior of the LD was also studied, as shown in Table 2. Considering the control accuracy of the in ow valve, the WLRR in Case A-2 was set to 3 mm/min. The WLRR in Case A-1 was set to three times faster than that in Case A-2. Moreover, the WLRRs in Case A-2 and Case B-1 were the same. Considering the changes in the WLRR of the Tangjiashan LD shown in Fig. 2, the water level in Case B-2 was set to increase at 3 mm/min to 160 mm and then increased at 1 mm/min to the crest of the model LD.
In the tests, the LD materials were rst prepared by mixing dry silica sands with different grain sizes to reach the prescribed grain size distribution curves, as shown in Fig. 3. Then, the bricks were stacked in the shape shown in Figure 4, and the surface layer of bricks was mortared with a mixture of water-cementclay at a mass ratio of 1:2:2 to reduce its permeability. The other parts of the model LD were prepared by compacting every layer to a 40-mm thickness according to the sizes and compaction degrees shown in Fig. 4, Table 1 and Table 2. After that, an initial water level of 50 mm was added on both sides of the model LD, and the measurement system was started. After 1 h, the data from the piezometers were stable, and all the measuring devices were set to zero, which was regarded as the beginning of the test. Next, the water level was raised gradually on the left side of the model LD according to the rate listed in Table 2. Simultaneously, all the measuring devices were started to record the acquired data at regular intervals. Except for Case B-2, the water level upstream was set to increase continually until the model LD failed by overtopping. Considering the situation that the model LD did not fail after overtopping, the in ow valve was closed after the water level reached the crest of the model LD in Case B-2.

Test results
Since the model LD was under a plane strain state, the sensors, which were set perpendicular to the side of the model LD (e.g., T1 and W1), measured the values of the same point. Consequently, the relation between the suction and water content of a measuring point can be obtained from the values of the tensiometer and water content sensor at that point. According to the test results, the model LD failed by overtopping in Case A-1, Case A-2 and Case B-1, while the failure mode of the model LD in Case B-2 was sliding failure. Moreover, the development processes of the seepage in the dams were different among the different cases. Unfortunately, some sensor data were missing due to technical issues. The available results of each case are introduced in detail as follows.
The results of Case A-1 are shown in Fig. 7 and Photograph 6. As shown in Fig. 7(b), with the increase in water level, the suction at all measuring points remained constant at rst. Then, the suction at T1 decreased signi cantly at T= 29 min, while the suction at other measuring points began to decrease at almost the same time (T=32 min). Meanwhile, it was also observed in Fig. 7(c) that the increase in water content at W1 was earlier than that at other measuring points. The reason for this phenomenon is that the model LD began to breach at T=31 min due to overtopping; the failure process is shown in Photograph 6. The soil at T1-W1 was saturated by the seepage ow from the upstream before overtopping, while the soils at other measuring points were saturated by water from the crest after overtopping. Since the in ltration path was signi cantly shortened, the soils at T2-W2, T3-W3, T4-W4 and T5-W5 almost became saturated at the same time. Thus, with a WLRR of 10 mm/min, the phreatic line developed to T1-W1 before overtopping, and the whole dam saturated quickly after overtopping.
The results of Case A-2 are shown in Fig. 8 and Photograph 7. As shown in Figs. 8(a) and (e), with the increase in the water level, the excess pore water pressure (EPWP) at P1 and P2 increased gradually, and their increments were approximately the same. This phenomenon indicates that sand can be easily saturated because of its high permeability, and the phreatic line in the upstream slope of the model LD was approximately horizontal. In addition, Figs. 8(b) and (c) show that the suction at T1-W1 started to decrease at T=118 min and that the phreatic line developed to T1-W1 at T=130 min. With the saturation of the upper layers in the dam, the vertical displacement at DS1 increased gradually, as shown in Fig. 8(f). Since the vertical displacement at DS2 was mainly contributed by dam breaching, the vertical displacement at DS2 was greater than that at DS1 at the end of the test. Dam breaching was triggered by overtopping at T=151 min, and the failure process is shown in Photograph 7.
For Case B-1, the WLRR was similar to that of Case A-2, as shown in Fig. 9(a). However, compared with that in Case A-2, the model LD in Case B-1 had a thinner brick layer and thicker soil layers, indicating that the soil of the well-graded mixture layer exhibited different initial states. Compared to the other layers, the well-graded mixture layer was closer to the initial water level, so the soil at T3-W3 and T5-W5 had lower initial suctions and higher initial water contents, as shown in Figs. 9(b) and (c). In addition, with the increase in water level, the suction at both T3-W3 and T5-W5 decreased gradually, and the water content increased gradually, which means that the phreatic line was approaching these two measuring points. As a result, the phreatic line in Case B-1 was able to develop downstream further than those in the other cases. Fig. 9(f) also shows that the vertical displacement of DS1 at the end of Case B-1 was greater than that of Case A-2. Therefore, the increase in the thickness of the soil layers might have led to greater vertical displacement at the dam crest during the increase in the water level. The decrease in the height of the dam crest could lead to earlier overtopping failure. In other words, the increase in the thickness of soil layers reduced the stability of the LD, a comprehensive result. The failure process of the model LD in Case B-1 is shown in Photograph 8.
In Case B-2, unlike in other cases, the water level decreased gradually to 337 mm after reaching the dam crest, as shown in Fig. 10(a). Fig. 10(b) shows that the suction at T3 and T5 decreased during the whole test, and the suction at T2 and T4 also began to decrease after T= 210 min. Moreover, the suction at all the measuring points reached 0 kPa at the end of the test. Hence, although the permeability of the gapgraded mixture was the lowest, the gap-graded mixture layer could also become saturated under seepage conditions. Interestingly, as shown in Fig. 10(d), although the s-q relation curves of the gap-graded mixture and the well-graded mixture are different, they are similar to the water retention curves obtained from water retention tests on corresponding materials (Ma et al., 2020). Finally, Fig. 10(e) shows that the vertical displacement at DS2 suddenly decreased at T= 1024 min, indicating the occurrence of a landslide. The landslide mass was composed of only the gap-graded mixture, as shown in Photograph 9.
According to the measurement data and the observed phenomenon, the development of the phreatic line in the model LD was estimated, as shown in Fig. 11. The ume test results indicate that the stability of the LD under seepage loading is signi cantly in uenced by the different debris materials in its heterogeneous strata. To obtain more detailed information, numerical simulations of the ume tests were conducted utilizing a suitable saturated/unsaturated constitutive model. Notably, according to the s-q relation for the model LD, as shown in Fig. 9(d) and 10(d), the water content could increase under constant suction conditions, which means that a deformation-dependent WRC model is necessary for accurate numerical simulation.

Saturated/unsaturated constitutive model and WRC model
The results of triaxial tests on three typical LD materials (Ma et al., 2020) indicate that the LD materials have a highly developed soil structure, as shown in Fig. 13. To accurately describe the behavior of LD materials, an unsaturated soil constitutive model (Xiong, 2020) was selected. To consider the structure of the soils formed in the natural deposit process, Xiong (2020) incorporated the concept of superloading (Asaoka et al., 2000) into a present constitutive model proposed by Zhang and Ikariya (2011), which utilizes skeleton stress and degree of saturation as independent variables. There are some studies on the application of the constitutive model (Zhang and Ikariya, 2011) in calculating the stabilities of an unsaturated slope and LD, and reasonable results have been observed-e.g., Xiong et al. (2014) and Xiong et al. (2018).
In addition, to properly consider the in uence of deformation on the degree of saturation, a WRC model proposed by Xiong (2020) was also utilized. Based on the experimental evidence, in this model, it is assumed that if subjected to only volumetric deformation, the shape of the intrinsic WRC at any void ratio does not change, but the curve shifts parallel to the S r axis in S r -s space. Coupled with the constitutive model, the WRC model has been used to describe the hydraulic/mechanical behavior of completely decomposed granite, and satisfactory accuracy has been observed (Xiong, 2020). Thus, it is reasonable to select this WRC model to simulate the development of the phreatic line in the model LD.

Parameter determination
In this research, to ensure the accuracy of the boundary value problem calculations, all the material parameters involved in the constitutive model of the saturated/unsaturated soil are determined via the element tests. The parameters involved in constructing the WRCs for the three LD materials were determined from the water retention tests (Xiong, 2020;Ma et al., 2020), and their values are listed in Table 3. The material parameters were essentially determined by the saturated triaxial test results (Ma et al., 2020) considering the available physical quantities obtained by the water retention tests. The calibrated material parameters of the three LD materials are listed in Table 4.  Tables 3 and 4. Generally, the calculated results agree well with the experimental data. Thus, the WRC model coupled with the unsaturated soil constitutive model (Xiong, 2020) can describe both the skeleton curves and scanning curves of LD materials with satisfactory accuracy.  Tables 3 and 4. Good agreement between the calculated and test results is observed, although some discrepancies exist. Throughout the shearing stage, the stress difference of the calculated and test results of the sand increases, while the stress differences of the well-graded mixture and gap-graded mixture rst increased and then gradually decrease. Moreover, the calculated EPWP is able to capture the overall trend of the tested EPWP. Therefore, numerical method used in this paper is able to properly describe the strain-stress-dilatancy relation of LD material in triaxial tests.
In the calculations, the saturated permeability of the LD materials was set according to the results of the permeameter tests conducted by Shi et al. (2015b). The unsaturated permeability of the LD materials was set by the formula of van Genuchten (1980).
Considering the high strength of the brick, its stress-strain relation was described by an elastic model. The material parameters of the brick in the numerical calculations were set as follows: Young's modulus of 30 MPa, Poisson's ratio of 0.25, speci c gravity of 2.8 and permeability of 8.54×10 -6 m/s.

Numerical model and initial conditions
The ume tests of the model LD were simulated with the soil-water-air coupled FEM program SOFT (Xiong et al., 2014), using the nite element-nite difference (FE-FD) scheme for soil-water-air three-phase coupling problems. There are ve cases in the numerical calculation, as shown in Table 5. In Case A-1, Case A-2, Case B-1 and Case B-2, the values of the soil parameters were the same as those listed in Tables 3 and 4. To investigate the in uence of the soil deformation on the development of the phreatic line, numerical calculation Case A-0 was conducted, which had the same WLRR and strata as Case A-2 but different soil parameters. In Case A-0, the parameter c e of the WRC model that controls the in uence of deformation on the degree of saturation was set to 0.0 in all the LD materials.
The nite element mesh used in the simulation was built according to Fig. 4, of which the size is the same as that of the model LD under plane strain conditions, as shown in Fig. 14. The nite element mesh is composed of 1351 nodes and 1255 4-node isoparametric elements. Fig. 15 shows the boundary conditions in the settings of the numerical model. Fig. 15(a) shows that drainage and vented boundary conditions are exactly the same as those for the ume tests. For the displacement boundary condition, as shown in Fig. 15(b), the bottom surface is xed in the x and y directions, and the other surfaces are free in both directions.
In addition to the boundary conditions, the setting of the initial conditions is also crucial in numerical simulation and will greatly affect the accuracy of the calculation. Since the model LD was prepared via the layered compaction method, the materials were heavily compacted at the initial condition. Mean effective stress of 12 kPa was added to the whole area to consider the compaction effect by comparing the simulated and test results. Fig. 16 shows the initial mean effective stress eld of the calculations.
For the suction and degree of saturation, the initial valves were essentially determined using the data measured during the ume tests and the WRC curves of the LD materials (Xiong et al., 2018;Ma et al., 2020), as shown in Table 5. Notably, the initial suction in the numerical calculations, 4 kPa or 6 kPa, is lower than the data measured during the ume tests. Xiong et al. (2018) and Ma et al. (2020) conducted water retention tests on LD materials with maximum grain sizes of 40 mm and 2 mm, respectively. The test results show that LD materials with a maximum grain size of 40 mm have a much lower suction at the residual degree of saturation state. Hence, the conclusion can be drawn that particles with grain sizes greater than 2 mm also in uence the water retention properties of LD materials. Considering that the length of a tensiometer probe is only 15 mm, the suction measured during the ume tests could be the local suction of small particles in the LD materials. Therefore, it is reasonable to set the initial suction in the numerical calculations as shown in Table 5.
A water head of 50 mm was loaded on the upstream slope of the model LD in the initial stage of the numerical simulations. After 1 h, the physical quantities at this moment were regarded as the initial state of the model LD (T=0 min). Subsequently, a prescribed increment of water head is applied to simulate the variation in water level. The WLRR of Case A-0 was set to be the same as that of Case A-2. Except for Case B-2, for all the other cases, after the water head peaked, a 420 mm high water head was quickly applied to the crest and the downstream slope of the model LD to simulate the overtopping process. At this time, in Case B-2, the upstream water level was unloaded to 337 mm over the following 600 min.

Flume test simulations
To demonstrate the performance of the numerical method used in this research, the calculation results for Case A-0 and Case A-2 are rst compared with the experimental results. Fig. 17 and Fig. 18 show the distributions of the degree of saturation at different times in Case A-0 and Case A-2, respectively. When considering the in uence of the deformation on the WRC, the phreatic line can clearly develop further downstream. Compared with the tested results shown in Fig. 11(b), the calculation result for Case A-2 coincides better with the tested development of the phreatic line. Therefore, the numerical method used in this paper is proved to have an accuracy that can adequately describe the development of phreatic lines in ume tests. Fig. 19 shows the test and simulated results for the mechanical quantities of Case A-2. Fig. 19(a) shows that the rising rates of the EPWP at P1 and P2 are approximately the same, and the maximum EPWP at P2 is greater than that at P1; these results are in good agreement with the test results. As shown in Fig.  19(b), the magnitude and the trend of the vertical displacement obtained from the test and the simulation are essentially the same. Although the breaching process of the LD cannot be completely simulated using the FEM, the simulated vertical displacement at DS2 decreased suddenly at T=149 min, which can indicate the occurrence of dam breaching. Moreover, once the overtopping breaching of an LD occurs, the entire breaching process will continue and cannot stop (Yang et al., 2015). The time of the occurrence of dam breaching is called the longevity of the LD in this research. Therefore, the numerical method proposed in this research can be utilized to estimate the occurrence of dam breaching and the longevity of the LD. Fig. 20 shows the distribution of the degree of saturation at different times in Case A-1. Compared with Case A-2, although the WLRR tripled in Case A-1, the phreatic lines at the same water level almost coincided. In addition, the maximum vertical displacement at DS1 in Case A-1 was 11.75 mm and slightly smaller than that in Case A-2, as shown in Fig. 19(b) and 21(b). The longevity of the model LD in Case A-1 was 34 min, which was one-third shorter than that in Case A-2. Therefore, although the WLRR has little in uence on the development of the phreatic line in the model LD, its increase can reduce the stability of the model LD.
For Case B-1, the phreatic line developed much further downstream, as shown in Fig. 22. Since the model LD in Case B-1 had thicker soil layers and a lower initial suction than those in Cases A-1 and A-2, the LD materials were more easily saturated by seepage. Fig. 23(a) shows that the EPWP at P3 increased to be greater than that at P2 in both the tested and simulated results. In addition, according to the simulated vertical displacement at DS2 shown in Fig. 23(b), the longevity of the model LD in Case B-1 was 156 min, which was longer than that of the model LD with the same WLRR in Case A-2. The reason for this phenomenon is considered that after overtopping, it took more time to completely saturate a model LD with thicker soil layers.
For Case B-2, the phreatic line developed to the downstream slope surface of the model LD, and the wellgraded mixture layer was totally saturated at the end of the test, as shown in Fig. 24. Fig. 25 shows good agreement between the calculated and test results in terms of the changes in EPWP and vertical displacement. Although it is di cult to estimate the occurrence of a landslide by the simulated vertical displacement at DS2, the displacement vectors of the downstream slope in Case B-2 are much greater than those in other cases, as shown in Fig. 26, which can also be used as a basis for determining the occurrence of a landslide. Moreover, Fig. 26(d) shows that the displacement is mainly concentrated in the gap-graded mixture layer. According to the triaxial test results of the LD materials (Ma et al., 2020), the strength of the gap-graded mixture is the lowest, and it is most intolerable to the seepage loading, which can explain why the landslide mass was composed only of the gap-graded mixture, as shown in Photograph 9.
A comparison between the test and the calculated results of all the test cases shows that based on a rational constitutive model and a deformation-dependent WRC model, the numerical method proposed in this study offers satisfactory accuracy for describing the stability of an LD under seepage loading conditions. Moreover, the hydraulic/mechanical behaviors of the LD materials and the heterogeneous strata of the LD play an important role in its stability. The height of the dam crest decreased with increasing water level due to the saturation of the LD materials. Hence, increasing the thickness of the soil layers can reduce the stability of an LD. In addition, owing to the low strength of the gap-graded mixture, sliding failure of the Tangjiashan LD could occur in the gap-graded mixture layer when the LD lake had a high water level. However, due to the small height-to-length ratio of the Tangjiashan LD, the phreatic line in the dam body is di cult to develop on the surface of the downstream slope before overtopping breaching. As a result, regardless of whether the spillway is constructed, with a continuously increasing upstream water level, even though its rising rate is small, Tangjiashan LD would fail by overtopping after formation. Therefore, the construction of a spillway on the Tangjiashan LD only caused the earlier occurrence of overtopping and did not increase the stability of the dam body. Considering the different debris materials in the heterogeneous strata of the Tangjiashan LD, it is essential to increase the stability of the dam body by controlling the water level of the LD lake. For example, excavation of a drainage tunnel through the Hongshiyan LD successfully reduced the dam breaching risk (Shi et al., 2017), and this LD was then rebuilt into a large-integrated water conservancy complex . If engineering measures would have been taken to control the water level below the gap-graded mixture layer, the Tangjiashan LD could have maintained stability, and the water resources of the LD lake could have been reasonably controlled and su ciently utilized.

Conclusions
In this paper, taking the Tangjiashan LD as a reference example, a systematic ume test program was conducted with three typical LD materials subjected to rising water levels. Particular attention was paid to the in uence of different debris materials in the heterogeneous strata on the LD stability. Meanwhile, soilwater-air coupled nite element analyses, based on a constitutive model for saturated/unsaturated soils and a nite deformation-dependent WRC model (Xiong, 2020), were conducted to simulate the corresponding ume tests as boundary value problems. The following conclusions were obtained.
1. From the geological and hydrological information of the Tangjiashan LD, three important characteristics were summarized, suggesting that proper selection of the debris materials for the model LD and careful zoning of the strata are necessary for clarifying the hydraulic/mechanical behaviors of LDs subjected to rising water levels.
2. The ume tests of the model LDs showed that the stability of an LD under seepage loading conditions is signi cantly in uenced by the different debris materials in its heterogeneous strata.
The model LDs failed by overtopping in Case A-1, Case A-2 and Case B-1, while the failure mode of the model LD in Case B-2 was sliding failure, and the landslide mass was limited to the layer of gapgraded mixture. From the measured data and observed phenomena, the development of the phreatic line in the model LD was successfully identi ed, e.g., it developed much further downstream in Case B-1 and Case B-2 than in the other cases.
3. A comparison of the results from the ume tests and the corresponding numerical calculations indicates that the proposed numerical method can describe the behavior of an LD subjected to a rising water level with satisfactory accuracy. It is clari ed for the rst time that the calculations that adopt the deformation-dependent WRC model can describe the observed phreatic line much better than those that do not adopt the deformation-dependent WRC model, showing the possibility of using the proposed numerical method to estimate the occurrence of dam breaching and the risk of LD failure. 4. As the water level increased, settlement of the dam crest occurred due to the saturation of the LD materials. Hence, the thicker the soil layer is, the lower the stability of the LD. In addition, owing to the low strength of the gap-graded mixture, sliding failure of the Tangjiashan LD might occur in the gapgraded mixture layer under seepage conditions. 5. Due to the small height-to-length ratio of the Tangjiashan LD, it is di cult for the phreatic line in the dam body to reach the surface of the downstream slope before overtopping breaching. As a result, with a continuous increase in the upstream water level, even though its increasing rate is small, the Tangjiashan LD will fail by overtopping without any treatment. The construction of a spillway on the Tangjiashan LD only caused the overtopping breaching to occur earlier but did not increase the stability of the dam body. Therefore, if it is possible to construct a diversion tunnel to avoid overtopping of LD, it might a wiser way not to destroy an LD whenever it is formed. From an engineering viewpoint, it is possible to utilize a naturally formed LD instead of destroying it. For example, the dam breaching risk can be reduced by excavation of a diversion tunnel on the condition that the dam stability is carefully estimated based on not only accurate hydraulic and geological data but also construction conditions.

Declarations
Tables       Figure 1 Longitudinal section of the Tangjiashan LD along the river Figure 2 Elevation of water level at the Tangjiashan LD lake (data from Zhicheng water gauging station (2008), close to the upstream of the Tangjiashan LD)

Figure 3
Grain size distribution curves of the LD materials used in the ume tests Side and top views of the tensiometers set in the model ground whose thickness is half the width of the ume and supported with a side plastic plate  Initial mean effective stress eld (unit: kPa) Figure 17 Calculated variation in saturation due to a rising water level for Case A-0

Figure 18
Calculated variation in saturation due to a rising water level in the Case A-2 ume test Page 29/33

Figure 19
Test and calculated results of the mechanical quantities in the Case A-2 ume test Simulated variation in saturation due to a rising water level in the Case A-1 ume test Page 30/33

Figure 21
Simulated results for the mechanical quantities in the Case A-1 ume test Figure 22 Simulated variation in saturation due to a rising water level in the Case B-1 ume test Page 31/33

Figure 23
Test and simulated results of the mechanical quantities in the Case B-1 ume test Figure 24 Simulated variation in saturation due to a rising water level in the Case B-2 ume test

Figure 25
Test and simulated results of the mechanical quantities in the Case B-2 ume test Figure 26 Distributions of the displacement vectors at the end of the ume tests (unit: m)