Abstract
A shallow-water equation (SWE) is used to simulate earthquake-induced water waves in this study. A finite-difference method is used to calculate the SWE. The model is verified against the models of Sato and of Demirel and Aydin with three kinds of seismic waves, and the numerical results of earthquake-induced water waves calculated using the proposed model are reasonable. It is also demonstrated that the proposed model is reliable. Finally, an empirical equation for the maximum water elevation of earthquake-induced water waves is developed based on the results obtained using the model, which is an improvement on former models.
1. Introduction
Strong earthquakes shake the banks of reservoirs and can consequently generate large water waves that may result in a geohazard in mountainous earthquake-prone areas [1]. Hydrodynamic pressures in dam reservoirs during earthquakes were first rigorously analyzed by Westergaard in 1933 [2]. In his study, the hydrodynamic equations were derived in terms of the theory of elasticity of solids, and an analytical solution was obtained. Later, Sato modeled a seismic wave as a simple sine wave and proposed a formula for an earthquake-induced water wave using analytical calculations based on linearized shallow-water equations (SWEs) [3]. Chopra and coworkers provided a revised solution to account for the effects of water compressibility on pressure distribution and to investigate dam-water interaction using a finite-element method [4, 5]. They found that the effects of water compressibility for most dams should be considered in their earthquake-response analysis. Hung and Chen analyzed the variation of the water surface and the hydrodynamic pressures of dam-reservoir systems during the El Centro Earthquake, solving Euler equations by a finite-element method [6]. Chen et al., solving three-dimensional coupling Euler equations and incompressible continuity equations with a complete three-dimensional (3D) finite-difference scheme, calculated the hydrodynamic pressures and water-level increase for nonlinear conditions, including the nonvertical dam face and the nonhorizontal reservoir bottom. They concluded that the geometry effects of reservoir shape on hydrodynamic pressures were significant and the effects of surface wave and convective acceleration on hydrodynamic pressure were about the same [7–9]. C.-H. Zee and R. Zee assumed that water was incompressible and modeled the water flow in reservoirs with a Laplace equation. They studied the hydrodynamic pressure on dams during a horizontal earthquake with a computer program [10]. Demirel and Aydin developed a wave-absorption filter for the far-end boundary of semi-infinite large reservoirs and computed the earthquake-excited surface waves and nonlinear hydrodynamic pressures in a dam-reservoir system based on solving the Navier-Stokes equations and a depth-integrated continuity equation. They used the numerical model to compute the maximum wave run-up height on the dam face during earthquakes. Based on the data obtained from numerical simulations for nine selected earthquake ground motions, they presented a relationship between dimensionless run-up height and the Froude number that was defined in terms of maximum ground velocity and water depth [11, 12].
Earthquake-excited reservoir systems have been modeled as a shaking-wave tank by many researchers. Faltinsen studied shaking-wave-tank problems with the boundary-element method, solving the Laplace equation based on potential-flow theory [13]. Cooker used linearized SWEs to model a horizontal rectangular wave tank. The free-surface boundary conditions were linearized by potential-flow theory [14]. They found that time-periodic solutions existed for the system. Liu and Lin developed a numerical model to study 3D nonlinear liquid sloshing with broken free surfaces under 6-degree-of-freedom excitations [15]. Godderidge et al. developed an inhomogeneous multiphase model for the computational-fluid-dynamics (CFD) analysis of violent sloshing flows [16]. Chen et al. solved the Reynolds-averaged Navier-Stokes (RANS) equations in both water and air regions to study the sloshing-tank problems [17]. Du et al. approached the characteristics of seismic surges using a shaking-water-tank experiment. They proposed two empirical relationships between the maximum water wave height and the peak ground acceleration for low- and high-frequency seismic waves, respectively [18].
Given the above, investigations of earthquake-induced water waves may apply three methods: analytical calculations, numerical simulations, and experimental studies. In the present study, for an inviscid and incompressible fluid, SWE is developed to simulate earthquake-induced water waves based on Chen’s hydrodynamic model [9]. A finite-difference method is used to discretize the SWE. The accuracy and applicability of our method were verified by the results of several classical models, which include analytical solutions of Sato [3] and numerical solutions of Demirel and Aydin [12]. Finally, we propose an empirical equation to predict the maximum water elevation of earthquake-induced water waves based on the calculated results of the proposed model.
2. Mathematical and Numerical Model
2.1. Basic Equations
Figure 1 is a sketch of the dam-reservoir system studied in this paper. As seen in the figure, , , and are space coordinates; , , and are the length, width, and height of the reservoir, respectively; is the water depth; is the still water depth; and is the free-surface displacement. The 3D equations of motion of fluid in this dam-reservoir system during an earthquake, proposed by Chen and Yuan [9], can be written as
The corresponding continuity equation of an incompressible fluid iswhere , , and are the fluid velocity components in the , , and directions, respectively; , , and are the instantaneous ground-velocity components the in the , , and directions, respectively; is the fluid density; is the gravity acceleration; and is the hydrodynamic pressure.
The 3D hydrodynamic model proposed by Chen focuses on the nonlinear hydrodynamic pressures on the banks during earthquakes and is calculated by a complete 3D finite-difference scheme. This work described in this paper focused on the maximum water elevation on the reservoir bank during earthquakes, which can be used for the critical hydrologic condition of geohazards such as a moraine-lake burst with overflow [19]. In order to simplify the calculation, we operated under the following hypotheses.
For the external seismic excitation, we assume that there is only the horizontal vibration without the vertical seismic motion, so
As there is only the horizontal vibration, the vertical acceleration of the fluid is far smaller than that of the gravity acceleration, so the assumption of nearly uniform vertical distribution of the hydrostatic pressure distribution is reasonable [20]; then,
The reservoir banks are vertical and rigid. There is no wave absorption and the water waves are totally reflected at the reservoir banks.
Based on the above assumptions, by depth integration of (1) and (2), the basic equations of the SWE model for simulation of earthquake water waves can be written as
2.2. Numerical Implementation
In this study, a finite-difference method was used to calculate (5). As shown in Figure 2, the parameters , , and are defined in the cell. The continuity equation and the momentum equations in the and directions, respectively, are discretized in the center of the cell at the th time step. The derivative terms of (5) can be discretized as shown in the following.
(1) Time-Derivative Terms
(2) Spatial-Derivative Terms
(3) The Middle Termswhere superscript denotes the time step, represents the time-step size, subscript denotes the spatial node, and and represent the spatial step sizes in the and directions, respectively. Considering (6)–(8), the continuity equation and the momentum equations can be discretized as follows:
2.3. Boundary Conditions
In this paper, we have assumed rigid and totally reflective boundary conditions. In order to treat these boundary conditions, one extra ghost grid is introduced outside the computational space domain [21]. The parameters of the boundary cell, namely, , , and , are defined at the center of the boundary cell. The boundary conditions in the direction are then updated usingand the boundary conditions in the y direction are updated by usingwhere superscript denotes the time step.
2.4. Stability Condition
The stability and accuracy of the numerical results significantly depend on the spatial step sizes and the selected time step. In this paper, the continuity equation and the momentum equations are discretized by an explicit difference scheme using a sequential iteration. Thus, the Courant-Friedrichs-Lewy (CFL) stability criterion is used [22], and the joint stability conditions for the time step are as follows:
2.5. Computational Cycle
With a chosen finite-difference mesh, and given initial conditions and ground acceleration for each time step, the complete cycle for the explicit iterative procedure can be summarized as follows.
Compute all the variables at the th time step, including variables at the center of the cell, such as , , and , and variables at the faces of the cell, such as , , and .
Compute the continuity equation (9) and the momentum equations (10) and (11) at the th time step to obtain the new water depths and the new water velocities and .
Update the boundary conditions with (12) and (13).
Now all of the new variables at the th time step have been obtained, so repeat Steps – and enter the next cycle.
3. Numerical Validation
To demonstrate the accuracy and applicability of the present model for earthquake-induced water waves, the calculated results of two classical models, which include analytical solutions of Sato and numerical solutions of Demirel and Aydin, will be used to facilitate a comparison to the results obtained using the present model. Three different seismic waves, namely, a Wolong wave, Mexicali wave, and Tongmai wave, having varied peak ground acceleration (PGA), are used in the comparison of calculated results.
3.1. Comparative Models
Two models of earthquake-induced water waves are used to facilitate the comparison to the present model.
(1) Sato’s Model (1967). Sato modeled the seismic wave as a simple sine wave and assumed that the earthquake-induced water waves were generated by horizontal movement of a vertical wall, and then proposed a formula for the earthquake-induced water wave using analytical calculations as follows:where is nondimensional, is the peak ground acceleration of the earthquake, and is the dominant frequency of the earthquake.
(2) Demirel-Aydin’s Model (2016). Demirel and Aydin used a numerical simulation method to compute breaking surface waves in the reservoir and to predict the maximum height of an earthquake-induced water wave. On the basis of these numerical simulation results, they presented the following equation for predicting the maximum wave height:where is the earthquake Froude number and is the peak ground velocity (PGV).
3.2. Earthquake Excitations
Three kinds of earthquake excitations are used in the comparison of calculated results. A Wolong seismic wave was recorded in the Wenchuan Ms 8.0 earthquake on 12 May 2008 at Wolong Station, as shown in Figure 3; a Mexicali seismic wave was recorded in the California Ms 7.2 earthquake on 4 April 2010 at Mexicali, as shown in Figure 4; and Tongmai seismic waves were the synthetic waves that were based on rock-response spectra (supplied by the Institute of Crustal Dynamic China and CCCC First Highway Consultants Co., Ltd., 2005), as shown in Figure 5. Details of the seismic waves and corresponding PGA and PGVs are shown in Table 1.
3.3. Calculated Results and Analysis
Before the numerical validation of the present model, the reservoir system is assumed to be a rectangular tank with a length of 300 m, a width of 150 m, and an initial water depth of 30 m. The simulations of earthquake-induced water waves based on the present model were then carried out with varied seismic waves at different output times, as seen in Figures 6(a)–6(f).
(a) Wolong wave, PGA = 0.5 g and s
(b) Wolong wave, PGA = 0.5 g and s
(c) Mexicali wave, PGA = 0.25 g and s
(d) Mexicali wave, PGA = 0.25 g and s
(e) Tongmai wave, PGA = 0.4 g and s
(f) Tongmai wave, PGA = 0.4 g and s
Figures 6(a) and 6(b) show the numerical water waves at different output times ( and 60 s) generated by a Wolong seismic wave with a 0.5 G horizontal acceleration. It can be seen in Figure 3 that the peak value of the acceleration of the Wolong seismic wave appeared at approximately s and the main quake ended at approximately s. Figures 6(c) and 6(d) show the numerical water waves at different output times ( and 30 s) generated by a Mexicali seismic wave with a 0.25 G horizontal acceleration. It can be seen in Figure 4 that the peak value of the acceleration of the Mexicali seismic wave appeared at approximately s and the main quake ended at approximately s. Figures 6(e) and 6(f) show the numerical water waves at different output times ( and 20 s) generated by a Tongmai seismic wave with a 0.25 G horizontal acceleration. It can be seen in Figure 5 that the peak value of the acceleration of the Tongmai seismic wave appeared at approximately s and the main quake ended at approximately s.
By comparing the water waves generated by the three kinds of seismic waves, it can be seen that, compared to the other two seismic waves, the wavelength of water waves generated by a Wolong seismic wave is shorter, the amplitude of water waves is smaller, and the water waves generated by a Wolong seismic wave calm down faster. These numerical results are reasonable. On one hand, the dominant frequency of a Wolong seismic wave is higher than the other two. During earthquakes, water waves are generated by the shaking bank, which is a kind of forced-vibration problem, and the frequency of water waves depends on the dominant frequency of the seismic wave. On the other hand, the dominant frequencies of all of the seismic waves are higher than the natural frequencies of reservoir water, which are very low. Then, for higher excitation frequency, the amplitude of water waves is smaller and decays faster [23].
3.4. Comparison to Former Models
The value of the present model, which will be used for comparison with former models, is the maximum water elevation at the right-hand side of the rectangular tank, as shown in Figure 6. In order to obtain the calculated results of the other models, the data of seismic waves in Table 1 and (15) and (16) are used. The initial water depth varies from 20 to 40 m at intervals of 5 m. Figure 7 depicts the comparison of the maximum water elevations under varied conditions obtained in the present model, in Sato’s model, and in the Demirel-Aydin model. As seen in Figure 7, the calculated values of the present model are all between the calculated results of the other two models [except in Figure 7(l)]. The results indicate that the calculated results of the present model are within the scope of the results of previous models, so the present model is reliable. The calculated results of the present model are closer to the results of the Demirel-Aydin model in most cases. This is because the present model and the Demirel-Aydin model both consider the time-history curve of an earthquake wave’s acceleration and use the seismic wave data more sufficiently. However, Sato’s model takes the seismic wave as a simple sine wave, may ignore more information, and could become more sensitive to the dominant frequency of an earthquake and to the PGA.
(a) Wolong, 0.1 g
(b) Wolong, 0.2 g
(c) Wolong, 0.3 g
(d) Wolong, 0.4 g
(e) Wolong, 0.5 g
(f) Mexicali, 0.05 g
(g) Mexicali, 0.1 g
(h) Mexicali, 0.15 g
(i) Mexicali, 0.2 g
(j) Mexicali, 0.25 g
(k) Tongmai, 0.05 g
(l) Tongmai, 0.25 g
(m) Tongmai, 0.4 g
3.5. Prediction of the Maximum Water Elevation
By comparing the calculated results of the three models, it can be seen that the calculated results of the present model are closer to those of the Demirel-Aydin model in most cases. However, while the dominant frequency of an earthquake is lower, the calculated results of the present model are also close to the results of Sato’s model. The Demirel-Aydin model considers that the PGV and the initial water depth are the main factors that influence the maximum water elevation of earthquake-induced water waves, while ignoring the influence of the PGA and the dominant frequency of the earthquake. Sato’s model simplified the seismic wave, which may increase the influence of the PGA and the dominant frequency of an earthquake.
In this paper, the PGA, PGV, dominant frequency of an earthquake, and initial water depth are considered to be the main factors that influence the maximum water elevation of earthquake-induced water waves. As there is a significant correlation among the PGA, PGV, and the dominant frequency of an earthquake [24], the PGA, PGV, and initial water depth are finally selected as the factors influencing the maximum water elevation. The three factors can be combined into two dimensionless parameters that can be calculated as follows:
The dimensionless wave height for the maximum water elevation of earthquake-induced water waves is defined by
In Section 3.4, we described the 65 group tests with varied conditions for earthquake-induced water waves that were carried out using the present model; the calculated results for the maximum water elevations are given in Table 2.
An empirical equation can then be developed for the dimensionless wave height by means of nondimensional multiple linear regression analysis. The relationship between the dimensionless wave height and the two dimensionless factors is expressed as follows:In Figure 8, we compare the results calculated using the present model and those using (19) for the 65 group tests. The correlation coefficient between the two sets of calculated results is 0.9045. Therefore, (19) can be used to estimate the maximum water elevation of earthquake-induced water waves if no better method is available.
4. Conclusions
In this study, a 2D SWE model was used to simulate earthquake-induced water waves. We used two former classical models to verify the accuracy and reasonability of our new model and drew the following conclusions from the results of our study.
It can be seen from the numerical results of the present model that, for a higher dominant frequency of an earthquake, the wavelength of water waves induced by an earthquake is shorter and their amplitude is smaller and decays faster.
The results calculated using the present model are within the scope of the results of previous models, which indicates that the present model is reliable. Furthermore, they are closer to results calculated using the Demirel-Aydin model in most cases. This is because the present model and the Demirel-Aydin model both consider the time-history curve of the earthquake wave’s acceleration and use the seismic wave data more sufficiently.
An empirical equation for the maximum water elevation of earthquake-induced water waves was developed based on the results calculated using the present model. The empirical equation selects the PGA, PGV, and initial water depth as factors influencing the maximum water elevation, which is an improvement to both Sato’s and the Demirel-Aydin model.
Ultimately, we have provided a new numerical model to calculate earthquake-induced water waves, and its calculation method is reliable. An empirical equation for the maximum water elevation of earthquake-induced water waves was developed that is an improvement to the former models. Moreover, our calculation method can greatly improve preventive engineering and provide references for disaster prevention and mitigation engineering.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (Grant no. 41571004), National Program on Key Research Project of China (Grant no. 2016YFC0802206), and Key Laboratory Foundation of Mountain Hazards and Earth Surface Processes, Chinese Academy of Sciences (Grant no. KLMHESP-17-03).