Rapid forecasting of tsunami runup heights from 2-D numerical simulations

We propose a method to compute tsunami runup heights that is based on an integration of numerical, 2-D shallow-water modelling and an analytical, 1-D long-wave runup theory. This approach provides a faster forecast of tsunami runup heights than a complicated coastal inundation model. Through simulations of potential tsunami scenarios, this approach can also be applied to long-term tsunami prediction. We tested the model by simulating the historical event in the East (Japan) Sea and found that the estimates of runup heights agreed well with the available observations.


Introduction
Numerical modeling of the propagation of tsunami waves is an important tool for forecasting tsunami heights and the risks posed to coastal populations.Previous analyses of tsunami characteristics have used the shallow-water model with a no-flux boundary condition at a depth of 5-10 m on the shore (Titov and Synolakis, 1998;Sato et al., 2003;Choi et al., 2003;Ioualalen et al., 2007;Schuiling et al., 2007;Zaitsev et al., 2009;Beisel et al., 2009).In this model, runup heights are determined using simplified formulae of the 1-D, analytic theory of long-wave runup for a fixed shape of an incident wave, whether a sine wave or a solitary wave (Choi et al., 2002a;Ward and Asphaug, 2003).The runup stage is included in the numerical model by taking into account the various assumptions concerning the hydraulic properties (roughness) of the dry land, an important calculation for tsunami risk assessment (Glimsdal et al., 2006;Lovholt et Correspondence to: E. Pelinovsky (pelinovsky@hydro.appl.sci-nnov.ru) al., 2006;Dominey-Howes and Papathoma, 2007;Gonzalez et al., 2009;Dall'Osso et al., 2010;Gayer et al., 2010).The direct calculation of the propagation of the tsunami wave from its source to the coastal zone using a single numerical model results in low accuracy.As a partial solution, various nested methods with different mesh resolutions in the open sea and the coastal zone have been developed (Titov andSynolakis, 1995, 1998;Choi et al., 2003;Wei et al., 2008;Roger and Hebert, 2008;Roger et al., 2010;Titov et al., 2011).Near the coast, accurate computing of tsunami wave dynamics requires small grid steps of 10-100 m.As a result, such numerical models are difficult to use in operational practice.For instance, in the East (Japan) Sea, tsunamis usually arrive at the eastern Korean coasts within 100 min.Thus, the runup height in coastal zones must be predicted well in advance, and it is for this reason that the Korea Meteorological Administration (KMA) uses the vertical wall approximation in its model of wave propagation.The model permits rapid warning, but it also yields low accuracy in its estimates of tsunami runup height.
In practice, the formula describing the runup of sine and solitary waves on a planar beach is used to approximate the tsunami runup heights; see, for instance, Choi et al. (2002a) and Ward and Asphaug (2003).It has been shown recently that the formulae for wave runup with various shapes can be parameterised and unified for incident waves of various shapes if the characteristic wavelength of the incident wave is determined on the level 2/3 from maximal wave height (Didenkulova and Pelinovsky, 2008;Didenkulova, 2009).Defining the incident wave near the coast, where the bottom profile can be approximated by a planar beach, is not straightforward if the computation is performed in a domain with a non-reflected boundary condition.
The primary purpose of this study is to develop a model that combines numerical simulations for wave dynamics far from the coast with analytical solutions for the wave runup.The novelty of our approach, as contrasted to previous studies (Choi et al., 2002a;Ward and Asphaug, 2003), lies in its original expression for the runup ratio.This new approach facilitates the rapid prediction of a tsunami's characteristics upon its arrival at the coast to help prevent tsunami-related disasters.Section 2 describes an analytical approach to compute the runup height using the wave height at the wall; this approach is based on a rigorous solution of a set of 1-D shallow-water equations for waves above a planar beach.Section 3 presents a shallow-water numerical model for simulating tsunami propagation in the East (Japan) Sea.The combined model is tested in Sect.4, which describes a numerical simulation of the 1993 event and calculations of the expected runup heights.Estimations of the runup heights along the Eastern Korean Coast are in close agreement with observations.
2 Analytical approach to compute the runup height through the height of a tsunami wave at the wall In the analytical theory of long-wave runup, rigorous solutions of 1-D, nonlinear, shallow-water equations can be obtained for a planar beach using the Carrier -Greenspan transformation (Carrier and Greenspan, 1958) for the various shapes of the incident wave (Pedersen and Gjevik, 1983;Synolakis, 1987;Kaistrenko et al., 1991;Synolakis, 1991;Pelinovsky and Mazova, 1992;Tadepalli and Synolakis, 1994;Carrier et al., 2003;Kânoglu, 2004;Tinti and Tonini, 2005;Kânoglu and Synolakis, 2006;Didenkulova et al., 2006Didenkulova et al., , 2008;;Antuono and Brocchini, 2007, 2008, 2010;Didenkulova and Pelinovsky, 2008;Madsen andFuhrman, 2008, Didenkulova, 2009;Briganti and Dodd, 2009;Dobrokhotov and Tirozzi, 2010).The important result here is that the linear and nonlinear theories predict the same maximal values for the runup height if the incident wave is determined far from the shore, and for this reason, linear theory can be applied to analyse the runup process.Within the framework of linear theory, rigorous solutions can be found for various bottom profiles, not only the planar beach.A popular approximation of the bottom profile is a planar beach combined with a flat bottom (a configuration usually applied in laboratory modeling).In this case, the incident and reflected waves are easily separated on the flat bottom, and the runup height of the monochromatic wave is (Shuto, 1972) where R is the maximal runup height, A is the amplitude of the incident wave, L is the shelf width, k is the wave number of the incident wave, and J 0 and J 1 are Bessel's functions.Using Fourier transformations, similar formulae (at least in integral form) can be derived for an incident wave of arbitrary shape (Synolakis, 1987;Pelinovsky and Mazova, 1992).The amplification factor (1) has been used to estimate the runup heights observed on the Korean Coast during the 1993 tsunami (Choi et al., 2002a).When applied to tsunamis, the bottom profile varies in the coastal zone and cannot be represented as flat.A classic numerical model includes the wall at depth h and distance L from the shore (Fig. 1).Numerical simulation can then capture the resulting oscillation of water level, η(L,t).Selecting the incident and reflected waves from the wave records near the wall is not simple; here, we apply another approach based on rigorous solutions of the wave equation.
The shallow-water equations in the linear approximation can be reduced to the wave equation for the water displacement, where g is the acceleration due to gravity, and h(x) describes the bottom profile.The general solution of Eq. ( 2) can be found using the Fourier transformation The Fourier transform η(x,ω) is the solution of the ordinary differential equation following from Eq. ( 2), We assume that h(x) =αx everywhere and α is constant.For this case, an elementary solution of (4) is presented through the Hankel functions Using the asymptote of the zero-order terms of the Hankel function, it can be shown easily that Eq. ( 5) corresponds to the incident wave with amplitude A (propagated onshore), and the first corresponds to the reflected wave with amplitude B (propagated offshore).We assume that characteristics of the incident wave are known, and therefore, A(ω) can be found from the Fourier presentation of the incident wave.
We first analyse the wave runup on the vertical wall (Fig. 1).The boundary condition on the wall is which leads to where T is the travel time of the wave from the wall on depth h to the shore, Equation ( 5) for x = L and ( 7) can serve as a system of equations to find the amplitudes of the incident and reflected waves, It is easy to simplify these equations, taking into account that the denominator is the Wronskian of the Hankel functions [see Eq. (9.1.17)in Abramowitz and Stegun (1964)]: iπ ωT H (1) Thus, the equations (10) can be used to compute the spectral amplitudes of the incident and reflected waves through the Fourier spectrum of the water oscillations on the wall.We now consider the wave runup on the same planar beach with no vertical wall, assuming that the incident monochromatic wave has the same amplitude A(ω) as before.In this case, the bounded (on the shore) solution to (5) is Using ( 11) and (9) to eliminate A(ω), we can calculate the wave field on a planar beach compared to the wave oscillation on the wall.In particular, the spectral amplitude of the water oscillations on the shore is η(0,ω) = iπ ωT H (1) The expression ( 12) can be considered the runup ratio for a sine wave, which is discussed in Harbitz and Pedersen (1992) and Kaistrenko et al. (1999) on a planar beach.Using the Fourier transformation of (12), we can compute η(0,t) through η(L,t).The inverse Fourier transformation in this asymptotic case is expressed by the simple integral Equation ( 13) can be integrated by parts and written as In Eqs. ( 13) and ( 14), t =0 corresponds to the wave approaching the vertical wall.In the initial moment, it is assumed that η(x = L,t =0) = dη(x = L,t =0)/dt = 0.The input function η(x = L,t) gives the water wall oscillations computed from the 2-D model.The maximum of η(x =0,t) in linear theory gives the maximal runup height of tsunami waves on the coast in nonlinear theory.Thus, considering the two geometries of a nearshore planar beach (in the runup study) and a planar beach with a vertical wall at a fixed depth (the equivalent boundary condition), we have shown that the runup height can be expressed through the characteristics of the water oscillations on the vertical wall.This procedure is rapid and easily realised on computers.As a result, the numerical simulation of tsunami waves can be used for estimations of the runup heights of tsunami waves.
Several limitations of the proposed approach should be mentioned.First, the 1-D analytical runup theory is applied, whereas the numerical simulation of tsunami waves in real basins assumes 2-D geometry.This distinction means that 2-D effects (refraction, focusing, etc.) nearshore are not included in the model.The wave field on the wall generally contains the approaching waves propagated onshore and the trapped waves propagated alongshore.It is evident that the first waves in the time series near the wall are "onshore" waves (trapped waves approach later), and thus the interval of integration in Eqs. ( 13) or ( 14) should be bounded by the few first waves.The second limitation concerns the use of the approximation of the nearshore bottom profile and the vertical wall as a planar beach with the same slope.Real bathymetry near the wall can be approximated by the piecewise linear function (Kânoglu and Synolakis, 1998) or by a power function (Didenkulova et al, 2009;Didenkulova and Pelinovsky, 2010).It is evident that near-shore bathymetric features influence the runup process and may affect the results described above.Another limitation is the ignoring of breaking effects in tsunami waves.Finally, the Eqs.( 13) or ( 14) contain linear approximation.However, as noted earlier, the maximal runup height computed by linear and nonlinear theories is the same, suggesting that this last limitation is less important.The appropriateness of these assumptions can be tested through modeling historic events.In the next section the developed analytical approach that will be used in estimations of the 1993 tsunami runup heights in the East (Japan) Sea.A general description of this event is given in (Choi et al., 1994(Choi et al., , 2003)).

Numerical model
A finite-difference model (Choi et al., 2003) was developed to simulate the generation and propagation of the tsunami.The numerical model uses the linear shallow-water equation with a spherical coordinate system covering the East Sea area with mesh dimensions 959 × 1118, a mesh size of 1 angular min, and a 2-s time step (Fig. 2).The basic equations are In the equations above, φ and χ are latitude and longitude, P and Q are discharge per unit width in the direction of φ and χ , respectively, R is the radius of the earth, and f is the Coriolis parameter.We compiled all of the available topographic data with a digitisation of hydrographic charts, including Chinese and  Russian navigation charts and Korean hydrographic charts.These data were used to produce a 1-min gridded topography and bathymetry dataset for the region (117 • E to 143 • E longitude and 24 • N to 52 • N latitudes) (Fig. 2) with a refined 1 arc-second topography and bathymetry west of 135 • East longitude (Choi et al., 2002b).The Hokkaido-Nansei-Oki earthquake occurred at 13:17 UTC on 12 July 1993.The epicentre of this earthquake was in the East (Japan) Sea at 42.8 • N, 139.2 • E, about 15-30 km off the small offshore island of Okushiri along the west coast of Hokkaido (Fig. 3).The magnitude of the earthquake was 7.8.Over 230 people were killed, and a majority of deaths were caused by the tsunami.Property losses were approximately $600 million.For the 1993 Hokkaido tsunami, the fault parameters suggested by Takahashi et al. (1994) are in Table 1.Here L is fault length, W is fault width, U is dislocation, h f is focal depth, θ is strike angle, δ is dip angle and λ is slip angle.The initial surface profile at the source was determined by the method of Manshinha and Smylie (1971).
The boundary conditions near the coast are "no-flux": the normal component of the particle velocity or flow discharge to the boundary is zero, corresponding to a vertical wall at the last sea points.We previously used this numerical shallow-water model in our earlier studies of tsunami waves in the East (Japan) Sea (Choi et al., 2002a(Choi et al., , 2003)).
Figure 4 shows the distribution of the maximum water elevations in the East Sea for the 1993 tsunami and the general direction of the tsunami's propagation.We can infer that the direction of the tsunami wave's propagation was slightly affected by bathymetry, because the general direction of propagation was to the north and to the Russian coasts.This movement focused the tsunami's energy away from the eastern edge of the East (Japan) Sea and the mid-eastern coast of Korea.

The distribution of tsunami runup heights along the eastern Korean coast
The tsunami that occurred on 12 July 1993 was well documented along the eastern Korean coast (Choi et al., 1994;Choi et al., 2002a).These observations were used to check the accuracy of the combined model.Figure 5 displays the time series of the water displacement along the vertical wall as computed by the finite-difference model (dashed line).The displacement of the water on the shore (solid line) was calculated using the integral in Eq. ( 14) for the coastal locations shown in Fig. 6.The amplification factor of the tsunami waves nearshore is approximately 2 to 2.2.The observed maximal runup heights in these locations are given by numbers and by dash-point lines.The computations and observations agree well at all except two (the Hosan and Bugu Beaches) of the eight coastal locations.

Conclusions
Combining 2-D numerical models and 1-D analytical runup theory, we proposed a novel and rapid method for forecasting tsunami runups on the coasts.In the first step, the 2-D numerical simulations of tsunami generation and propagation are performed assuming impermeable boundary conditions of a 5-10 m depth at the last sea points (equivalent to a wall).
Then the time series of the water oscillations on the wall are used to calculate the runup heights using the analytical integral expression following from 1-D theory.The accuracy of this approach was validated against the historic tsunami event in the East (Japan) Sea; the observed runup heights are in close agreement with the predicted values.We have demonstrated that the proposed approach is a more reasonable and rapid method of forecasting than complicated coastal inundation models.This new method can also be applied to other regions.

Fig. 2 .
Fig. 2. Topography and bathymetry (contour unit: metre) of the East (Japan) Sea.The 8 dots along the Korean coast indicate beaches where comparisons between the observed and calculated runup heights were made (see Fig. 6).

Fig. 3 .
Fig. 3.The location and initial water elevations of the 1993 event.Solid and dashed lines indicate the positive and negative values, respectively.The contour interval is 0.5 m.

Fig. 4 .
Fig. 4. The distribution of maximum wave heights (cm) in the East (Japan) Sea for the 1993 tsunami.

Fig. 6 .
Fig. 6.Comparison of the computed and observed runup heights on the eastern Korean coast.

Table 1 .
Fault parameters of the 1993 Hokkaido earthquake.