# Numerical and Statistical Analyses of Tsunami Heights with the L-Moments Method
## 2 Numerical Simulation
### Governing Equations  
We numerically modeled an offshore tsunami, where the linear shallow-water theory could be applied because the water surface displacement is relatively small compared with the water depth [5]. As a tsunami can propagate a long distance, the physical frequency dispersion, which significantly affects the propagation, must be considered. Therefore, we used linear Boussinesq equations that can consider physical frequency dispersion as the governing equations, following [PERSON] et al. [6]. The propagation model, except for the nonlinear terms and the Coriolis force, consisted of the continuity equation and momentum equations as follows:  
\[\frac{\partial\zeta}{\partial t}+\frac{\partial p}{\partial x}+\frac{\partial Q }{\partial y}=\ 0 \tag{1}\]  
\begin{table}
\begin{tabular}{c c c c c} \hline \hline
**Region** & **A** & **B** & **C** & **D** \\ \hline Boundary condition & Fully reflected & Fully reflected & Fully reflected & Fully reflected \\ Boundary condition & Free transmission & Dynamic linking & Dynamic linking & Dynamic linking \\ \hline \hline \end{tabular}
\end{table}
Table 2: Boundary conditions for each region.  
Figure 3: Virtual earthquake sources located near the western coast of Japan.  
\[\frac{\partial P}{\partial t}+gh\frac{\partial\zeta}{\partial x} = 0 \tag{2}\] \[\frac{\partial Q}{\partial t}+gh\frac{\partial\zeta}{\partial x} = 0 \tag{3}\]  
where \(\zeta\) is the free surface displacement, \(h\) is the still-water depth, \(P(=uh)\) and \(Q(=vh)\) are the depth-integrated volume fluxes in the \(x\) and \(y\) directions, respectively, and \(g\) is the gravitational acceleration. Equation (1) is the continuity equation with no truncation error. Equations (2) and (3) are the momentum equations that do not consider the hydrodynamic pressure.  
Using a modified leap-frog scheme technique, we employed the following propagation model to generate numerical dispersion consistent with the physical dispersion term of the linear Boussinesq equations, without changing the lattice spacing and the computation time interval, even when the water depth changes [6].  
\[\frac{\zeta_{i,j}^{n+1/2}-\zeta_{i,j}^{n-1/2}}{\Delta t}+\frac{P_{i+1/2,j}^{n}- P_{i-1/2,j}^{n}}{\Delta x}+\frac{Q_{i,j+1/2}^{n}-Q_{i,j-1/2}^{n}}{\Delta y}= 0 \tag{4}\]  
\[\begin{array}{c}\frac{P_{i+1/2,j}^{n+1}-P_{i+1/2,j}^{n}}{\Delta t}+gh_{i+1/2,j}\frac{\zeta_{i+1,j}^{n+1/2}-\zeta_{i,j}^{n+1/2}}{\Delta x}+\frac{\alpha}{12 \Delta x}\ gh_{i+1/2,j}\left[\zeta_{i+2,j}^{n+1/2}-3\zeta_{i+1,j}^{n+1/2}+3 \zeta_{i,j}^{n+1/2}-\zeta_{i-1,j}^{n+1/2}\right]\\ +\frac{\gamma}{12\Delta x}gh_{i+1/2,j}\left[\left(\zeta_{i+1,j+1}^{n+1/2}-2 \zeta_{i+1,j}^{n+1/2}+\zeta_{i+1,j-1}^{n+1/2}\right)-\left(\zeta_{i,j+1}^{n+1/ 2}-2\zeta_{i,j}^{n+1/2}+\zeta_{i,j-1}^{n+1/2}\right)\right]= 0\\ \frac{Q_{i,j+2}^{n+1/2}-Q_{i,j+1/2}^{n}}{\Delta t}+gh_{i,j+1/2}^{n}\frac{ \zeta_{i+1/2}^{n+1/2}-\zeta_{i,j}^{n+1/2}}{\Delta y}\\ +\frac{\alpha}{12\Delta y}gh_{i,j+1/2}\left[\zeta_{i,j+2}^{n+1/2}-3\zeta_{i,j+ 1}^{n+1/2}+3\zeta_{i,j}^{n+1/2}-\zeta_{i,j-1}^{n+1/2}\right]\\ +\frac{\gamma}{12\Delta y}gh_{i,j+1/2}\left[\left(\zeta_{i+1,j+1}^{n+1/2}-2 \zeta_{i,j+1}^{n+1/2}+\zeta_{i-1,j+1}^{n+1/2}\right)-\left(\zeta_{i+1,j}^{n+1 /2}-2\zeta_{i,j}^{n+1/2}+\zeta_{i-1,j}^{n+1/2}\right)\right]= 0\end{array} \tag{6}\]  
In [PERSON] et al. [6], the variance correction factors \(\alpha\) and \(\gamma\) are determined as follows, and the numerical results are the same as those of the linear Boussinesq equation.  
\[\alpha = \frac{4h^{2}+gh(\Delta t)^{2}-(\Delta x)^{2}}{(\Delta x)^{2}}, \ \gamma\ =\ \alpha\ +\ 1 \tag{7}\]  
Detailed information on the numerical scheme, including the leap-frog scheme, can be found in the literature [1; 4].