Numerical Investigation of the Landslide and Its Surge: A Case Study of the Gongjiafang Landslide in the Three Gorges Reservoir Area

Once the slope in the near bank area enters the state of failure, a geological disaster chain caused by landslide and its generated surge occurs very likely. In this study, a two-layer depth-averaged model was used to describe the disaster chain. The HLL (Harten-Lax-van Leer) finite volume method was used for numerical simulation and analysis. Meanwhile, the linear interpolation technique was employed to obtain second-order accuracy. The numerical results of the analytical examples reflect the movement characteristics of the two-layer fluid and verify the correctness of the numerical model. On the basis of numerical verification, the Gongjiafang landslide and its surge were simulated and analyzed. In the early stage, the deformation of the underwater geomaterials disturbed the water, forming the prototype of the surge, and then the landslide movement promoted the rapid development of the surge. After the landslide was deposited, the surge continued to travel forward and formed the largest form near the opposite bank. The numerical simulation is applicable to complex terrain and reveals some mechanisms and characteristics of the disaster chain. Compared with empirical methods, the numerical model adopted could reproduce the process of disasters chain more accurately and effectively and then improve the understanding of the disaster chain. It is feasible that the proposed numerical model can be applied under approximate plane strain conditions but is no longer applicable under 3D conditions. This work can provide reference for further research on disaster chain caused by landslide.


Introduction
The landslide and generated surge caused by earthquakes, rainfall, and reservoir filling are devastating geological disasters. A large landslide in the Vajont reservoir in Italy in 1963 is a representative case of this kind. A dam break and about 2,000 casualties were caused by the landslide and its surge [1,2]. These disasters threaten the navigation of vessels, the stability of dams, and the life and property safety along the banks [3,4]. This kind of disasters with great destructiveness has aroused great concern of scholars and has become a hot issue in the field of geological disaster research in recent years.
To deeply study the mechanism and characteristics of landslide and its surge, a lot of research have been carried out in the field observation [5], model experiment [6], theoretical analysis [7], empirical method [8], and numerical simulation [9][10][11]. Among them, the field observation pro-vides the most reliable data, but due to the unpredictability of disasters and the limitations of field data, it cannot provide sufficient support for research. Experiments can be designed for specific research goals, and the observation and experimental data provided by them are important references for research but limited by experimental conditions and testing methods; so, the experiment methods cannot meet all the needs of research. Theoretical analysis is the foundation of research, but the analytical analysis cannot be applied to complex situations. Empirical methods are simple and easy to apply, but their applicability is often limited. The numerical simulation can realize the back analysis and prediction of the disaster process under complex conditions. Through numerical simulation, the mode of slope failure [12,13] and the effect of protective structure [14] can be studied. The numerical simulation can further solve complex problems such as energy transfer in disasters [9] and the impact of surge on protective structures [10]. The discussion of "drag along," "lift up," and "push ahead" effects in the disaster chain has deepened the understanding of the interaction between landslide and water [11]. With the development of modern computing science, numerical simulation has become one of the indispensable means to study the geological disaster. Nonetheless, it must be noted that the numerical investigation should be verified by field observations, experimental data, and analytical solutions.
In the existing numerical simulation of landslide and its surge, the movement model of landslide is mainly divided into two types: rigid body model and deformed body model [15]. For the rigid body model, the landslide is considered as a rigid body or a bottom boundary varying with time, and the actual deformation of the landslide during motion process is generally ignored. For the rigid model, the kinematics and semiempirical formulas need to be combined to define the motion of the rigid body, with the consideration of the effects of gravity, buoyancy, friction, flow resistance and inertial force, etc. [16,17]. The rigid body models are generally suitable for the landslides with higher stiffness. Numerical and experimental studies have shown that the deformation of landslide has a great impact on the landslide and its surge [18,19]. Therefore, to make a more accurate calculation of landslide and its surge, the effects of deformation and rheological properties should be considered in the numerical simulation. The landslide can be regarded as a deformed body, and a fluid is commonly used to describe the landslide. For the deformed body model, the landslide and the surrounding water can be regarded as twophase fluid, which is distinguished by the volume fraction [20]. It is also possible to treat the landslide body in the lower layer and the water body in the upper layer as two separate layers, each layer having its own rheological properties. It is worth mentioning that due to the complexity of material composition and motion, the landslide generally adopts the depthaveraged models. The main difference of these models is that different rheological models are used [21][22][23][24].
The models currently used to simulate the surge process are mainly the depth-averaged models and the models based on N-S (Navier-Stokes) equations. The depth-averaged models are mainly represented by shallow water equation and Boussinesq equation, which has three main forms, namely, typical Boussinesq equation [25], extended Boussinesq equation [26][27][28][29], and multilayer Boussinesq equation [30]. The shallow water equation is currently the main choice for most scholars to simulate the surge process because of its strong applicability and high computational efficiency [25]. Boussinesq equation can obtain higher order precision, mainly focusing on the nonlinearity and dispersivity of surge wave [29], but the calculation cost is relatively high. For more complex surge process, it is necessary to solve the N-S equations [31]. Since the N-S equations require large storage, it also needs to cooperate with the corresponding interface tracking technology. As these solution methods of N − S equations are computationally inefficient, it is difficult to solve the N-S equations for large-scale engineering. The landslide and its surge models are coupled by phase-to-phase coordination and forces. The meshing methods such as finite difference method [32], finite element method [33], and finite volume method [15], as well as the meshless methods such as smoothed particle hydrodynamics method [34], moving particle semi-implicit method [35], and discrete element method [31,36], are used for the numerical solution of the model of landslide and its surge. These numerical methods have their own advantages and are widely used in the numerical simulation of landslide and its surge.
The goal of this study is to achieve an effective simulation of the landslide and its surge, so as to analyze the mechanism and characteristics of the disaster chain combining with typical case. For the purpose of this study, a numerical investigation of the Gongjiafang landslide and its surge in the Three Gorges region, China, was carried out using the depthaveraged theoretical model and the finite volume method. The main contribution of this work is to fully demonstrate the feasibility and reliability of the proposed numerical scheme for simulating the disaster chain, thus providing a good numerical scheme option for the numerical investigations of similar disasters. At the same time, some understanding of the disaster chain is also helpful for future research.
In the study of this case, the applicability of numerical and empirical methods was analyzed through comparison, so as to provide reference for further research of similar disasters.

Numerical Model
Since the density of the landslide is relatively large, it moves under the water, and the following subscript m = 1, 2 indicates the variables of the water and the landslide, respectively, as shown in Figure 1.
The mass and momentum equations of incompressible medium can be expressed as Equation (1): where t represents the time, ρ is the density of the medium, u is the velocity vector, and T is the Cauchy stress tensor. The water and the landslide are, respectively, subjected to averaging processing in the depth direction, and the variables are defined as averages along the depth, which can be expressed by Equation (2): where h m is the depth of flow, u i,m is the flow velocity, σ ij,m is the stress component, and u i,m and σ ij,m represent the velocity and the stress component after average along the depth, respectively. By depth-averaging the mass and momentum equations of incompressible continuous medium, the model equation used to simulate the landslide and its surge can be expressed by Equation (3) where g is the gravitational acceleration, k is the lateral stress coefficient, z b2 and z b1 represent the elevation of the sliding surface and the two-phase interface, respectively, r = ρ 1 /ρ 2 is the density ratio, and τ xz is the shear stress.
The shear stress τ xz (z b1 ) between the two layers in Equation (3) can be expressed as Equation (4): where κ is a positive constant. The bottom shear stress τ xz (z b2 ) of the landslide is affected by the pore pressure and can be expressed as Equation (5): where φ b is the bottom friction angle, and μ b is the underwater friction coefficient. In this study, Equation (3) was discretized by the finite volume method of the HLL scheme (see the appendix for the scheme). The HLL scheme used in this study to calculate the interface flux is relatively simple in solving the Riemann approximation problem. In order to improve the accuracy of the numerical solution, the second-order linear interpolation technique is needed to reconstruct the variables on the left and right sides of the interface. The code was originally developed based on the model theory and numerical methods described above and used in the numerical investigation of landslide-induced disaster chain in this study. The appendix provides a detailed introduction to the numerical scheme.

Numerical Results and Analysis
3.1. Dam-Break Flow. In order to verify the effectiveness of the numerical scheme in calculating the discontinuous transient flow, the numerical simulation of dam-break flow was carried out. The initial condition of the two-layer dambreak flow is shown in Figure 2.
In order to facilitate the comparison with the analytical solution, the same density was applied to the two-layer fluid in the numerical simulation (this means that the density ratio r is equal to 1), and the two cases shown in Table 1 were calculated and compared. The dam break is located at x = 10 m. Under the condition that the upstream and downstream water depths are determined, the analytical solution can obtain a free surface at a certain moment [37]. In order to be consistent with the assumptions of the analytical solution, the numerical calculation does not consider the resistances, and the lateral stress coefficient k is set to 1.
The numerical results of case 1 and case 2 at 0.8 s are shown in Figure 3.
In the calculations of case 1 and case 2, since the total water depths in the upstream and downstream are identical (i.e., h L = h R ), the free surface of the calculated dam-break wave should be theoretically the same. It can be seen from the results of the dam-break flow that the transient free surfaces of the two cases both have good consistency with the analytical solution. This difference does not seem obvious, indicating that the numerical scheme used is correct and effective. At the same time, since the upper and lower layers have the same density in the two cases, a consistent contact discontinuity obtained after exchanging the initial depths of the upper and lower layers can confirm that the calculation is reasonable. The results show that the contact discontinuity obtained by the two cases has almost the same profile. In despite of some differences from the analytical solution, they are generally consistent with the analytical solution and show good numerical stability, without numerical oscillation. Due to the limitation of space, only the two-layer dam-break flow with the same density was simulated in this study, while the dam-break flow with different densities in the upper and lower layers have not been fully discussed. This is because the numerical solution under the condition that the two layers have the same density can facilitate the comparison with the analytical solution and then verify the correctness. Obviously, in the case of two layers with different densities, the advantages of the two-layer fluid model can be fully exploited. To deeply apply the two-layer fluid model, it is necessary to further study the engineering cases under the condition that the upper and lower layers have different physical properties. Nevertheless, the numerical simulation of the dam-break flow show that the two-layer flow can be transformed into a single-layer flow with the same density, and the calculation can capture the flow characteristics of the two-layer dam-break flow, which can be further applied to the study of the landslide and its surge. The mountain where the landslide is located is steep, and the relative height is about 600 m. The slope is generally eastsouth, with an original average incline of 53°. The landslide is bounded by two gullies, and the rear edge elevation is about 400 m. The rock of the landslide has a fragmented layer structure. The landslide occurred on November 23, 2008. On that day, the Three Gorges Reservoir was in the storage stage. In the Gongjiafang section, the water level was 172 m, the width of the Yangtze River was about 480 m, and the minimum elevation of the riverbed was about 40 m. After the landslide occurred, the investigation revealed that the exposed fresh rock surface appeared as an isosceles trapezoid with an upper width of 45 m and a width of 194 m at the water surface [8]. The longitudinal length of the landslide was 293 m, the height was 210 m, the upper slope was 64°, and the lower slope was 44°. The water entry angle of the landslide was about 40°. The comparative analysis of the topography before and after the landslide shows that the average thickness was 15 m, and the volume was 380,000m 3 . According to the image at the time of the slope failure and the underwater measurement analysis made by the Chongqing Three Gorges Geological Disaster Office, the depth of the landslide after entry into the water was about 50 m. According to the field investigation, the surge caused the cables of several ships docked to be broken and the bottom of the ships to be damaged [38]. Fortunately, no ships have been overturned, and no casualties have been caused. The regional location map and site pictures of the Gongjiafang landslide are shown in Figure 4, and the main information of the Gongjiafang landslide is shown in Table 2.
The Gongjiafang landslide occurred when the water level of the Three Gorges Reservoir reached 172 m. For the internal induced factors related to the geological structure, the mountain where the Gongjiafang landslide is located was cut by two structural planes, and the rock on the slope surface was broken as a fragmented and scattered structure [8]. These are the internal controlling factors of the collapse. At the same time, the left and right parts of the landslide are topographically bounded by the gullies without lateral constraints, and the slide along joint surface was caused by the level uplift of groundwater. As for the external cause, the water level of the Three Gorges Reservoir rose, the front edge of the landslide was hollowed out with the entrainment, and then the support at the slope foot was lost. At the same time, the groundwater level in the landslide rose rapidly, resulting in the continuous decrease of the strength of the slip zone, and the landslide eventually occurred under the influence of internal and external factors. The surge caused by the rapid impact on the water after the slope failure is an important secondary disaster. As the spreading scope of the surges continued to expand, sometimes, the consequences of disasters are much more serious than the landslides themselves. Therefore, it is necessary to conduct a more in-depth study on the generation and propagation of surge in combination with landslide movement.
According to the actual situation of the disaster, the range of the horizontal distance was 0-1200 m, and the elevation was 0-450 m in the numerical simulation. In this model, the elevation of the rear edge of the landslide was 403 m, the maximum elevation of the opposite bank of the landslide was 397 m, and the initial water level of the river was 172 m, as shown in Figure 5.
The calculation parameters must be calibrated first in order to better represent the propagation process of the Gongjiafang landslide and its surge. In consideration of unavailability of enough information, final deposition of the landslide, and the maximum surge height, the two pieces of the most important known information about disaster were used for back analysis of the calculated parameters. According to the existing data, the density ratio r was 0.4 (the specific gravity of water and geomaterials was 1 and 2.65, respectively) in calculation. Considering the movement characteristics of Gongjiafang landslide fluidization, the lateral stress coefficient k was 1, so that the underwater friction coefficient μ b became the most important parameter for back calculation. After back calculation and analysis, the final result obtained at μ b = 0:42 complies with the actualities. See Figure 6 for the landslide deposition and the maximum surge height obtained by the numerical calculation.
It can be concluded from Figure 5 that the deposition obtained by numerical simulation is basically consistent with actualities under the given calculation parameters. Since the deposition of actual landslide expands along the river in the deceleration zone, the deposition calculated on the vertical section is larger than that in actualities, which also leads to certain difference between the numerical results of the front and rear edges of the depositions and those in actualities, even if the overall deposition obtained by the numerical simulation is consistent with that in actualities. All these differences also exert certain influence on the formation and spread of surge (mainly on the front edge of landslide). To be specific, the maximum surge height calculated is 32.65 m (the result of t = 15:1 s), slightly higher than the actual maximum surge height 31.8 m (the front edge of the landslide obtained by numerical calculation moves longer, which pushes water more   Case 4 Geofluids significantly). The differences can be ignored for the calculation of a large scale; so, the differences between the numerical result and the actual maximum surge height are believed limit and acceptable. Based on the comparison between numerical calculation and actual disaster, the numerical simulation is indeed correct and effective.
The characteristics of the Gongjiafang landslide and its surge were further studied through numerical simulation on the basis of verifying the correctness and effectiveness, as shown in Figures 7 and 8.
The occurrence and development process of the disaster chain caused by the landslide are shown in Figure 6, where   5 Geofluids some characteristics of the process are also fully demonstrated. The overall forward movement of landslide is insignificant when the slope becomes unstable but a few waves form on the water surface, which is caused by the deformation of unstable geomaterials. At this time, no obvious displacement of landslide front edge can be found. This indicates that the formation of the surge is closely related to the deformation of the landslide; although, the landslide movement may be the more dominant cause of the formation of the surge (Figure 7(a)). Subsequently, the unstable geomaterials moved downward along the sliding surface under gravity, and the front edge of the landslide moves with a faster speed. The landslide pushes the water to move towards the opposite bank, and several waves spread to the opposite bank under the pushing effect of the geomaterials. During this period, the movement of the front and rear edges of the landslide is similar, and the movement is mainly in the form of a whole (Figure 7(b)). The landslide moves slower after the front edge of the landslide enters the bottom of the river bed, and the waves caused by the landslide continue to spread forward and will enter the shallow water area. The shape of the landslide changed significantly during this period, and it became longer due to the much slower movement of the rear edge. (Figure 7(c)). Wave height increases along with the spreading of the surge to shallow water area. The difference in the veloc-ities of several waves causes the spacing between the wave crests to gradually decrease. Although the front edge of the landslide moves forward at a slow speed, its overall impact on the surge has been weakened (Figure 7(d)). The front edge of the landslide reaches its final position along with time passing by, but the geomaterials have not yet reached the balance state, and the main movement occurs at the rear edge of the landslide. Meanwhile, several waves superimpose in the shallow water area, which boosts the formation of a higher surge gradually. At this time, the impact of the landslide on the surge is negligible (Figure 7(e)), when the surge will start run-up movement upon its arrival at the opposite bank. Meanwhile, the landslide almost does not move any longer, and the deposition adjusts its final shape mainly by deformation (Figure 7(f)). Figure 8 shows the results of several critical moments on the same figure, which clearly shows the dynamic process of maximum surge formation. The results of surge runup under vertical section calculation are inaccurate (in the numerical model, the resistance of the surge in the run-up process is much greater than the resistance along the river. Therefore, the surge spreads to upstream and downstream after reaching the opposite bank but its run-up to the opposite bank is extremely limited). Therefore, the run-up and the possible further back spread process of surge were not simulated in the study. Nevertheless, it could be seen from the above analysis  There are indeed some studies based on landslide and its surge, which have played a huge role in promoting the study of disaster chains. However, the author did not use experiments to verify the numerical model because the experimental data and information we can obtain are not enough, or some experiments are not suitable for the numerical model of this study. The model equations in this study can simulate landslides with fluidized motion, which is significantly different from the disaster chain caused by landslides that move in rigid body form as a whole. The HLL finite volume scheme was used to numerically discretize the model equation, and this numerical method has good conservation. Since the linear reconstruction technique was also used, the second-order accuracy can be achieved, which is sufficient for large-scale disasters. The numerical scheme used also has the advantage that the method can handle discontinuities as well as wet and dry interfaces. Through case analysis, it is confirmed that the adopted method is effective and feasible, thus providing an applicable choice for the study of landslide and its surge.
To verify validity of the existing empirical formula, the results obtained by the empirical formula were further   Geofluids compared with that in actualities and analyzed. In the disaster chain caused by the landslide, the maximum surge height plays a key role in evaluating the scale of the disaster chain, which may be affected by factors such as geomaterials volume, landslide speed, and water depth. At present, some formulas for calculating the maximum surge height have been proposed [8,40,41], as shown in Table 3: Using the above empirical formulas to calculate the surge caused by the Gongjiafang landslide, in addition to the coefficients contained in the formulas, the calculation parameters of the Gongjiafang landslide and water used in the formula calculation are the same as the numerical simulation. The calculation results of the maximum surge height caused by the Gongjiafang landslide using the above empir-ical formulas are shown in Table 3. The results of the empirical formulas in Table 4 are provided by the literature [8].
The calculation results show that the calculated values of the S&V method and the IWHR method are relatively small, and the calculation errors are the largest; so, these two methods are not suitable for the calculation of the surge caused by the Gongjiafang landslide. The formulas of the Noda method and Pan method have similar expressions, and the calculated results range from 17.90 m to 23.80 m, which are closer to the actual situation. Compared with the formula calculation method, the maximum surge height obtained by numerical simulation is the most consistent with the actual situation, which is 32.65 m. However, since the numerical simulation does not consider the surge propagation out of the plane l is the length of the landslide, s is the thickness of the landslide, b is the average width of the landslide, h is the water depth, ρ is the density of the landslide, ρ w is the density of water, u is the sliding velocity, and a is the surge height.
IWHR a max = ku 1:85 v 0:5 /2 g k is the comprehensive influence coefficient, which takes the average value 0.12; v is the volume of the landslide, and u is the velocity of the landslide when entering the water.
Noda a max = 1:32hu/ ffiffiffiffiffi ffi gh p a max is the maximum surge height, u is the sliding velocity of the landslide, and h is the water depth.

Pan
When the landslide enters the water mainly in horizontal motion: ξ 0 = 1:17hu/ ffiffiffiffiffi ffi gh p When the landslide enters the water mainly in vertical motion: ξ 0 = h•f u′/ ffiffiffiffiffi ffi gh p The ξ 0 function relationship can be expressed in sections as follows: (i) When 0 < u′/ ffiffiffiffiffi ffi gh p ≤ 0:5, ξ 0 = hu′/ ffiffiffiffiffi ffi gh p (ii) When 0:5 < u′/ ffiffiffiffiffi ffi gh p ≤ 2, f u′/ ffiffiffiffiffi ffi gh p changes in a curve (iii) When u′/ ffiffiffiffiffi ffi gh p > 2, ξ 0 = h ξ 0 is the surge height, u is the horizontal sliding velocity, u′ is the vertical sliding velocity, and h is the water depth. 9 Geofluids direction, the obtained result is larger than the field observation value. The results in Table 3 show that the actual surge process caused by landslides has certain complexity. Using empirical formulas to calculate the maximum surge height does not yield very reliable results, because these formulas are obtained under various theoretical assumptions and simplifications, as well as empirical summaries, with limited applicability. Compared with the formula method, numerical simulation has stronger applicability, and reliable results can still be obtained under complicated conditions.

Discussion
The analytical example used in this study is simple, which not only verifies the correctness and effectiveness of the numerical model but also confirms the ability of the method to handle discontinuous problems. Obviously, the two-layer model can be reduced to a single-layer model under certain conditions; so, it is an effective extension of the single-layer model.
During the numerical investigation of the case, some mechanisms were revealed. The front edge of the landslide was first deformed significantly, and the deformation disturbed the water, thus forming the prototype of the surge. When the landslide moved downward under gravity, the landslide pushed the water forward and promoted the development of the surge. At this time, the landslide movement became the main reason for the promotion of the surge. Affected by the irregular shape and complex movement process of the landslide, several waves with different frequencies were formed. After the landslide movement stopped, the waves continued to propagate forward, and several waves superimpose and form the largest surge in the shallow water area close to the opposite bank. In the case, the above mechanism is representative to a certain extent, but due to the extremely complex process of surge caused by landslide, the disaster chain mechanism in more complex situations still needs to be further studied.
In the investigated case, the numerical deposition is in good agreement with the actual situation, but because the numerical method cannot consider the lateral movement of the mass on the riverbed, the numerical deposition is larger than the actual situation. Although the lateral movement is only a minor factor in this case, ignoring it makes the landslide move further distances and also makes the numerical maximum surge larger than the actual one, which is not a significant difference. Based on the above analysis, the numerical model can make numerical deposition and maximum surge larger than actual ones. The calculation in this study does not discuss more about the maximum run-up but focuses on the maximum surge height. This is because when the surge approaches the opposite bank, it must not only make a run-up movement towards the bank but also spread to the upstream and downstream of the river.
The relatively high resistance to the run-up movement forces the surge to spread in the direction along the river when approaching the bank. As a result, the calculated run-up height is obviously larger, which is also a major shortcoming of the vertical section model. However, the maximum surge does not generally occur in the near bank area; so, the calculation of the maximum surge height is not much influenced by the numerical model, and the calculation of the maximum surge is more in line with the actual situation.
In addition, the model used in this study still cannot effectively simulate some special phenomena in the process of landslide and its surge, such as the process of mixing geomaterials with water under high-speed motion on top surface of the landslide and the " hydroplaning" [42] caused by water wedging into the bottom of the landslide (which greatly enhances the sportiness of landslides). Obviously, the study of these complex processes is our future work. However, in order to simulate these phenomena, it is necessary to add experiments and theoretical research on the basis of existing ones and build more complex and accurate 3D models to complete the further expansion of the theoretical and numerical models.
Based on the above discussion, the method in this study is mainly limited by the application environment. For complex disaster processes with obvious three-dimensional effects or special mechanisms, this method is no longer applicable. In future research, the method can be extended to three dimensions, and the equations and solutions for special mechanisms can be supplemented to better cope with the application limitations of the current work. Although the numerical model in this study cannot solve all the problems in the disaster chain, it has the advantages of simple numerical scheme and high efficiency, can still effectively simulate the disaster with plane strain, and quickly realize the back analysis and prediction. Therefore, it is an important and effective numerical model for studying the characteristics and development of disasters.

Conclusions
In this study, a two-layer depth-averaged continuous fluid model and a HLL finite volume method were used to simulate the landslide and generated surge. The main conclusions are as follows: (1) Numerical results of analytical and practical examples show that the adopted numerical model can reflect the deformation and deceleration of the lower layer due to the blockage of the upper layer and also reflect the fluctuation of the upper layer driven by the lower layer. The abovementioned movement characteristics of the two-layer fluid correspond to the landslide and its surge process; so, the effective

Appendix
The elevation of the sliding surface and the conservation variable U were both defined in the center of the grid. This treatment avoids the creation of semidry and semiwet element, and Green's theorem was applied to obtain the discrete equation as Equation (A.1): where Δt is the time step, Δx is the space step, F i+1/2 and F i−1/2 are the numerical fluxes at the interface of x i+1/2 and x i−1/2 , respectively, and S i is the source term. The HLL scheme is suitable when dealing with dry and wet elements. The scheme for flux can be expressed as Equation (A.2):  where ∂U/∂x is the intraelement variable gradient. The minmod limiter was used to avoid nonphysical numerical fluctuations, as shown in Equation (A.10) [43]:

Data Availability
The numerical calculation data used to support the findings of this study are available from the corresponding author upon request.