Effective coastal boundary conditions for tsunami wave run-up over ...

Nonlin. Processes Geophys., 21, 987–1005, 2014 www.nonlin-processes-geophys.net/21/987/2014/ doi:10.5194/npg-21-987-2014 © Author(s) 2014. CC Attribution 3.0 License.

Effective coastal boundary conditions for tsunami wave run-up over sloping bathymetry W. Kristina1 , O. Bokhove1,2 , and E. van Groesen1,3 1 Department

of Applied Mathematics, University of Twente, Enschede, the Netherlands of Mathematics, University of Leeds, Leeds, UK 3 LabMath-Indonesia, Bandung, Indonesia 2 School

Correspondence to: W. Kristina ([email protected]), O. Bokhove ([email protected]), and E. van Groesen ([email protected]) Received: 10 February 2014 – Published in Nonlin. Processes Geophys. Discuss.: 21 March 2014 Revised: 19 July 2014 – Accepted: 28 July 2014 – Published: 23 September 2014

Abstract. An effective boundary condition (EBC) is introduced as a novel technique for predicting tsunami wave runup along the coast, and offshore wave reflections. Numerical modeling of tsunami propagation in the coastal zone has been a daunting task, since high accuracy is needed to capture aspects of wave propagation in the shallower areas. For example, there are complicated interactions between incoming and reflected waves due to the bathymetry and intrinsically nonlinear phenomena of wave propagation. If a fixed wall boundary condition is used at a certain shallow depth contour, the reflection properties can be unrealistic. To alleviate this, we explore a so-called effective boundary condition, developed here in one spatial dimension. From the deep ocean to a seaward boundary, i.e., in the simulation area, we model wave propagation numerically over real bathymetry using either the linear dispersive variational Boussinesq or the shallow water equations. We measure the incoming wave at this seaward boundary, and model the wave dynamics towards the shoreline analytically, based on nonlinear shallow water theory over bathymetry with a constant slope. We calculate the run-up heights at the shore and the reflection caused by the slope. The reflected wave is then influxed back into the simulation area using the EBC. The coupling between the numerical and analytic dynamics in the two areas is handled using variational principles, which leads to (approximate) conservation of the overall energy in both areas. We verify our approach in a series of numerical test cases of increasing complexity, including a case akin to tsunami propagation to the coastline at Aceh, Sumatra, Indonesia.

1

Introduction

Shallow water equations are widely used in the modeling of tsunamis, since their wavelengths (typically 200 km) are far greater than the depth of the ocean (typically 2–3 km). Tsunamis also tend to have a small amplitude offshore, which is why they generally are less noticeable at sea. Linear shallow water equations (LSWE) therefore often suffice as a simple model of tsunami propagation (Choi et al., 2011; Liu et al., 2009; Kâno˘glu and Synolakis, 1998). On the contrary, it turns out that the lack of dispersion is a shortcoming of shallow water modeling when the tsunami reaches the shallower coastal waters on the continental shelf, and thus dispersive models are often required (Madsen et al., 1991; Horrillo et al., 2006). Numerical simulations based on these linear models are desirable because they involve a small amount of computation. However, as the tsunami approaches the shore, shoaling effects cause a decrease in the wavelength and an increase in the amplitude. Here, the nonlinearity starts to play a more important role, and thus the nonlinear terms must be included in the model. To capture these shoaling effects in more detail, a smaller grid size will be needed. Consequently, longer computational times are required. Some numerical models of tsunamis use nested methods with different mesh resolutions to preserve the accuracy of the solution near the coastal area (Titov et al., 2011; Wei et al., 2008), while other models employ an impenetrable vertical wall at a certain depth contour as the boundary condition. Obviously, the reflection properties of such a boundary condition can be unrealistic. We therefore wish to alleviate

Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.

988

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry

this shortcoming by an investigation of a so-called effective boundary condition (EBC) (Kristina et al., 2012), including run-up. In one horizontal spatial dimension, an outline of the desired mathematical modeling is sketched in Fig. 1. In the deep ocean for x ∈ [B, L] with horizontal coordinate x and seaward boundary point x = B, denoted as the simulation area, we model the wave propagation numerically using a linear model. In the coastal zone for x ∈ [xs (t), B] with shoreline position xs (t) < B, denoted as the model area, we model the wave propagation analytically using a nonlinear model by approximating the bathymetry as a planar beach. We calculate the run-up heights at the shore and the reflection caused by the slope. The reflected wave is then influxed back into the simulation area using the EBC. The coupling between the numerical and analytic dynamics in the two areas is handled using variational principles, which leads to (approximate) conservation of the overall energy in both areas. Following Kristina et al. (2012), an observation and influx operator are defined at x = B to measure the incoming wave signal and to influx the reflected wave, respectively. The shoreline position and wave reflection in the model area (sloping region) are determined using an analytical solution of the nonlinear shallow water equations (NSWE) following the approach of Antuono and Brocchini (2010) for unbroken waves. The decomposition of the incoming wave signal and the reflected one is also described in Antuono and Brocchini (2007, 2010) for the calculation of the shoreline and the wave reflection. Nevertheless, the method in their paper is applied by determining the incoming wave signal with the solution of the Korteweg–de Vries (KdV) equation. The novelty of our approach is the utilization of an observation operator at the boundary x = B to calculate the incoming wave elevation towards the shore from the numerical solution of the LSWE in the simulation area. For any given wave profile and bathymetry in the simulation area, the numerical solution can be calculated, and the signal arriving at x = B can be observed. Afterwards, the data are used to calculate the analytical solution of the NSWE in the onshore region and the reflected waves. A rapid method for estimating tsunami run-up heights is also suggested by Choi et al. (2011, 2012) by imposing a hard-wall boundary condition at x = B. Giving the water wave oscillations at this hard wall at x = B, the maximum run-up height of tsunami waves at the coast is subsequently calculated separately by employing a linear approach. It is claimed that the linear and nonlinear theories predict the same maximal values for the run-up height if the incident wave is determined far from the shore (Synolakis, 1987). In contrast, Li and Raichlen (2001) show that there is a difference in the maximum run-up prediction between linear and nonlinear theory. In addition to calculating only the maximum run-up height as in Choi’s method, our EBC also includes the calculation of reflected waves. The pointwise wave height in the whole domain (offshore and onshore area) is thus predicted accurately. For the inundation predicNonlin. Processes Geophys., 21, 987–1005, 2014

Figure 1. At the seaward boundary x = B, we assign (η, u) data, and we want to find a solution of the NSWE on the sloping region near the shoreline.

tion, we have verified that the method introduced by Choi et al. (2011, 2012) performs as well as our EBC method, while the reflection wave comparisons show larger discrepancies due to the usage of a hard-wall boundary condition. The interaction between incoming and reflected waves needs to be predicted accurately, since subsequent waves may cause danger at later times. Stefanakis et al. (2011) investigate the fact that resonant phenomena between the incident wavelength and the beach slope are found to occur. The resonance happens due to incoming and reflected wave interactions, and the actual amplification ratio depends on the beach slope. It explains why in some cases it is not the first wave that results in the highest run-up. Determination of the location of the seaward boundary point x = B is another issue that must be addressed. Choi et al. (2011) put the impermeable boundary conditions at a 5–10 m depth contour. In comparison, Didenkulova and Pelinovsky (2008) show that their run-up formula for symmetric waves gives optimal results when the incoming wave signal is measured at a depth that is two-thirds of the maximum wave height. We determine the location of this seaward boundary as the point before the nonlinearity effect arises, and examine the dispersion effect at that point as well. Considering the simple KdV equation (Mei, 1989), the measures of nonlinearity and dispersion are given by the ratios  = A/ h and µ2 = (kh)2 for the wave amplitude A, water depth h, and wavenumber k. Provided with the information of the initial wave profile, we can calculate the amplification of the amplitude and the decrease in the wavelength in a linear approach, and thereafter estimate the location of the EBC point. The EBC in this article will be derived in one spatial dimension for reasons of simplicity and clarity of exposure. The numerical solution in the simulation area is based on a variational finite element method (FEM). In order to verify the EBC implementation that employs this asymptotic closed-form solution, we also numerically simulate the NSWE in the model area using a finite volume method www.nonlin-processes-geophys.net/21/987/2014/

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry (FVM). Both cases are coupled to the simulation area to compare the results. We also validate our approach against the laboratory experiment of Synolakis (1987). In Sect. 2, we introduce the linear variational Boussinesq model (LVBM) and shallow water equations (SWE), both linear and nonlinear, from their variational principles. The coupling conditions required at the seaward boundary point are also derived here. The solution of the NSWE using a method of characteristics is shown in Sect. 3, which includes the solution of the shoreline position. In Sect. 4, the effective boundary condition is derived. It pinpoints the newly derived coupling conditions between the finite element simulation area and the model area. Numerical validation and verification are shown in Sect. 5, and we conclude in Sect. 6. 2

Water wave models

Our primary goal is to model the water wave motion to the shore analytically, instead of resolving the motion in these shallow regions numerically. We therefore introduce an artificial open boundary at some depth, and wish to determine an effective boundary condition at this internal boundary. Consider motion in a vertical plane normal to the shore, with an offshore coordinate x. The artificial boundary is then placed at x = B, while the real (time-dependent) boundary lies at x = xs (t) with xs (t) < B. For example, at rest, land starts at x = 0, where the total water depth h(0, t) is zero. This water line is time dependent, as the wave can move up and down the beach. We will restrict our attention to the dynamics in a vertical plane with horizontal and vertical coordinates x and z, respectively. Nonlinear potential flow water waves are succinctly described by variational principles as follows (Luke, 1967; Zakharov, 1968; Miles, 1977):

0=δ

The approximation for the velocity potential 8 in Eq. (1) can be of various kinds, but all are based on the idea of restricting the class of wave motions to a class that contains the wave motions one is interested in van Groesen (2006),Cotter and Bokhove (2010), and Gagarina et al. (2013). Following Klopman et al. (2010), we approximate the velocity potential as follows:

L [φ, 8, η, xs ] dt 0 ZT ZL

for a function F = F (z). Its suitability is determined by insisting that F (η) = 0, such that φ is the potential at the location z = η of the free surface and satisfies the slip flow condition at the bottom boundary z + hb (x) = 0. The latter kinematic condition yields ∂z 8 + ∂x 8∂x hb = 0 at z = −hb (x). For a slowly varying bottom topography, this condition is approximated by (∂z 8)z=−hb (x) = F 0 (−hb ) ψ = 0 .

=δ 0 xs

Zη −

 1 |∇8|2 dz dxdt 2

(1)

−hb

with velocity potential 8 = 8(x, z, t) and surface potential φ(x, t) = 8(x, z = η, t), where η = h − hb is the wave elevation and h = h(x, t) the total water depth above the bathymetry b = −hb (x), with hb (x) the rest depth. Time runs from t ∈ [0, T ]; partial derivatives are denoted by ∂t , etc., the gradient in the vertical plane by ∇ = (∂x , ∂z )T , and the acceleration of gravity by g.

www.nonlin-processes-geophys.net/21/987/2014/

(3)

Substitution of Eq. (2) into Eq. (1) yields the variational principle for Boussinesq equations, as follows (Klopman et al., 2010): ZT L [φ, ψ, η, xs ] dt

0=δ

0 ZT ZL 

=δ 0

 1 1  φ∂t η − g (h + b)2 − b2 − (η + hb ) |∂x φ|2 2 2 xs  1 1 2 2 ˘ ˘ x ψ| − γ˘ ψ dxdt , − β∂x ψ∂x φ − α|∂ (4) 2 2

˘ where functions β(x), α(x), ˘ and γ˘ (x) are given by Zη

−hb

  1  φ∂t η− g (h + b)2 − b2 2

(2)

8(x, z, t) = φ(x, t) + F (z)ψ(x, t)

˘ β(x)=

ZT

989

Zη F dz, α(x)= ˘ −hb

2

Zη

F dz, γ˘ (x)=

(F 0 )2 dz .

(5)

−hb

The shallow water equations (SWE) are derived with the assumption that the wavelengths of the waves are much larger than the depth of the fluid layer, so that the vertical variations are small and will be ignored. In this case, there is no dispersive effect. The velocity potential is approximated over depth by its value at the surface, such that F (z) = 0. Hence, when β˘ = α˘ = γ˘ = 0 in Eq. (4), the nonlinear shallow water equations are obtained as a limiting system. We a priori divide the domain into two intervals: x ∈ [B, L], where we model the wave propagation linearly, and x ∈ [xs (t), B], where we keep the nonlinearity. To be precise, in the simulation area x ∈ [B, L], we linearize the equations, and thus the wave propagation in this domain is modeled by the linear shallow water shallow water equations or the linear yet dispersive Boussinesq model. In the model area x ∈ [xs (t), B], we only consider depth-averaged shallow water flow. The non-dispersive and nonlinear shallow Nonlin. Processes Geophys., 21, 987–1005, 2014

990

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry

water equations are thus used to model the wave propagation in this region. Hereafter, we write φ˘ and η˘ for the linear variables. Consequently, by applying the corresponding approximations to variational principle (4), the (approximated) variational principle becomes

=

ZT hZL 0


˘ + ∂x (β∂ ˘ x ψ) ˘ δ φ˘ ∂t η˘ + ∂x (hb ∂x φ)

B

 ˘ + ∂x (β∂ ˘ x φ) ˘ − γ˘ ψ)δ ˘ ψ˘ dx − (∂t φ˘ + g η)δ ˘ η˘ + (∂x (α∂ ˘ x ψ) ˘ x ψ)δ ˘ φ| ˘ x=B + (α∂ ˘ x φ)δ ˘ ψ| ˘ x=B + (hb ∂x φ˘ + β∂ ˘ x ψ˘ + β∂

ZT h i ˘ ψ, ˘ η, ˘ φ, η, xs dt 0 = δ L φ,

+

(6a)


∂t η + ∂x ((η + hb ) ∂x φ) δφ − (∂t φ + gη

xs

0 ZT

hZL  1 1 ˘ t η˘ − g η˘ 2 − hb |∂x φ| ˘ 2 − β∂ ˘ x ψ∂ ˘ x φ˘ =δ φ∂ 2 2 0 B  1 1 ˘ 2 − γ˘ ψ˘ 2 dx ˘ x ψ| − α|∂ 2 2 B Z   1  + φ∂t η − g (h + b)2 − b2 2 xs  i 1 − (η + hb ) |∂x φ|2 dx dt . 2

 1 + ∂x2 φ)δη dx − (η + hb )∂x φδφ|x=B 2 i dxs + (φδη)|x=xs − (φ∂t η)|x=xs δxs dt, dt

Z0 α˘ = α(x) ˘ =

F 2 dz =

(6b)

˘ β˘ = β(x) =

γ˘ = γ˘ (x) =

(9a)

˘ + ∂x (β∂ ˘ x ψ) ˘ = 0, ∂t η˘ + ∂x (hb ∂x φ) ˘ + ∂x (β∂ ˘ x φ) ˘ − γ˘ ψ˘ = 0, ∂x (α∂ ˘ x ψ)

(9b) (9c)

and for x ∈ [xs (t), B], we get the nonlinear equations of motion 1 ∂t φ + gη + ∂x2 φ = 0, 2 ∂t η + ∂x ((η + hb ) ∂x φ) = 0.

(10a) (10b)

The last two terms in Eq. (8b) are the boundary terms at x = xs . They can be rewritten as follows:

0

ZT 

2 F dz = − hb , 3

=

−hb

Z0

∂t φ˘ + g η˘ = 0,

 ZT  dxs − (φ∂t η) |x=xs δxs dt (φδη) |x=xs dt

8 hb , 15

−hb

Z0

(8b)

where we used the endpoint conditions δη(0) = δη(T ) = 0, and no-normal through-flow conditions at x = L and h (xs (t), t) = 0. Since the variations are arbitrary, the linear equations emerging from Eq. (8b) for x ∈ [B, L] are as follows:

We choose a parabolic profile function F (z; hb ) = 2z/ hb + z2 / h2b , in which the x dependence is considered to be parametric when the total water depth h is sufficiently slowly varying. The coefficients in Eq. (5) simplify to their linearized counterparts in the simulation area where the linear Boussinesq equations hold (while these coefficients disappear in the model area where the nonlinear depth-averaged shallow water equations hold)

−φ∂x (η+hb ) 0

(F 0 )2 dz =

4 . 3hb

(7)

−hb

The variations in Eqs. (6) yield

1 0 = lim →0 

ZB 

0

Nonlin. Processes Geophys., 21, 987–1005, 2014

(8a)

(11)

since the total depth is h(xs , t) = η(xs , t) + hb (xs ) = 0 at the shoreline boundary. We therefore have the relation 0 = δh (xs , t) = δh + ∂x hδxs = δη + ∂x (η + hb ) δxs . Substituting Eq. (10b) into Eq. (11), the boundary condition at the shoreline is dxs = ∂x φ at x = xs (t) ; dt

ZT h ˘ ψ˘ + δ ψ, ˘ η˘ + δ η, L φ˘ + δ φ, ˘ φ + δφ ,

i h i ˘ ψ, ˘ η, η + δη, xs + δxs −L φ, ˘ φ, η, xs dt

  dxs − φ∂t η δxs dt, dt x=xs

(12)

i.e., the velocity of the shoreline equals the horizontal velocity of the fluid particle. The underlined terms in Eq. (8b) apply at the seaward point, where we want to derive the coupling of the effective boundary conditions. To derive the condition for the linear model, the goal is to write these terms www.nonlin-processes-geophys.net/21/987/2014/

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry ˘ Because the depth-averaged using the variations δ φ˘ and δ ψ. shallow water equations are considered, we have

991

(13)

in which h0 is the still water depth at the seaward boundary, and u0 , l0 , and t0 are defined below as s s γ gh0 h0 γ g ? h0 , (19) , l0 = ? , t0 = ? u0 = g γ γ g?

where the last equality arises from approximation (2) for the velocity potential. The variation of δφ thus becomes

where g = 1 and γ = 1 are dimensionless gravity acceleration and beach slope, respectively. The NSWE in dimensionless form are then given by

¯ φ(x, t) = 8(x, t) =

Z0

1 hb

8(x, z, t)dz = φ˘ +

β˘ ψ˘ , hb

−hb

δφ = δ φ˘ +

β˘ δ ψ˘ . hb

(14)

Substituting this into Eq. (8b), we get the coupling condition at x = B for the linear model as follows: ˘ x ψ˘ = h∂x φ hb ∂x φ˘ + β∂ β˘ ˘ x φ˘ = h∂x φ . α∂ ˘ x ψ˘ + β∂ hb

(15a) (15b)

To derive the condition for the nonlinear shallow water model, we use the approximation for the velocity potential (2) again. Since F (z = η) = 0 at the surface, we have ˘ From Eq. (8b), the coupling condiφ = φ˘ and thus δφ = δ φ. tion for a nonlinear model is given by ˘ x ψ˘ . h∂x φ = hb ∂x φ˘ + β∂

(16)

Note that the coupling conditions (15)–(16) are used to transfer the information between the two domains. Coupling condition (15) gives the information of φ˘ and ψ˘ in the simulation area, provided the information of φ from the model area is given. Meanwhile, coupling condition (16) gives the information of φ in the model area, provided the information of φ˘ and ψ˘ from the simulation area is given.

3

Nonlinear shallow water equations

3.1

Characteristic form

We will start with the NSWE in the shore region. Using η = −hb + h and velocity u = ∂x φ, we may rewrite Eq. (10) as follows (starred variables are used here for later convenience): ∂t ? h? + ∂x ? h? u ?

?

?


?

=0 ?

∂t ? u + u ∂x ? u = −g ∂x ?

(17a) −h?b + h?




.

(17b)

The dimensionless form of Eq. (17) for a still water depth h?b = γ ? x ? (where γ ? = tan θ is the beach slope) is obtained by using the scaling factors (Brocchini and Peregrine, 1996) h=

h? u? x? t? ,u = ,x = ,t = , h0 u0 l0 t0

www.nonlin-processes-geophys.net/21/987/2014/

(18)

∂t h + ∂x (hu) = 0 ∂t u + u∂x u = gγ − g∂x h.

(20a) (20b)

The asymptotic solution of this system of equations for wave propagation over sloping bathymetry has been given for several initial-value problems using a hodograph transformation (Carrier and Greenspan, 1958; Synolakis, 1987; Pelinovsky and Mazova, 1992; Carrier et al., 2003; Kâno˘glu, 2004), and also for the boundary-value problem (Antuono and Brocchini, 2007; Li and Raichlen, 2001; Madsen and Schaffer, 2010) that will be used in this article. Since the system is hyperbolic, it has the following characteristic forms: dα dx = 0 on = u−c dt dt dβ dx = 0 on = u + c, dt dt √ in which c = gh, α = 2c − u + gγ t, and β = 2c + u − gγ t.

(21a) (21b)

(22)

Variables α and β are the so-called Riemann invariants, since they do not change their value along the characteristic curves in Eq. (21). Assuming the flow to be subcritical (that is, |u| < c), the first characteristic curves with u − c < 0 are called “incoming”, since they propagate signals towards the shore. The second ones with u + c > 0 are called “outgoing”, since they move towards the deeper waters (carrying information on the wave reflection over the sloping region). 3.2

A trivial solution of the characteristic curve

In the trivial case of no motion (u = η ≡ 0), as well as in the dynamic case presented later, we focus on the incoming characteristic curve. In the rest case, it is given by dx √ = − gγ x. dt

(23)

√ For x 6= 0, substituting y = gγ x results in the general solution for variable y as follows: 1 y = − gγ t + C2 , 2

(24)

Nonlin. Processes Geophys., 21, 987–1005, 2014

6 992

Kristina et al.: Effective Coastal Boundary Conditions f

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry

with a constant C2 . When the curve intersects x = B at time t = τ , with √ h0 the depth at x = B such that h0 = γ B and y(B) = gγ B = c0 , the particular solution is given by 2c0 − gγ (t − τ ) . 2

(25)

In the case of no motion, the boundary data α = α0 (τ ) and β = β0 (τ ) are as follows: α0 = 2c0 + gγ τ, β0 = 2c0 − gγ τ .

(26)

Transforming back to the x variable while using these expressions, we get the incoming characteristic curve x=

1 gγ (2ω − (t − τ ))2 (gγ t − α0 )2 = 4gγ 4

(27)

with ω = c0 /(gγ ). Along this characteristic curve, the Riemann invariant is constant. Figure 2 shows the characteristic curves of the dimensionless NSWE over sloping bathymetry b(x) = −x for x ∈ [0, 1], and LSWE over flat bathymetry h0 = 1, B = 1 for x ∈ [1, 2]. As in our previous paper (Kristina et al., 2012), the405 characteristic curves of the LSWE are given by dx/dt = ±c0 . The “incoming” and “outgoing” characteristic curves are shown by the solid and dashed lines, respectively. For each characteristic curve (27), the location of the shoreline can be determined by looking for the τ = τs for410 which the characteristic reaches the shoreline position, here x = 0, at time t. It is given by the condition ∂x = 0, so that τs = t − 2ω . ∂τ

(28)

As displayed in Fig. 2, the incoming characteristic curves415 that reach the shoreline at time t intersect x = B = 1 at time τ = t − 2 (ω = 1 in this case). Since u equals zero in the rest case, the boundary condition (12) is of course satisfied. 3.3

Boundary value problem (BVP)

Li and Raichlen (2001) and Synolakis (1987) combine linear and nonlinear theory to reduce the difficulties in the assign-420 ment of the boundary data for solving the BVP problem in the NSWE. Later, it is shown that the proper way to solve the assignment problem without using linear theory at all is not given in terms of η or u (Antuono and Brocchini, 2007), but in terms of the incoming Riemann variable α. This article fol-425 lows the approach of Antuono and Brocchini (2010), who use this incoming Riemann variable as boundary data and solve the dimensionless NSWE by direct use of physical variables instead of using the hodograph transformation introduced by Carrier and Greenspan (1958). We do, however, clarify the 430 mathematics of the boundary condition at the shoreline. Given the data of η and u at the seaward boundary x = B, for all time t, we want to find a solution of the NSWE in the sloping region to the shoreline, including the reflected 435

Nonlin. Processes Geophys., 21, 987–1005, 2014

0.5

440

x

y=

0

1 1.5 2 0

of η( trivia B is data α chara along Eq. (

α=α 1

2

3

4

5

t Figure Plotof of the in the motion Fig. 2. 2.Plot the characteristic characteristiccurves curves in case caseofofnono motion 445 (η = u = 0) for the dimensionless NSWE over sloping bathymetry (η = u = 0) for the dimensionless NSWE over sloping bathymetry b(x) = = −x bathymetry h0 = b(x) −x for for xx∈∈[0, [0,1]1]and andLSWE LSWEover overflatflat bathymetry h01,= 1, B = 1 for x ∈ [1, 2]. The “incoming” and “outgoing” characteristic B = 1 for x ∈ [1, 2]. The ”incoming” and ”outgoing” characteristic curves are shown by the solid and dashed lines, respectively. The curves are shown by the solid and dashed lines, respectively. The shoreline x = 0 can be seen as the envelopes of the characteristic shoreline x = 0 can be seen as the envelope of the characteristic curves themselves. curves themselves. 450

waves traveling back into the deeper waters. If the sea is as-

As in our previous paper (Kristina et al., 2012), the sumed in the rest state during the initial condition, the characdata teristic curve of the LSWE are given by dx/dt = ±cpre0 . The are η(B, t) = u(B, t) = 0 for t < 0. In accordance to the “incoming” and “outgoing” characteristic curves are shown vious trivial case, the initial time where a characteristic meets 455 by and dashed respectively. x =the B solid is labeled as τ , andlines we write x = χ (t, τ ), so we have curveτ )+gγ (27), τthe location of the theFor dataeach α = αcharacteristic along the incom0 ≡ 2c(B, τ )−u(B, ing characteristic and β = ≡ 2c(B,for τ)+ u(B, shoreline can be curves determined byβ0looking the τ =τ )τ− s for gγ τ along the outgoing characteristic curves. We can thenhere which the characteristic reaches the shoreline position, (21)t.as xrewrite = 0, atEq. time It is given by the condition β−3α0 ∂x α ==α00 on that χt =u−c= +gγ t, (29a)(28) 460 socurves that τsuch s = t − 2ω. 4 ∂τ

3β0 −α β =displayed β0 on curves such 2, thatthe χt =u+c= +gγ t, (29b) As in Fig. incoming characteristic curves 4

that reach the shoreline at time t, intersect x = B = 1 at time that boundary by the τwhich = t −means 2 (ω = 1 the in this case).values Since are u =carried 0 in the restincase, coming and outgoing characteristic curves. To be concise, we the boundary condition (12) is of course satisfied.

write χt = ∂t χ and χτ = ∂τ χ. Our aim is to obtain a closed 465 equation for the dynamics, and we(BVP) focus on the incoming 3.3 Boundary Value Problem characteristic by fixing α = α0 . We can rewrite Eq. (29a) as follows: Li and Raichlen (2001) and Synolakis (1987) combine linear

and nonlinear theory to reduce the difficulties in the assignβ = 3α + 4 (χt − gγ t) . (30) ment of0 the boundary data for solving the BVP problem in the NSWE. Later, is shown that the proper to solve the Here, β = β(χ , t), itsince we are moving alongway an incoming assignment problem without using linear theory at is not 470 characteristic curve. By taking the total t derivative of all β, we given obtainin terms of η or u (Antuono and Brocchini, 2007) but in   α. This article folterms of the incoming Riemann variable dβ β − 3α0 lows=the and+Brocchini βt approach + βx χt = of βt Antuono + gγ t βx (2010) who use dt incoming Riemann variable 4 as boundary data and solve this = 4 (χtt − gγ ) NSWE , (31) the dimensionless by direct use of physical variables instead of using the hodograph transformation introduced by 475 in which the last equality comes from Eq. (30). In addition, Carrier and Greenspan (1957). We do, however, clarify the the τ derivative of Eq. (30) gives mathematics of the boundary condition at the shoreline. boundary x = B, ∂βGiven the data of η and u at the seaward 3α˙ 0 + 4χ tτ = βx χτ = 3α˙ 0 + 4χtτ ⇒ βx = , (32) for all time t, we want to find a solution of the NSWE in ∂τ χτ the sloping region to the shoreline including the reflected 480 waves traveling back into the deeper waters. If the sea is assumed in www.nonlin-processes-geophys.net/21/987/2014/ the rest state during the initial condition, the data

β=β

whic comi we w close ing c as fol

β=3

Here chara obtai

dβ = dt

in wh the τ

∂β = ∂τ

We s taine βt +

Co tial e

2χτ ( with

χ|t

χτ |τ =

The cond β=− 4c = the ri sisten χtt 6=

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry in which α˙ 0 = ∂τ α0 . We still need an explicit expression for βt , which can be obtained by rewriting Eq. (21b) in the following way:   3β − α0 βt + + gγ t βx = 0 . (33) 4 Combining Eqs. (31)–(33), we get the following differential equation for the incoming characteristic curves:

993

We then make the change of variables ν = −(2ω −t +τ ) and ξ = τ , and Eq. (37) becomes   (1) (1) (1) ν 2ϒνξ − ϒνν − 2ϒν(1) + ϒξ = 0 . (38) Denoting the Fourier transform F(·) with respect to ξ , ρ

(1)

Z∞   (1) (ν, s) = F ϒ (ν, ξ ) (s) = ϒ (ν, ξ ) e−isξ dξ , (39) −∞

2χτ (χtt −gγ ) = (4χtτ +3α˙ 0 ) (gγ t−α0 −χt ) for t > τ, (34a) with boundary conditions (34b) (34c)

χ |t=τ = B χτ |τ =τs = 0 .

The second boundary condition is the shoreline boundary condition. We have 4c = α + β from Eq. (22), which implies β = −α at the shoreline c = 0. Using Eq. (30), we note that 4c = α0 + β = 4(α0 + χt − gγ t) = 0 at the shoreline. Hence, the right-hand side of Eq. (34a) is zero, such that for consistency, χτ must be zero at the shoreline, since generally, χtt 6 = gγ . 3.3.1

we obtain from Eq. (38) a differential equation related to a Bessel equation:   (1) ν 2isρν(1) − ρνν − 2ρν(1) + isρ (1) = 0 , (40) which has the general solution  h i ρ (1) (ν, s) = eisν A1 (s) J0 (sν) − iJ1 (sν)  + A2 (s) [iY0 (sν) + Y1 (sν)] ,

(41)

with J0,1 and Y0,1 the Bessel functions of the first and second kinds. To recover ϒ(ν, ξ ), we just need to take the inverse Fourier transform of Eq. (41), and by using ϒ (1) = χ (1) + να0,1 /2, we get

Perturbation expansion (1)

1 (ν, ξ ) = 2π

Z∞

 h i eis(ν+ξ ) A1 (s) J0 (sν) − iJ1 (sν)

Due to the nonlinearity in χ , we use a perturbation method to solve Eq. (34). We expand it in a perturbation series around the rest solution (27), with the assumption of small data at x = B. Using the linearity ratio  = A/ h0 (A is the wave amplitude), we say a wave is small if   1, and expand as follows:

with ξ = τ ≤ t.

α0 = α0,0 + α0,1 + O( 2 ),

(35a)

3.3.2

χ = χ (0) + χ (1) + O( 2 ),

(35b)

2

(35c)

In order to calculate the unknown functions A1 (s) and A2 (s), we need to assign the boundary conditions (34). In (ν, ξ ) space, t = τ corresponds to ν = −2ω, and by imposing the first boundary condition, we get
−F α0,1 ωe2isω = A1 (s) [J0 (2sω) + iJ1 (2sω)] + A2 (s) [iY0 (2sω) − Y1 (2sω)] . (43)

τs = τ0 (t) + τ1 (t) + O( ),

in which α0,0 = 2c0 +gγ τ is the incoming Riemann invariant in case of no motion, χ (0) is given by Eq. (27), and τ0 = t − 2ω. By substituting Eq. (35) into Eq. (34), we obtain at first order in :     (1) (1) (1) (2ω − t + τ ) χtt + 2χtτ − χτ(1) − χt − α0,1 3 + (2ω − t + τ ) α˙ 0,1 = 0, 2 (1) χt=τ = 0,

(36b)

(1) χτ(0) τ (t, τ0 ) τ1 + χτ (t, τ0 ) = 0.

(36c)

(36a)

χ

−∞

h i ν + A2 (s) iY0 (sν) + Y1 (sν) ds − α0,1 , 2

(42)

Boundary value assignment

The second boundary condition is given by Eq. (36c), in which (1)

χτ(1) = −χν(1) + χξ i = 2π

Z∞

 h J1 (sν) i eis(ν+ξ ) A1 (s) sJ0 (sν)−isJ1 (sν)− ν

−∞

By letting ϒ (1) equal χ (1) − (2ω − t + τ )α0,1 /2, we can rewrite Eq. (36a) as (1)

(1)

(1)

(2ω − t + τ ) (ϒtt + 2ϒtτ ) − ϒτ(1) + ϒt

= 0.

www.nonlin-processes-geophys.net/21/987/2014/

(37)

h Y1 (sν) i + A2 (s) isY0 (sν) + sY1 (sν) + ds iν α0,1 ν α˙ 0,1 + − , 2 2

(44)

Nonlin. Processes Geophys., 21, 987–1005, 2014

994

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry

is evaluated at τ = τ0 ; i.e., ν = 0 needs to be finite. Evaluating Eq. (44) at ν = 0 gives us convergence when the coefficient A2 (s) is zero, which avoids an unbounded result. Hence, from the first boundary condition (36b), coefficient A1 (s) is given by F(α0,1 )ωe2isω A1 (s) = − . J0 (2sω) + iJ1 (2sω)

of Eq. (49) in Eq. (6) yields

0 =δ

0

1 1 − Bkl ψk φl − Akl ψk ψl − Gkl ψk ψl 2 2 ZB   i 1 1 + φ∂t η − gη2 − h (∂x φ)2 dx dt 2 2

(45)

The solution of the incoming characteristic curves at first order is thus given by

Z∞ −∞


J0 (sν) − iJ1 (sν) ds J0 (2sω) + iJ1 (2sω) (46)

The shoreline position must satisfy χτ |τ =τs = 0, and in the first-order approximation, it is given by h i xs (t) =χ (0) (t, τ0 ) +  χτ(0) (t, τ0 ) τ1 + χ (1) (t, τ0 )   + O  2 . (47) Since τ = τ0 corresponds to ν = 0 and ξ = t − 2ω, we get  
ω −1 xs (t) = −F F α0,1 . (48) J0 (2sω) + iJ1 (2sω)

4.1

Finite element implementation

The region x ∈ [B, L] will be approximated using a classical Galerkin finite element expansion. We use first-order spline polynomials on N elements with j = 1, . . . , N +1 nodes. The variational structure is simply preserved by substituting the expansions φ˘ h (x, t) = φj (t)ϕj (x), ψ˘ h (x, t) = ψj (t)ϕj (x) , and η˘ h (x, t) = ηj (t)ϕj (x) (49a) into Eq. (6) for x ∈ [B, L] concerning N elements and (N + 1) basis functions ϕj . We used the Einstein summation convention for repeated indices. To ensure continuity and a unique determination, we employ Eq. (13) and substitute β˜ ψ1 (t)ϕ1 (x) and hb

η(x, t) = η(x, ˜ t) + η1 (t)ϕ1 (x),

(49b)

with ϕ1 the basis function in element 0 for x ∈ [xs , B], and ˜ with φ(B, t) = η(B, ˜ t) = 0. For linear polynomials, the use Nonlin. Processes Geophys., 21, 987–1005, 2014

ZB 

+


1 ∂t η + ∂x (h∂x φ) δ φ˜ − (∂t φ + gη + ∂x2 φ)δ η˜ dx 2

xs

+ (φδη) |x=xs +

ZB 

dxs − (φ∂t η) |x=xs δxs dt

 1 (∂t η + ∂x (h∂x φ))ϕ1 δφ1 − (∂t φ + gη + ∂x2 φ)ϕ1 δη1 dx 2

xs

− h∂x φ|x=B δφ1 −

i β˜ h∂x φ|x=B δψ1 dt, hb

(50b)

where we introduced mass and stiffness matrices Mkl , Skl , Akl , Bkl , and Gkl , and used endpoint conditions δηk (0) = ˜ δηk (T ) = 0, connection conditions δ η(B, ˜ t) = δ φ(B, t) = ˜ δ ψ(B, t) = 0, and no-normal through-flow conditions at x = L. The matrices in Eq. (50) are defined as follows:

Effective boundary condition

˜ φ(x, t) = φ(x, t) + φ1 (t)ϕ1 (x) +

0

− (Akl ψl +Bkl φl + Gkl ψl ) δψk

ν − α0,1 . 2

4

h (Mkl η˙ l −Skl φl −Bkl ψl )δφk −(Mkl φ˙ k +gMkl ηk )δηl

= eis(ν+ξ +2ω) ωF α0,1

(50a)

xs

ZT

χ (1) (ν, ξ ) 1 =− 2π

ZT h 1 1 Mkl φk η˙ l − gMkl ηk ηl − Skl φk φl 2 2

ZL

ZL

h∂x ϕk ∂x ϕl dx,

ϕk ϕl dx, Skl =

Mkl =

B

B

ZL Akl =

ZL α∂ ˘ x ϕk ∂x ϕl dx, Bkl =

B

˘ x ϕk ∂x ϕl dx, β∂

B

ZL and Gkl =

γ˘ ϕk ϕl dx.

(51)

B

Provided we let the size of the zeroth element go to zero such that the underlined terms in Eq. (50b) vanish, the equations arising from Eq. (50) are Mkl η˙ l − Skl φl − Bkl ψl − δk1 (h∂x φ) |x=B − = 0 Mkl φ˙ k + gMkl ηk = 0  β˜  Akl ψl + Bkl φl + Gkl ψl − δk1 h∂x φ |x=B − = 0, hb

(52a) (52b) (52c)

www.nonlin-processes-geophys.net/21/987/2014/

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry with Kronecker delta symbol δkl (one when k = l, and zero otherwise) and Eq. (10) for x ∈ [xs , B] with boundary condition (12). Taking this limit does not jeopardize the time step, as this zeroth element lies in the continuum region, in which the resolution is infinite. The time integration is solved using ode45 in MATLAB, which uses its internal time step. From Eq. (52), we note that we need the depth h and the velocity u from the nonlinear model at x = B, whose values are given at time t = τ in the characteristic space. The definitions (22), while using α = α0 and β in Eq. (30) with expansions up to first order, yield h = c2 /g =

1 (α0 + β)2 16g (0)

995

solve for the perturbation expansion in the nonlinear area, but we do not perturb the incoming wave data. The values η˘ and u˘ in Eq. (56) are obtained from the simulation area [B, L]. In this region, we only have the values of ˘ and ψ. ˘ The depth-averaged velocity u(B + , t) is deterη, ˘ φ, mined by using the approximation (13) as follows: u˘ = ∂x φ˘ +

β˘ ∂x ψ˘ at x = B + , hb

which is the limit from the right at node 1. The solutions of η = h − hb and u in Eq. (55) account for the reflected wave, so we may define η = ηI + ηR and u = uI + uR ,

(1)

2

= α0,0 + χt − gγ t +  α0,1 + χt /g gγ t−α0,0 (1)

2 = α0,0 + −gγ t+ α0,1 +χt /g 2 1 (1)

2 /g (53a) = c0 + gγ (τ − t) +  α0,1 + χt 2
1 (1)
(53b) u = gγ t + β − α0 =  α0,1 + 2χt . 2 Note that for  = 0, we indeed find the rest depth hb (x) = (1) γ x. The function χt follows from evaluation of Eq. (46), and since t = τ is equivalent to ν = −2ω, we immediately obtain (1)

χt |t=τ ≡ χν(1) (−2ω, ξ ) Z∞
i J1 (2sω) =− eisξ F α0,1 ds 4π J0 (2sω) + iJ1 (2sω) −∞

α0,1 . − 2

(54)

The solutions of h and u at t = τ are thus given as follows: h(B, t) = hb + η  
c2 c0 J0 (2sω) = 0 +  F −1 F α0,1 g g J0 (2sω) + iJ1 (2sω)  
iJ1 (2sω) . u(B, t) = −F −1 F α0,1 J0 (2sω) + iJ1 (2sω)

(55a) (55b)

In order to calculate the solutions for h and u at x = B and the shoreline position, we need the data of the incoming Riemann invariants at the first order as follows: α0,1 ≈ α − α0,0  p p ˘ − gγ B |x=B + − u| ˘ x=B + , = 2 g(γ B + η)

which is obtained by disregarding higher-order terms in Eq. (35a). This expression is actually the incoming Riemann invariant in LSWE (Kristina et al., 2012). By imposing the effective boundary condition (EBC) and choosing the location x = B before the nonlinearity arises, we thus do actually www.nonlin-processes-geophys.net/21/987/2014/

(58)

where ηI and ηR are the wave elevations of the incoming and reflected waves, respectively, at x = B. This superposition is also described in Antuono and Brocchini (2007, 2010), and is actually in line with our EBC concept, since linearity holds in the simulation area. To obtain the expression for the reflected wave, we need to know the incoming one. Using the knowledge of the incoming and outgoing Riemann invariants in the LSWE as derived in Kristina et al. (2012), the observation operator is given by O = hu˘ + cη˘ = 2cηI ,

(59)

which is calculated using approximation (57). We can thus calculate the incoming wave elevation for any given wave signal at x = B. Implementation of this observation operator allows us to have any initial waveform at the point of tsunami generation, and to let it travel over the real bathymetry to the seaward boundary point x = B. From Eq. (55), the expressions for the reflected wave are as follows: ηR = M(ηI ) h i c0 J0 (2sω) = F −1 F(α0,1 ) − ηI g J0 (2sω) + iJ1 (2sω) uR = M(uI ) h
= −F −1 F α0,1

i iJ1 (2sω) − uI , J0 (2sω) + iJ1 (2sω)

(60a)

(60b)

where the Fourier transform and its inverse for any incoming wave signal are evaluated using the FFT and IFFT functions in MATLAB. The influxing operator is defined as the coupling condition in Eq. (52) to send the NSWE result to the simulation area. It is shown that we need the value of h∂x φ, and hence I = h∂x φ = (hb + η) u.

(56)

(57)

(61)

In order to verify the EBC implementation, we perform numerical simulations with a code that couples the LSWE in the simulation area to the NSWE in the model area (Bokhove, 2005; Klaver, 2009). For numerical simulation of the LSWE, we use a finite element method, while for the NSWE, we use a finite volume method. The implementation of the finite volume method is explained in Appendix A. Nonlin. Processes Geophys., 21, 987–1005, 2014

10 996

Kristina et al.: Effective Coastal Boundary Conditions forfor Tsunami over sloping Slopingbathymetry Bathymetry W. Kristina et al.: EBC tsunamiWave waveRun-Up run-up over

Solitary 55.1Study case wave

745

750

755

760

765

770

775

Nonlin. Processes Geophys., 21, 987–1005, 2014

η(x,t)

0.98

(a)

0

5

0

5

0

5

0

5

0

5

x

10

15

20

10

15

20

10

15

20

10

15

20

10

15

20

η(x,t)

1.04 1.02 1 0.98

(b)

x

1.08

η(x,t)

1.06 1.04 1.02 1 0.98

(c)

η(x,t)

740

1

1.1 1.08 1.06 1.04 1.02 1 0.98

(d)

x

x

1.04 1.02 η(x,t)

735

1.02

1 0.98 0.96

(e)

x

Fig. 3.3.Run-up waves Figure Run-upofofa solitary a solitary waveover overa acanonical canonicalbathymetry bathymetryat t= 40,40, (c)(c) t =t = 50,50, (d)(d) t =t 60, (d)and t =(e) 70.t The attimes times(a) (a)t = t =30, 30,(b)(b) t= = 60, = solid is LSWE EBCwith implementation at x = 10,atthe 70. Theline solid line iswith LSWE, EBC implementation x dashed = 10, anddashed dotted-dashed lines are lines the coupling of LSWE NSWE the and dotted-dashed are the couplings of and the LSWE model respectively, and the symbols laboratory data of Synoand NSWE models, respectively, and theare symbols are the laboratory lakisof(1987). data Synolakis (1987). 2 1.5

not appear. Figure 3 shows the time evolution of this profile 1 for scaled time t = 30–70 with 10 increments. It can be seen that the0.5EBC implementation and the coupling of numerical solutions0 agree well with the laboratory data. The comparison of−0.5 the shoreline movement between the simulation with EBC implementation and the coupling of numerical solutions −1 60 till the 80 100 is shown 0in Fig. 4.20For the 40 simulation scaled physical t time t = 100, the computational time for the coupled numerFig.solutions 4. The shoreline a solitary wave ical in both movement domains isof 10.9 s, while theintroduced computa-in Synolakis LSWE model to the NSWE is shown tional time(1987). of theThe simulation withcoupled the EBC implementation by the dashed line, 18 while solid time. one is the shoreline movement only takes about % the of that Hence, we notice thatof thesimulation LSWE model with EBC reduces implementation. the with theanEBC the computational time xs (t)

730

The run-up of a solitary wave is studied byone means the wellThree test cases are considered. The first is aofsynthetic known case of Synolakis (1987). A solitary centereditat one concerning a solitary wave, such that wewave can compare x = other x0 at results. t = 0 has the followingwe surface profile: with Subsequently, consider periodic wave influx as the second case to test the robustness of the tech2 η(x, 0) = Athere sechis κ(x − x0 ). interaction between the (63) nique when continuous incoming and reflected waves. The third case is a more realistic concerningthe tsunami propagation and and run-up on Weone benchmark EBC implementation thebased coupling simplified bathymetry the Aceh coastline.data of Synolakis of numerical solutionsat with experimental The location EBC point determined the (1987) providedof attheNOAA Centeris for Tsunami from Research linearity condition  = A0 / h0  1. FromSolitary linear theory, the (http://nctr.pmel.noaa.gov/benchmark/). wave runwave amplification over depth is given by the ratio A up over a canonical bathymetry is considered with the scaled 0= p √ 4 h/ h0 , where A and hand areκthe wave amplitude Aamplitude 3A/4 = 0.1178. Theand iniA = 0.0185 = initial depth. Hence, the EBC point must be located at depth tial condition is centered at x = 37.35 over the bathymetry 0 with q a constant slope γ = 1/19.85 for x < 19.85. The EBC 5 4 h/ 4 . hpoint located at x = 10 such that the domain is divided(62) into 0  isA the model area for x ∈ [−5, 10] and the simulation area for Since a dispersive model is also used in the simulation x ∈ [10, 80]. The spatial grid size is ∆x = 0.25 in the simulaarea, we will discuss the dispersion effect at this EBC point tion area and ∆x = 0.0125 in the model area for the numerias well. The non-dispersive condition is given by µ2 = cal solution of NSWE. In all cases, several spatial resolutions (k0 h0 )2  1, where k0 = 2π/λ0 is the wavenumber and λ have been applied to verify numerical convergence. For the is the wavelength. In linear wave theory, the wavelength detime integration, we use the fourth order ode45 solver that creases with the square root of√the depth when running in uses its own time step in MATLAB. shallower water, that is λ0 = λ h0 / h. Using this relation, The simulations with the EBC implementation and the we can thus investigate the significance of the dispersion coupling of numerical solutions are only presented for the given the information of the initial condition and bathymetry LSWE model in the simulation area since the initial condiprofile. tion has long wavelength and thus dispersion effects will not appear. Figurewave 3 shows the time evolution of this profile for 5.1 Solitary scaled time t = 30 − 70 with 10 increments. It can be seen thatrun-up the EBC the by coupling The of aimplementation solitary wave is and studied means of of numerical the wellsolutions well with(1987). the laboratory compariknown caseagree of Synolakis A solitarydata. waveThe centered at the shoreline movement between the simulation with xson = xof at t = 0 has the following surface profile: 0 EBC implementation and the coupling of numerical solutions η(x, 0) = A sech24.κ(x x0 ).simulation till the scaled physical (63) is shown in Fig. For−the timebenchmark t = 100, the time for the numerWe thecomputational EBC implementation andcoupled the coupling solutionssolutions in both domains is 10.9 s. While computaofical numerical to experimental data oftheSynolakis tional provided time of the simulation with theforEBC implementation (1987) at the NOAA Center Tsunami Research only takes about 18 % of that time. Hence, we wave noticerunthat (http://nctr.pmel.noaa.gov/benchmark/). Solitary theover simulation withbathymetry the EBC reduces the computational time up a canonical is√considered with the scaled significantly, to approximately 82 %,=compared withinithe amplitude A =up 0.0185 and κ = 3A/4 0.1178. The computational in theatwhole domain. tial condition is time centered x0 = 37.35 over the bathymetry, to slope show γthe dispersion we consider withIna order constant = 1/19.85 for effect, x < 19.85. The EBCa shorter wave with given Eq.domain (63) forisκdivided = 0.04, point is located at xthe = profile 10, such thatinthe x0 = m, and m. The bathymetry is given by into the150 model area A for=x0.1 ∈ [−5, 10] and the simulation area constant depth for x >grid 50 m, continued a constant for x ∈ [10, 80]. 10m The spatial sizes are 1x =by0.25 in the slope γ =area 1/5 and towards shore. spatial grid simulation 1x =the 0.0125 in A theuniform model area for the ∆x = 1msolution is used of in NSWE. the simulation area several and ∆xspatial = 0.015 numerical In all cases, res-m in the model area applied for the numerical solution of the NSWE. olutions have been to verify numerical convergence. Evaluating (62) for we  =use 0.02 1, the EBC ode45 point must be For the time Eq. integration, thefourth-order solver located 3.3 m. Accordingly, we choose this seaward 0  time that uses at itshown step in MATLAB. boundary point at hwith at theimplementation toe of the slope, 0 = 10 The simulations themEBC andthat theis at x = Bof=numerical 50 m. Therefore, weare divide domain for intothe the coupling solutions onlythe presented simulation for simulation x ∈ [50, 250] m and for LSWE modelarea in the area, sincethe themodel initial area condix ∈has [−5,a 50] m.wavelength, and thus dispersion effects will tion long

1.04

www.nonlin-processes-geophys.net/21/987/2014/

2

η(x,t) [m]

η(x,t) [m]

0.05 0 −0.05

1

0

(a)

0

100 150 x [m]

200

0.1 0 20

40

x [m]

60

80

60

t

80

100

η(x,t) [m]

40

η(x,t) [m]

20

0 −0.05

0

50

100 150 x [m]

0

20

40

0

50

100

200

250

60

80

100

150

200

250

0.05 0

(b)

−0.05 x [m]

0.1

0.05

0

50

100 150 x [m]

200

0.05 0 −0.05

250

x [m] Figure Theshoreline shorelinemovement movement of aa solitary in in (c) Fig. 4. 4.The solitarywave waveintroduced introduced 0.1 0.1 Synolakis toto thethe NSWE is shown Synolakis(1987). (1987).The TheLSWE LSWEmodel modelcoupled coupled NSWE is shown 0.05 0.05 by movement of of by the the dashed dashedline, line,while whilethe thesolid solidone oneisisthe theshoreline shoreline movement 0 0 the with implementation. −0.05 −0.05 the LSWE LSWE model withananEBC EBC implementation. Kristina etmodel al.: Effective Coastal Boundary Conditions for Tsunami Wave Run-Up Sloping Bathymetry 0 50 100 over 150 200 250 0 50 100 150 11 200 x [m] x [m] (d) η(x,t) [m]

η(x,t) [m]

0.1

xs(t) [m]

0.05 0

Kristina et al.: Effective Coastal Boundary Conditions for Tsunami Wave Run-Up over Sloping Bathymetry 0

50

100

x [m]

0.15

150

200

11

0.5

250

1 −1.5 In Fig. 5, we can see the initial profile of the solitary wave.

x (t) [m]

η(x,t) [m]

0 10 20 30 40 Fig. 5. A5. solitary wave initial condition for for thethe NSWE (dottedFigure A solitary wave initial condition NSWE (dotted- (a) t [s] 780 Comparisons between these two simulations at several time 805 0.1 −1 dashed (dashed line), andand thethe linear dashedline) line)coupled coupledtotothe thelinear linearmodel model (dashed line), linear −1.5 steps can be seen in Fig. 6 (left: LSWE, right: LVBM). Commodel with the EBC implementation (solid line) at x = 50 m. The model with the EBC implementation (solid line) at x = 50 m. The 0.05 −0.5 paring the left and right figures, we can see that the wave solid solidand anddashed dashedlines linesare areonontop topofofanother. one another. −1

is slightly 0dispersed in the LVBM. Because we have flat −0.5 bathymethy in this case, the dispersion ratio at the simula0.1−0.05 0.5 0 50 100 150 200 250 785 tion area is constant and given by µ2 = 0.39 < 1. Hence, it is 810 0.05 0.05x [m] 0 significantly, down by approximately 82 %, compared with 0 0 shown that1 the long waves propagate faster than the shorter the −0.05 computational time in the whole domain. −0.05 0 10 20 30 40 0.5 0 50 for 100 the 150NSWE 200 250 0 50 solitary 100 150 200 initial 250 Fig. 5. A wave condition (dottedones x [m] [s] 7, the shoreline move(a) In order to show (a) in LVBM simulations. In tFig. the dispersion effect,x [m]we consider a dashed line) coupled to the linear model (dashed line), and the linear ment caused by this solitary wave is shown. The paths of 0.2 0.05 1 −1.5 shorter wave with the profile given in Eq. (63) for κ = 0.04, 0 10 20 30 40 model with the EBC implementation 0.1 0 (solid line) at x = 50 m. The t [s]the shoreline are also presented (b) characteristic curves forming x0solid = 0150 m, and A = 0.1 m. The bathymetry is given by a −0.05another. and dashed lines are on top of −1 figure. We can see that the shoreline is formed by 815 0 20 40 of 6010 m 80 for 100 x > 50 0 m,20continued 40 60 80 constant depth by a100con- 790Fig.in7.this x [m] x [m] (b) The shoreline movement of a solitary wave for the linear the (a: envelope of the characteristic curves. The result with the stant 0.1slope γ = 1/5 towards the0.1 shore. A uniform spatial model −0.5 b: LVBM) coupled to the NSWE (dashed line) LSWE, 0.05 0.05 0.1 LVBM shows a lower run-up but higher (solid run-down with some grid 1x = 1 m is used in the simulation area, and 1x equals and the linear model with an EBC implementation line). Paths 0 0 0.05 0.05 0characteristic oscillations at later times. of the first-order curves are shown by the thin lines. 0.015 m in the model area for the numerical solution of the −0.05 −0.05 0 0 50 100 150 200 250 0 100 150 200 250 0 50 x [m] x [m] (c) For simulation till the physical time t = 40 s, the comNSWE. Evaluating Eq. (62) for  0.1=−0.05 0.02 50 1, 100 the EBC point 0.1−0.05 0 0.5 0 150 200 250 50 100 150 200 250 x [m] (a)0.05be located at xh[m]  3.3 m. Accordingly, 795 putational time for the coupled numerical solutions in both 820 0.05 must we choose this 0 0.05 0 0.2 domains is103.3 times 10 the physical seaward boundary point at h0 = 100 m at the toe of the slope, 20 time for30the LSWE 40and −0.05 0.1 −0.05 0 0 50 100 150 200 250 0 50 100 150 200 250 t [s] (b) that is at x = B = 50 m. We therefore divide the domain into x [m] x [m] (d) 0 −0.05 0 20 area 40 60 x 80 100 250] m, 0 20 40 model 60 80 area 100 the(b)simulation for ∈ [50, and the x [m] x [m] Fig. 7. movement of ofaasolitary solitarywave waveforforthethe linear Figure 7. The The shoreline movement linear Fig. 6. Free-surface profiles of solitary wave propagation are shown 0.1 for x ∈ 0.1 [−5, 50] m. model (a: LSWE, b: LVBM) coupled to the NSWE (dashed line) 2.2 times the physical time for the LVBM. While the commodel (a: LSWE; b: LVBM) coupled to the NSWE (dashed line), for the coupled linear model (left: LSWE, right: LVBM) with the 0.05 0.05 In Fig. 5, we can see the initial profile of the solitary wave. putational and the linear model with an EBC implementation (solid line). Paths NSWE (dashed and dotted-dashed lines), and for the linear model and the linear model with an EBC implementation (solid line). Paths time of the simulation using an EBC only takes 0 0 Comparisons between these(solid two line), simulations at several time 0.12 first-order characteristic curves shown by thin lines. −0.05 at times (a) t = 10 s, with an−0.05 EBC implementation of of thethe first-order characteristic curves areare shown thethe thin lines. times the physical time for the LSWE andby0.06 times 0 50 100 150 200 250 0 50 100 150 200 250 x [m] x [m] (c) steps can be seen in Fig. 6 (left: LSWE; right: Com(b) t = 200.1s, (c) t = 30 s, (d) t = 40 s. The0.1 solid andLVBM). dashed lines 800 for the LVBM. Hence, we notice that the simulation with are on top of another. paring the left and right figures, we0.05can see that the wave the EBC reduces the computational time significantly, up to 0.05 0 is slightly dispersed in the LVBM. 0Because we have flat approximately 97 %, compared with the computational time −0.05 −0.05 0 100 case, 150 200 250 0 50 100 150 simula200 250 bathymetry in50 this the dispersion ratio inx [m]the the shoreline are also presented forcharacteristic the numericalcurves modelsforming in the entire domain. The computax [m] (d) 2 = 0.39 < 1. Hence, it is tion area is constant and given by µ in this figure. We can see that the shoreline is formed tional time for the LSWE with an EBC is slower than the one by In Fig. 5, we can see the initial profile of the solitary wave. Fig. 6. Free-surface profiles of solitary wave propagation are shown shown that the long waves propagate faster than the shorter theLVBM envelope characteristic curves. time The step resultofwith Comparisons between these two simulations at several time 805 with and of an the EBC, because the internal the the 2.2 times the physical time for the LVBM. While the comfor the coupled linear model (left: LSWE, right: LVBM) with the ones in LVBM simulations. In Fig. 7, the shoreline moveLVBM shows a lower run-up but higher run-down, with some steps can be seen in Fig. 6 (left: LSWE, right: LVBM). Comode45 time step routine in MATLAB required a smaller time NSWE and dotted-dashed lines), and the putational the simulation usingthe anstability. EBC only takes ment caused this solitary wave is shown. paths of steposcillations attime later times. paring the(dashed left by and right figures, we can see for thatThe thelinear wavemodel ∆t (compared to of the LVBM) to preserve an EBC implementation (solid line), at times (a) t = 10 s, 0.12 times the physicaloftime for thecompare LSWE well and with 0.06 times is with slightly dispersed in the LVBM. Because we have flat The shoreline movement our result (b) t = 20 s, (c) t = 30 s, (d) t = 40 s. The solid and dashed lines 800 for the LVBM. Hence, we notice that the simulation with bathymethy in this case, the dispersion ratio at the simulathe one of Choi et al. (2011). We can see the comparison in are on top of another. www.nonlin-processes-geophys.net/21/987/2014/ Nonlin. Processes Geophys., 21, significantly, 987–1005, 2014 theThe EBC reduces time tion area is constant and given by µ2 = 0.39 < 1. Hence, it is 810 Fig. 8. solution of the Choicomputational et al. (2011) gives a higher pre- up to approximately 97 %, compared with the computational diction for the shoreline, but it cannot follow the subsequent time shown that the long waves propagate faster than the shorter the numerical the entire domain. The smallfor positive wave. It models may bein caused by neglecting the computareones in LVBM simulations. In Fig. 7, the shoreline movetional time for the LSWE with an EBC is slower than In Fig. 5, we can see the initial profile of the solitary wave. flection wave and nonlinear effects in their formulation. Wethe one ment caused by this solitary wave is shown. The paths of s

xs(t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

s

x (t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

785

η(x,t) [m]

η(x,t) [m]

780

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

0

Fig. 7. model ( and the of the fi

250

Fig. 6. Free-surface profiles of solitary wave propagation are shown Figure −1.5 6. Free-surface profiles of solitary wave propagation are for the for coupled linear model right: LVBM) with the shown the coupled linear (left: modelLSWE, (left: LSWE; right: LVBM) NSWE dotted-dashed lines), and for and the linear −1 (dashed with the NSWE and (dashed and dotted-dashed lines), for themodel linwithmodel an EBC implementation (solid line), at times (a) t(a) =t10 ear with an EBC implementation (solid line), at times = s, (b)s,−0.5 t(b) = 20 t =t 30 s, (d) t =(d) 40t s. and and dashed lines 800 10 t =s,20(c)s, (c) = 30 s, and = The 40 s.solid The solid dashed are onare another. lines onoftop of one another. 0top

0.15

−

0

100

0.1

−1 0

0.05

−0.05

250

0.2

0

(b)

−0.5

50

η(x,t) [m]

η(x,t) [m]

0.5

−0.05

−

0.1

1.5

xs (t)

onsider a κ = 0.04, given by constant atial grid 0.015 m NSWE. t must be seaward e, that is into the area for

(a)

xs(t) [m]

times (a) t = 30, (b) t = 40, (c) t = 50, (d) t = 60, (d) t = 70. The solid line is LSWE with EBC implementation at x = 10, the dashed Fig. 5. A solitary wave initial condition for the NSWE (dottedand dotted-dashed lines are the coupling of LSWE and NSWE dashed line) coupled to the linear model (dashed line), and the linear model respectively, and the symbols are laboratory data of Syno- model with the EBC implementation (solid line) at x = 50 m. The and dashed lines are on top of another. lakisKristina (1987). et al.: EBC for tsunami wave run-up over sloping solid W. bathymetry 997

η(x,t) [m]

d numercomputamentation otice that onal time with the

2.2 tim putatio 0.12 tim for the the EB approx for the tional t with LV ode45 t step ∆t The the one Fig. 8. diction small p flection also co Fig. 9. tion at accurac the enti for the become coastlin

20

30

t [s]

40

50

η(x,t) [m]

10

Fig. 8. Comparison of the shoreline movement of Choi et al. (2011) (c) (dashed line) the LSWE withRun-Up an EBC simulation (solid line) for Kristina et al.: Effective Coastal Boundary Conditions for and Tsunami Wave over Sloping Bathymetry wavetsunami case. wave run-up over sloping bathymetry W. Kristina et al.:solitary EBC for

12 998

0.05

−0.05

0.5

−0.05 500

−0.05

0

50 0

0 −0.05

0.05 0 −0.05

10050 150100 200150 250 200 x [m] x [m]

0.05 0 −0.05

0 50 100 150 200 250 0 −0.1 50 100 150 200 250 −0.1 x [m] x [m] (d)250 0 50 100 150 200 0 50 100 150 200 x [m] x [m] (b) Fig. 9. Free-surface profiles of solitary wavewave propagation are shown Figure 9. 0.05 Free-surface profiles of solitary propagation are 0.05 −0.5 0 10 20 30 40 50 0 0 for the for coupled LSWELSWE with the NSWE (dashed and dotted-dashed shown the coupled with the NSWE (dashed and dottedt [s] −0.05 −0.05 lines), the an EBC implementation (solid line), and dashedfor lines), for the with LSWE with an EBC implementation (solid −0.1 −0.1 LSWE −0.15 −0.15 Figure 8. Comparison of the shoreline movement of Choi et (2011) al. for theand LSWE the100method of (2011) Fig. 8. Comparison of the shoreline movement of Choi et al. line), for0 with the 50LSWE with the method et al. (2011) 150 200 Choi 250 et al.of 0Choi50(solid 100 line 150 with 200 x [m] x [m] (c) (2011) (dashed line) the LSWE with an simulation EBC simulation marker) at times (a)marker) t = 10 s,at(b) t =(a) 20ts,=(c) 30ts,=(d) t= (solid line with the “o” times 10ts,=(b) 20 s, (dashed line) and theand LSWE with an EBC (solid(solid line) for ’o’ line) for wave the solitary (c) s.t = 30solid s,0.2andand (d)dashed t = 40 s. Theare solid lines are on top 40 The lines on and top dashed of0.2another. solitary case. wave case.

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

825

830

845

250

0.1

50

100 150 x [m]

0

250 0.05

50

100 150 x [m]

200

250

η(x,t) [m]

η(x,t) [m]

200

0 0 Fig. 10. Free-surface profiles of periodic waves are shown for the

η(x,t) [m]

Using the same bathymetry profile in the previous case, we −0.05 −0.05 coupled linear model (left: LSWE, right: LVBM) with the NSWE (a) influx a0periodic right boundary (x L) with the 50 100wave 150 at 200the250 0 50 100 =150 200 250 x [m] (a) (dashed and dotted-dashed lines), and for thex [m]linear model with an profile: 0.05 EBC implementation (solid line),0.05at times (a) t = 20 s, (b) t = 40 s, 0 (c)0 t = 60 s, (d) t = 75 s. The solid and dashed lines are on top of −0.05

−0.05 η(L,another t) = A atsin (2πt/T ) −0.1 several plots. 0

(b)

50

100 150 x [m]

200

(64)

−0.1 250

0

50

100 150 x [m]

200

250

η(x,t) [m]

xs(t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

xs(t) [m]

xs(t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

in which A = 0.05 m is the amplitude and period T = 20 s. 0.05 0.05 0 0 A smoothened characteristic function until t = 10 s is used 2 −0.05 −0.05 in influxing this periodic wave. We −0.1 use uniform spatial grid −0.1 −0.15 −0.15 0 m in 50 1 100 simulation 150 200 250area and 0 50 = 1000.015 150 m 200 in250 ∆x = 1 the ∆x the x [m] x [m] (c) model0.2 area for the numerical solution of the NSWE. 0.2 0 (b) As0.1in the previous case, we 0.1also choose the seaward 0 0 boundary point −1 at h0 = 10 m at the toe of the slope, that Fig. 11. −0.1 −0.1 0 50 100 150 is 200 0 50 100 m. 150 Thus, 200 250 is the simulation area for 250 x∈ x [m] x [m] (d)at x = B = 50 model ( −2 [50, 250] m and the model area for x ∈ [−5, 50] m. Comparand for Fig. 10. Free-surface profiles ofofperiodic waves aretime shown forthe the Figure 10. Free-surface profiles periodicatwaves are shown for isons between simulations several steps can Paths of −3these two 20 30 LSWE; 40 right: 50 60withthe 80 coupled linear model (left: LSWE, right: LVBM) the70 NSWE coupled linear model (left: LVBM) with NSWE be seen in Fig. 10 (left: LSWE, right: LVBM). We can see lines. t [s] (a) and dotted-dashed lines), and for the linear model with an (dashed and dotted-dashed lines), and is forslightly the lineardispersed model with in the comparison that the wave inanthe 2 EBC implementation implementation (solidline), line),atattimes times(a) (a)t t==2020s,s,(b) (b)t = t = 40 (solid s, s, LVBM. The dispersion ratiosolid at the simulation area is 40 given (c) tt = = 60 s, (d) t = 75 s. The and dashed lines are on top 60 s, and (d) t = 75 s. The solid and dashed lines are on topof robustl by µ2 =at0.0986 < 1, which is less dispersive than the previanother several plots. of one another in1 several plots. and the ous case. In Fig. 11, the shoreline movement caused by the For periodic wave0 as well as the paths of characteristic curves putatio forming the shoreline are shown. Observing the results of −1 in the simulation area and 1x = 0.015 m in grid 1x =2 1 m this case, we can that the EBCoftechnique can deal 850 domain the model area for conclude the numerical solution the NSWE. 1 −2 As in the previous case, we also choose the seaward boundary 0point −3 at h0 = 10 m at the toe of the slope, that 20 30 40 50 60 70 80 is, at x = B area is thus for x ∈ (b) −1 = 50 m. The simulationt [s] As in the previous case, we also choose the seaward 5.2 Periodic wave boundary point at h0 = 10 m at the toe of the slope, that [50, 250] m and the model area is for x ∈ [−5, 50] m. Com5.2 Periodic wave Fig. −2 11. The shoreline of at periodic for the linear between these twomovement simulations severalwaves time steps is at xthe =B = 50 m. Thus,profile the simulation area is case, for x ∈ parisons Using same bathymetry as in the previous (a: in LSWE, b: LVBM) coupledright: to theLVBM). NSWE (dashed line) (left: LSWE; We can model be seen Fig. 10 [50, 250] m andbathymetry thewave model arearight for x [−5, 50] −3 the linear model with an EBC implementation (solid line). we influx a periodic atprofile the (x =m. L)Comparwith Using the same in boundary the∈ previous case, we and for 20 the comparison 30 40 that50the wave 60 can see in is 70 slightly80dist [s] isons two at several time steps the profile (a) Paths influx a between periodic these wave at thesimulations right boundary (x = L) with the can persed of the first-order characteristic curves aresimulation shown by the thin in the LVBM. The dispersion ratio in the be seen in Fig. 10 (left: LSWE, right: LVBM). We can see area lines. 2 profile: 2 is given by µ = 0.0986 < 1, which is less dispersive η(L, = A sin (2πt/T , wave is slightly dispersed(64) in thet)comparison that) the in the than the previous case. In Fig. 11, the shoreline movement η(L, t) = A (2πt/T (64) LVBM. The ratioamplitude at the simulation area given caused by1 the periodic wave as well as the paths of charin which Asin =dispersion 0.05 m )is the and T = 20 s isis the robustly withforming consecutive interactions between the incoming by µ2 =A0.0986 < 1, characteristic which is lessfunction dispersive 0curves period. smoothened untilthan t = the 10 spreviis acteristic the shoreline are shown. Observinous which A In = 0.05 m the shoreline amplitudemovement and periodcaused T = 20by s. andresults the reflected Fig.this 11,is the usedcase. in influxing periodic wave. We use a uniform spatial the ing the of this wave. case, we can conclude that the EBC −1 Aperiodic smoothened function t = 10 s is used For simulation till the physical time t = 80 s, the comwavecharacteristic as well as the paths until of characteristic curves informing influxingthe thisshoreline periodic wave. We use uniform spatial grid −2 are shown. Observing the results of putational time for the coupled numerical solutions in both Nonlin. Geophys., 21,and 987–1005, 2014 m in the www.nonlin-processes-geophys.net/21/987/2014/ ∆x = 1 mProcesses simulation 0.015 this case,inwethecan concludearea that the ∆x EBC=technique can deal 850 domains is 1.83 times the physical time for the LSWE and −3 model area for the numerical solution of the NSWE. 20 30 40 50 60 70 80 t [s] (b) As in the previous case, we also choose the seaward boundary point at h0 = 10 m at the toe of the slope, that xs(t) [m]

840

Fig. 10. coupled (dashed EBC im (c) t = 250 another

250

−0.1 0

5.2 (d) Periodic wave 0.05

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

xs(t) [m]

η(x,t) [m]

η(x,t) [m]

For0.1simulation till the physical time t = 40 s, the com0.2 0.05 0.1 putational time for the coupled numerical solutions in both 0 0 −0.05 domains times the physical time for50 the100LSWE and250 1 0 is 503.3 100 150 200 250 0 150 200 x [m] x [m] (a) times the physical (b) 2.2 time for the LVBM, while the com0.1 0.1 putational time of the simulation using 0.05 an EBC only takes 0.05 0.5 825 0 0 0.12 times the physical time for the LSWE and 0.06 times −0.05 −0.05 0 50 100 150 200 250 0 50 100 150 200 250 that with x [m] Hence, we notice x [m] (c) for0 the LVBM. (d) that the simulation the EBC reduces the computational time significantly, down Fig. 9. Free-surface 97 profiles of solitarywith wavethe propagation are shown by approximately %, compared computational −0.5coupled LSWE with the NSWE (dashed and dotted-dashed for 10 20 30 entire domain. 40 50 comtimethe for0 the numerical models in the The t [s] lines), for the LSWE withLSWE an EBC implementation (solid line), and putational time for the with an EBC is slower than 830 for the LSWE with theand method of Choi et al. (2011) (solid line with the one with LVBM an EBC, because the internal time Fig. 8. Comparison of the shoreline movement of Choi et al. (2011) ’o’ at times (a) tstep = 10 (b) tin =MATLAB 20 s, (c)(solid t= 30 s, (d) stepmarker) ofline) theand ode45 time routine required (dashed the LSWE with ans,EBC simulation line) fora t = 40 s. The solid and dashed lines are on top of another. smaller time step 1t (compared to the LVBM) to preserve solitary wave case. the stability. The shoreline movement of our result compares well with 835 0.1 0.2 5.2 Periodic the 0.05 one of Choi wave et al. (2011). We can see the comparison in 0.1 0 Fig. 8. The solution of Choi et al. (2011) gives a higher pre0 −0.05 0 the 50 same 100 150 200 250 0 50 the 100previous 150 200 case, 250 Using bathymetry profile in we x [m] xthe [m] subsequent diction for the shoreline, but it cannot follow (b) (a) 0.1 0.1 influx a periodic wave at theberight boundary (x = L)thewith small positive wave. It may caused by neglecting re- the 0.05 0.05 profile: flection wave and nonlinear effects 0in their formulation. We 0 840 −0.05 −0.05 for several time steps in also compare the 150 free-surface profile 0 50 100 200 250 0 50 100 150 200 250 x [m] x [m] (c)η(L, t) = A sin (d) (2πt/T ) of the hard-wall boundary condi-(64) Fig. 9. The implementation tion at x = B in the method of Choi et propagation al. (2011) causes inFig. 9. Free-surface profiles of solitary wave are shown which Ain= 0.05 m isthethe amplitude period T =in20 s. accuracies the prediction of the point-wise height forin the coupled LSWE with NSWE (dashedand andwave dotted-dashed lines), for thedomain. LSWEcharacteristic with an EBC (solid andused A function t = line), 10waves s is thesmoothened entire In this case,implementation the effect until of reflected 845 forin the LSWE with method of Choi We et al.use (2011) (solid withgrid for the shoreline movement prediction is small, butline it may influxing thistheperiodic wave. uniform spatial ’o’∆x marker) at times t = 10 s, (b)area t = and 20ofs,waves (c) t= =arrives 30 s, (d) = the become a compound at = 1 important m in the(a)when simulation ∆x 0.015 mtthe in 40coastline. s. The solid and dashed lines are on top of another. model area for the numerical solution of the NSWE. η(x,t) [m]

835

0

η(x,t) [m]

(c)

−0.1

830

0.1

0.1 Kristina et al.: Effective Coastal Boundary Conditions foranother. Tsunami Wave Run-Up over Sloping Bathymetry of one 0 0

12

825

0.2

(d)

0 −0.05

0

−0.1

−0.1

0.1 η(x,t) [m]

0

0.1

(b)250

0.05 η(x,t) [m]

η(x,t) [m]

0

0.05

0.2

10050 150100 200150 250 200 x [m] x [m]

0.1

s

x (t) [m]

0

(a) (a)

0

η(x,t) [m]

0

η(x,t) [m]

0.05

η(x,t) [m]

0.1

η(x,t) [m]

1

η(x,t) [m]

0.05

0

−0.05

−0.15

η(x,t) [m]

0

Fig. 11. The shoreline movement of periodic waves for the linear

250

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry

999

2 1 885

xs(t) [m]

n are shown tted-dashed d line), and id line with 0 s, (d) t =

for the LVBM. Obviously, we notice that the simulation with the EBC reduces the computational time up to approximately 97 %, compared with the computational time for whole do- 880 main simulation.

−1 −2

s case, we L) with the

−3 20

30

40

(a)

50 t [s]

60

70

80

T = 20 s. 0 s is used patial grid 5 m in the E. e seaward slope, that is for x ∈ . Compare steps can We can see sed in the ea is given the prevised by the tic curves results of e can deal 850

890

Fig. 12. The three piece-wise linear bathymetry profile.

Figure 12. The three piece-wise linear bathymetry profiles.

2 1

(64)

0

5.3

−1 −2 −3 20

30

40

50 t [s]

60

70

80

The b Figur zero (95.0 the so an ea (95.8 Prese

η(x, 0

10 895

5

(b)

et al. mono al. (2 show wave that t est w

0

R(t) [m]

200

xs(t) [m]

0

coupled linear model (left: LSWE, right: LVBM) with the NSWE (dashed and dotted-dashed lines), and for the linear model with an 855 EBC implementation (solid line), at times (a) t = 20 s, (b) t = 40 s, (c) t = 60 s, (d) t = 75 s. The solid and dashed lines are on top of another at several plots.

and t the p w0 =

0 −5

Fig. 11. movementofofperiodic periodicwaves waves linear Figure 11.The The shoreline movement forfor thethe linear −10 model (a: LSWE, b: LVBM) coupled to the NSWE (dashed line) model (a: LSWE; b: LVBM) coupled to the NSWE (dashed line), 0 100 200 300 400 500 600 700 andforforthethelinear linear model with EBC implementation (solid line). t [min] and model with anan EBC implementation (solid line). Pathsofofthethefirst-order first-order characteristic curves shown Paths characteristic curves areare shown by by the the thinthin Figure The run-up heightofofperiodic periodicwaves waveswith with initial initial amplitude ampliFig. 13.13. The run-up height lines. lines. tude 1 m frequency and frequency = 0.0009 radThe s−1 .solid The solid A =A 1m=and ω =ω0.0009rad/s. line isline theisrunthe height run-up calculated height calculated by employing the LSWE model in the up by employing the LSWE model in the simusimulation area with the EBC implementation. The dashed one is lation area with EBC implementation. The dashed one is the result technique can deal robustly with consecutive interactions berobustly with consecutive interactions between the incoming thecoupling result of the coupling themodel NSWE model inarea thewith model area to the in of NSWE in model LSWE model tween the incoming and reflected waves. and the reflected wave. LSWE model in the simulation area. simulation area.

For t =t 80 s, the comForsimulation simulationtilltillthethephysical physicaltime time = 80 s, the computational time for the coupled numerical solutions in putational time for the coupled numerical solutions both in both domains thethe physical time forfor thethe LSWE, andand domainsis is1.83 1.83times times physical time LSWE 2.01 times the physical time for the LVBM, while the computational time of the simulation using an EBC only takes860 0.07 times the physical time of the LSWE, and 0.06 times that of the LVBM. Obviously, we notice that in the simulation with the EBC, the computational time decreases down by approximately 97 %, compared with the computational time for whole domain simulation. 865 As mentioned in the Introduction, resonant phenomena were investigated by Stefanakis et al. (2011) for monochromatic waves on a planar beach. Subsequently, Ezersky et al. (2013) used three piece-wise linear profiles of unperturbed depths (see Fig. 12), akin to a real coastal bottom topography, to find the analytical run-up amplification870 due to resonance effects. We follow this bathymetry profile, with tan α = 0.0036, tan β = 0.0414, h0 = 2500 m, and h1 = 200 m. These choices roughly characterize the Indian coast bathymetry (Neetu et al., 2011). The EBC point is located at the edge of the last beach slope. We influx a periodic wave (64) with amplitude A = 1 m and ω = 2π/T = 0.0009 rad s−1 . As a result, we get 10.67 times amplification, as shown in the run-up height R(t) in Fig. 13, while the result of Ezersky et al. (2013) gives about 12 times amplification. It should be noted that they use a linear approximation to www.nonlin-processes-geophys.net/21/987/2014/

As mentioned in the Introduction, resonant phenomena were investigated by Stefanakis et al.waves. (2011)Infor calculate the amplification of periodic ourmonochroresult, matic waves on a planar beach. Subsequently, Ezersky et al. the NSWE model is employed in the last beach slope region. (2013) usedofthree piece-wise linear profiles The period this wave is approximately 2 h, of andunperturbed it coindepths (seetheFig. 12), akin to aatreal coastal bottom topogracides with observed tsunami the Makran coast, accordphy, find et theal.analytical amplification to aresing totoNeetu (2011). In run-up nature, one would not due expect tsunami effects. of a monochromatic train. The investigation of tan onance We follow wave this bathymetry profile with Stefanakis et al. for the 25h0October 2010and Mentawai Isα = 0.0036, tan(2011) β = 0.0414, = 2500m, h1 = 200m. lands tsunami that thecharacterizing period of the dominant mode These choicesshowed are roughly the Indian coast of the incident wave is within the resonant regime, and it ex- at bathymetry (Neetu et al., 2011). The EBC point is located plained factlast thatbeach the highest is not driven wave by the(64) the edgethe of the slope.run-up We influx periodic leading and highest waves. with amplitude A = 1 m and ω = 2π/T = 0.0009 rad/s. As a result, we get 10.67 times amplification as shown in the run5.3 height Simulation bathymetry up R(t) inusing Fig. simplified 13, while Aceh the result of Ezersky et al. (2013) gives about 12 times amplification. It should be noted that they use a near linear approximation calculate the ampliThe bathymetry Aceh, Indonesia, isto displayed in Fig. 14. fication of periodic In our result, NSWE(genmodel Figure 14a concerns waves. bathymetry data from the GEBCO eral bathymetric chart of the oceans), with a zero value on land. Figure 14b concerns the cross-section at (95.0278◦ E, 3.2335◦ N)–(96.6583◦ E, 3.6959◦ N) shown by the solid line. The 2004 Indian Ocean tsunami was a result of an earthquake of magnitude of Mw = 9.1 at the epicenter point (95.854◦ E, 3.316◦ N), shown by the symbol in Fig. 14a. Presently, we Nonlin. Processes Geophys., 21, 987–1005, 2014

(a)

(b

Fig. 1 (95.0 conce proxim

14 Kristina et al.: Effective Coastal Boundary Conditions for T W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry

1000

η(x,t) [m]

0.6 0.3 0 −0.3 −0.6 0 0.5 0 −0.5 −1 −1.5

(a)

−0.5 20

40

60

80 100 x [km]

120

140

160

Fig. 15.15. The initial Figure The initialcondition conditionforforthe theAceh Acehcase caseisisshown shown for for the the linear model coupled to the NSWE (dashed and dotted-dashed lines) lines) an EBC over implementation (solid line). and for linear model with line). Kristina et al.: Effective Coastal Boundary Conditions forthe Tsunami Wave Run-Up Sloping Bathymetry The solid and dashed lines are on top of another. one another.

(b)

(c)

0

3 η(x,t) [m]

0

14

η(x,t) [m]

0

0.6 0.3 0 −0.3 −0.6 0

η(x,t) [m]

η(x,t) [m]

0.5

2 1 0

160

Fig.Figure 15. The condition the(a)Aceh is shown for(b)the 14. initial Bathymetry near for Aceh and case the cross-section at ◦ E, 3.6959 ◦ N). The solid linear model◦ coupled to◦the NSWE (dashed and dotted-dashed lines) (95.0278 E, 3.2335 N)–(96.6583 line andconcerns for the linear model with anand EBC (solidthe line). the bathymetry data, theimplementation dashed line concerns ap- 910 Theproximation solid and dashed lines are on top of another. used in the simulations.

905

910

915

920

925

930

consider the following initial N -wave profile The location of the EBC point is also determined from  d 2 Eq.η(x, (62). For (x)/S,  = 0.02with  1,f (x)= the linear is0 )valid 0)=Af expmodel −(x−x /w0 2for 915 dx h0  25.1 m. Hence, we choose the EBC point at depth h0 = and S = max (f (x)) , (65) 41.4 m, which is located at x = B = 12.4 km. Thus, the simulation area isthe forinitial x ∈ [12.4, 162.4] km, where and where velocity potential is zero.we Wefollow take Athe = real0.4 bathymetry of Aceh calculate propagation. m, the position of thetowave profilethe x0wave = 107.4 km, and a 920 It is coupled the model area for x ∈ [−8.6, 12.4] km, width w0 =with 35 km. whereThe a uniform = 1/300 is used to location slope of thewith EBCgradient point is γalso determined from calculate the For reflection andshoreline position. We use an irEq. (62).  = 0.02 1, the linear model is valid for p regular depththewith /h ash0the h0 grid 25.1according m. Hence, to wethe choose EBCratio point ath0depth = 41.4 m,ofwhich is located at x = traveling B = 12.4 km. simuladecrease the wavelength when fromThe a deep re- 925 tion area is thus x ∈ [12.4, region 162.4] km, follow gion with depth h tofor a shallower withwhere depth we h0 in linthe real bathymetry of Aceh to calculate the wave propagaear wave theory. The grid size used in the simulation area is It ismcoupled to the model areanear for xx ∈=[−8.6, 12.4] km, ∆xtion. = 305 at the shallowest point B. This choice where a uniform slope with gradient γ = 1/300 is used to of spatial resolution is fairly close to a tsunami numerical calculate the reflection and shoreline position. We use a nonsimulation (Horrillo et al., 2006 use ∆x = 100√m offshore 930 depth, with ratio simulations). h0 / h as the anduniform ∆x = grid 10 maccording onshore tointhe one dimensional decrease in the wavelength when traveling from a deep For the numerical solution of the NSWE in the model area,rea gion with depth h to a shallower region with depth h 0 in linuniform grid with ∆x = 3 m is used. ear wave theory. The grid size used in the simulation area is In Fig. 15, we show the initial profile. Comparisons be1x = 305 m at the shallowest point near x = B. This choice tween these two simulations at several time steps can be seen 935 of spatial resolution is fairly close to other numerical tsunami in Fig. 16. In this case, the wave elevation measured at B simulations (Horrillo et al., 2006 use 1x = 100 m offshore has been deformed from its initial condition due to the reand 1x = 10 m onshore in one-dimensional (1-D) simulaflection the numerical bathymetry beforeofentering the model area, tions).from For the solution the NSWE in the model seearea, Fig.a16a and b. We hardly see any differences between uniform grid with 1x = 3 m is used. the LSWE and LVBM simulations because the wavelength is much larger than the depth. The dispersion ratio at the iniProcesses 987–1005, 2014 tialNonlin. condition is givenGeophys., by µ2 =21, 0.002  1, and at the EBC 2 −5 point is approximately µ = 7.5 × 10  1. Therefore, the dispersion effect is not significant in this case. In Fig. 17, the shoreline position is displayed. From this plot, it is shown

η(x,t) [m] η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

η(x,t) [m]

0 −0.3

(f)

0

Fig. 16. Fr model (left by the dash model with times (a) t (f) 5400 s.

−1

−0

−1.2 −0.8 www.nonlin-processes-geophys.net/21/987/2014/ −0.4 0 0.4

xs(t) [km]

140

0.3

−0

0

0

1

(a)

−1

−0 xs(t) [km]

120

η(x,t) [m]

80 100 x [km]

η(x,t) [m]

60

0 0.6

η(x,t) [m]

40

0.6 0.3 0 −0.3

(e)

η(x,t) [m]

20

η(x,t) [m]

0

η(x,t) [m]

905

η(x,t) [m]

−0.5

The location of the EBC point is also determined from Eq. (62). For  = 0.02  1, the linear model is valid for 20 Hence, 40 60 80 20 40 60 h0  25.1 m. we choose the EBC point at depth h080 = x [km] x [km] (a) 0.6 which is located at x = B = 0.612.4 km. Thus, the sim41.4 m, 0.3 0.3 ulation0 area is for x ∈ [12.4, 162.4] km, where we follow the 0 −0.3 real −0.3 bathymetry of Aceh to calculate the wave propagation. −0.6 −0.6 0 20 40 60 80 0 20 40 60 80 x [km]the model area for x ∈ [−8.6, x [km] 12.4] km, (b) It is coupled with 0.5 where0.50 a uniform slope with gradient 0 γ = 1/300 is used to −0.5 −0.5 calculate the reflection and shoreline position.p We use an ir−1 −1 −1.5 grid according to the depth−1.5 regular with ratio h150 /h20as the 0 5 10 15 20 25 0 5 10 25 x [km] x [km] (c) decrease of the wavelength when traveling from a deep re3 3 2 2 gion with depth h to a shallower region with depth h0 in lin1 1 ear wave theory. The grid size used in the simulation area is 0 0 0 5 10 15 20 25 0 5 10 15 20 25 ∆x = 305 m at the shallowest point near x = B. This choice x [km] x [km] (d) 0.6 of spatial resolution is fairly close0.6to a tsunami numerical 0.3 0.3 0 simulation (Horrillo et al., 2006 use0 ∆x = 100 m offshore −0.3 −0.3 and ∆x = 1020m onshore in60one dimensional simulations). 0 40 0 20 40 60 x [km] x [km] (e) For the numerical solution of the NSWE in the model area, a 0.6 0.6 uniform grid with ∆x = 3 m is used. 0.3 0.3 0 0 In Fig. 15, we show the initial profile. Comparisons be−0.3 −0.3 0 20 simulations 40 60 several 0 time steps 20 tween these two at can40be seen60 x [km] x [km] (f) in Fig. 16. In this case, the wave elevation measured at B Fig. 16.16. Free-surface profiles Acehcondition simulations with the linear Figure Free-surfacefrom profiles of simulations with linear has been deformed its of initial due tothe the remodel (left: LSWE, right:LVBM) LVBM) coupled NSWE shown (left: LSWE; right: coupled totothethe NSWE areare shown flection from the bathymetry before entering the model area, by dashed and dotted-dashed lines and ofdifferences simulationsbetween for aa linlinear dashed lines, for seethe Fig. 16a and anddotted-dashed b. We hardly seeand anyfor ear model with an EBC implementation, are shown by the solid model with an EBC implementation are shown by the solid line the LSWE and LVBM simulations because the wavelength is at at times (a)t = t =than 800s,the s,(b) (b)depth. 1600 s,The (c) 2700 4000 times (a) 800 1600 2700s,s,(d) (d)3200 3200s,at s,(e)(e) 4000 much larger dispersion ratio the ini-s, s, and5400 (f) 5400 s. The solid and dashed lines are on topofofanother. one another. 2 lines (f) s. The solid and dashed are on top tial condition is given by µ = 0.002  1, and at the EBC point is approximately µ2 = 7.5 × 10−5  1. Therefore, the dispersion effect is not significant in this case. In Fig. 17, the In Fig. position 15, we show the initial profile. beshoreline is displayed. From this Comparisons plot, it is shown tween these two simulations at several time steps can be seen that the −1.2 wave runs up 1 km in the horizontal direction in apin Fig. 16. In10 thismin, case, the wave elevation measured at B proximately roughly in the time interval from 50has to −0.8 changed from its initial condition due to the reflection 60 min. For the given slope, it corresponds with a 3.3 m from run−0.4 the bathymetry before entering the model area; see Fig. 16a up height. 0 andFor b. We hardly anyphysical differences LSWE and simulation see till the timebetween t = 120the min, the comLVBM simulations because the wavelength is much greater putational0.4time for the coupled numerical solutions in both than the 0.8 depth. dispersion ratiotime at the condition domains is 0.03The times the physical forinitial the LSWE and 2 = 0.002  1, and at the EBC point, it is apis given by µ 0.03 times the physical time for the LVBM. While the com1.2 2 proximately = 20 7.5 10 1.using The 80 dispersion effect is 0µ of 40−5  60 100only 120 putational time the ×simulation an EBC takes t [min]In Fig. 17, the shoreline (a) therefore not significant in this case. η(x,t) [m]

0

0.6 0.3 0 −0.3 −0.6 0

xs(t) [km]

η(x,t) [m]

900

0.6 0.3 0 −0.3 −0.6 0

xs(t) [km]

0.5

900

0

(d)

−0

0

0

1

(b)

Fig. 17. Sh (a: LSWE, the linear m the first ord

he comin both WE and he comly takes

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry

1001

−1.2 −0.8 xs(t) [km]

ons bebe seen ed at B the redel area, between ength is the inihe EBC ore, the . 17, the s shown n in apm 50 to 3 m run-

model with an EBC implementation are shown by the solid line at times (a) t = 800 s, (b) 1600 s, (c) 2700 s, (d) 3200 s, (e) 4000 s, (f) 5400 s. The solid and dashed lines are on top of another.

−0.4 0 0.4 0.8 1.2 0

20

40

60 t [min]

80

100

120

20

40

60 t [min]

80

100

120

(a) −1.2 −0.8 xs(t) [km]

deep re0 in linn area is s choice merical offshore lations). l area, a

−0.4 0 0.4 0.8 1.2 0

(b)

Fig. 17. 17. Shoreline movement in in thethe Aceh case forfor thethe linear model Figure Shoreline movement Aceh case linear (a: LSWE, b: LVBM) coupled to thetoNSWE (dashed line)line), and for model (a: LSWE; b: LVBM) coupled the NSWE (dashed andlinear for themodel linear with model EBC implementation the anwith EBCanimplementation (solid(solid line).line). Paths of Paths of order the first-order characteristic are shown by lines. thin lines. the first characteristic curvescurves are shown by thin

position is displayed. From this plot, it is shown that the wave runs up 1 km in the horizontal direction in approximately 10 min, roughly in the time interval from 50 to 60 min. For the given slope, it corresponds to a 3.3 m run-up height. For simulation till the physical time t = 120 min, the computational time for the coupled numerical solutions in both domains is 0.03 times the physical time for the LSWE and 0.03 times the physical time for the LVBM, while the computational time of the simulation using an EBC only takes 0.003 times the physical time for the LSWE and 0.004 times that for the LVBM. We again notice that the simulations using the EBC decrease the computational times down by approximately 92 % of the computational times with the coupled model in the entire domain. In this case, the simulation with the LSWE is faster, as expected, since the LVBM involves more calculations within the same time step. For the case when breaking occurs, we use the same profile, with an amplitude twice as high (A = 0.8 m). In Fig. 18, the shoreline position is presented. Compared to the numerical NSWE solution, it can be seen that the shoreline movement is well represented by the characteristic curves, while the shoreline position xs (t) given by Eq. (48) gives a less accurate result. Breaking occurs when two incoming characteristic curves intersect before reaching the shoreline. As can be seen in the right figure, the first breaking takes place at approximately t = 45 min. The corresponding free-surface profiles for several times before and after breaking are shown in Fig. 19. www.nonlin-processes-geophys.net/21/987/2014/

Figure 18. Shoreline movement (a) and an inset (b) of a breaking wave simulation. The linear model coupled to NSWE is shown by the dashed line, while the solid one is the shoreline movement of a linear model simulation with EBC implementation. Paths of the first-order characteristic curves are shown by thin lines. Breaking occurs when two incoming characteristic curves intersect before reaching the shoreline. It is indicated by the red oval at approximately t = 45 min.

6

Conclusions

We have formulated a so-called effective boundary condition (EBC), which is used as an internal boundary condition within a domain divided into simulation and model areas. The simulation area from the deep ocean up to a certain depth at a seaward boundary point at x = B is solved numerically using the linear shallow water equations (LSWE) or the linear variational Boussinesq model (LVBM). The nonlinear shallow water equations (NSWE) are solved analytically in the model area from this boundary point towards the coastline over linearly sloping bathymetry. The wave elevation at the seaward boundary point is decomposed into the incoming signal and the reflected one, as described in Antuono and Brocchini (2007, 2010). The advantages of using this EBC are the ability to measure the incoming wave signal at the boundary point x = B for various shapes of incoming waves, and thereafter to calculate the wave run-up and reflection from these measured data. To solve the tsunami wave run-up in the nearshore area analytically, we employ the asymptotic technique for solving the NSWE over sloping bathymetry derived by Antuono and Brocchini (2010), applied to any given wave signal at x = B.

Nonlin. Processes Geophys., 21, 987–1005, 2014

1002

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry

Figure 19. Free-surface profiles of a breaking wave simulation for the linear model coupled to the NSWE (dashed and dotted-dashed lines) and for the linear model with an EBC implementation (solid line) at t = 40–70 min. The solid and dashed lines are on top of one another.

An extension of this EBC technique to the case where the NSWE model is used both in the simulation and the model areas follows directly from the variational methodology. The analytical benchmark for this case is provided by Carrier et al. (2003) and Kâno˘glu (2004). The two-dimensional (2-D) extension of this technique can be formulated asymptotically using an approach by Ryrie (1983). For waves incident at a small angle to the beach normal, the onshore problem can be calculated using the analytical 1-D run-up theory of the nonlinear model, and independently, the longshore velocity can be computed asymptotically. By using a 2-D linear model in the open sea towards the seaward boundary line (i.e., in the simulation area), and by employing this approach in the model area, we can in principle apply the EBC method for this 2-D case as well. This will be approximately valid for 2-D flow with slow variations along the EBC line. The EBC formulation for the case when the shoreline is fronted by a vertical wall as presented by Kâno˘glu and Synolakis (1998) can be obtained by requiring the normal velocity at the shoreline wall boundary to be zero. Another characteristic of the outgoing or reflected waves must be derived (either for the LSWE model or the NSWE model), but the coupling between the numerical and analytical models remains the same as that derived in this article.

The EBC implementation has been verified in several test cases by comparing simulations in the whole domain (using numerical solutions of the LSWE/LVBM in the simulation area coupled to the NSWE in the model area) with ones using the EBC. We have also validated our approach with the laboratory experiment of Synolakis (1987) for the run-up of a solitary wave over a plane beach. The location of the boundary point x = B is considered before the nonlinearity plays an important role in the wave propagation. The comparisons between both simulations show that the EBC method gives a good prediction of the wave run-up as well as the wave reflection, based only on the information of the wave signal at this seaward boundary point. It is also shown that the EBC technique can capture the resonance effect that occurs due to the incoming and reflected wave interactions. The computational times needed in simulations using the EBC implementation show a large reduction compared to times required for the corresponding full numerical simulations. Hence, without losing the accuracy of the results, we could compress the time needed to simulate wave dynamics in the nearshore area.

Nonlin. Processes Geophys., 21, 987–1005, 2014

www.nonlin-processes-geophys.net/21/987/2014/

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry Appendix A: Finite volume implementation The conservative form of NSWE is given by ∂u ∂f(u) + = s, ∂t ∂x with     hu hu2 + 12 gh2 u= , f(u) = h hu

(A1)

(A2)

(A3)

The system (A1) is discretized using a Godunov finite volume scheme. First, the domain [A, B] with some fixed A < xs (t) is partitioned into N grid cells, with grid cell k occupying xk− 1 < x < xk+ 1 . The Godunov finite volume 2 2 scheme is derived by defining a space–time mesh with elements xk− 1 < x < xk+ 1 and tn < t < tn+1 , and integrating 2 2 Eqs. (A1) over this space–time element: xk+ 1

Z

Z

2

Ztn+1 f(u(xk− 1 , t))dt− f(u(xk+ 1 , t))dt 2

tn

2

− 12 gh2

(k− 12 )+

0

tn xk− 1

! ,

2

 = max hk+1 + bk+1 − bk+ 1 , 0 , 

h

k+ 12

s dxdt.

(A4)

+

2

(A11)

2

2

In the grid cells, we define the mean cell average Uk = Uk (t) as xk+ 1 2

u(x, t)dx,

(A5)

xk− 1 2

with cell length hk = xk+ 1 −xk− 1 . The function Uk is piece2 2 wise constant in each cell. A numerical flux F is defined to approximate the flux f: 1 ≈ 1t

and bk+ 1 = max(bk , bk+1 ).

tn xk− 1

F

s dxdt ≈ Sk = 1t

1 2 2 gh(k+ 1 )−

2

+




1

k+ 2

xk+ 1

Unk , Unk+1

x

Ztn+1 Zk+ 2

2

tn

Z

The topographic term s is discretized as

with the water depths h(k+ 1 )− and h(k+ 1 )+ chosen as fol2 2 lows, to ensure the non-negativity of these depths:   h 1 − = max hk + bk − bk+ 1 , 0 ,

2

1 Uk (t) := hk

(A9)

k+ 2 +

(A10) u(x, tn )dx=

Ztn+1

Ztn+1 Z

where the interface values are given by     u h k 1 k+ − 2  and Un 1  =  h 1  k+ 2 − k+ 2 −     u h k+1 1 + k+ 2 . Un 1  =  h 1  k+ 2 +

2

xk− 1

2

k+ 2

k+ 2

2

u(x, tn+1 )dx− xk− 1

which is a forward Euler explicit method. To ensure that the depth is non-negative and that the steady state of a fluid at rest is preserved, the approach of Audusse et al. (2004) is used. The numerical flux F is then defined as !
n n n n (A8) F Uk , Uk+1 = Fk+ 1 U 1 − , U 1 + , 2

and the topographic term   −gh db/dx s= . 0

xk+ 1

1003

Ztn+1 f(u(xk+ 1 , t))dt.

(A6)

2

tn

By using Eqs. (A5)–(A6), expression (A4) then becomes


1t Un+1 =Unk − F Unk , Unk+1 − F Unk−1 , Unk k hk x

1 + hk

1

Ztn+1 Zk+ 2 s dxdt, tn xk− 1 2

www.nonlin-processes-geophys.net/21/987/2014/

(A7)

(A12)

The discretization of the shallow water equations thus reads
1t  Un+1 = Unk − F 1 Un 1 − , Un 1 + k hk k+ 2 k+ 2 k+ 2
 1t − Fk− 1 Un 1 − , Un 1 + + Sk . (A13) 2 hk k− 2 k− 2 The Harten–Lax–van Leer (HLL) flux (Harten et al., 1983; Toro et al., 1994) is used as the numerical flux. It is given by   if 0 < SL FL SR FL −SL FR +SL SR (UR −UL ) HLL Fk+ 1 = if SL ≤ 0 ≤ SR . SR −SL  2  FR if 0 > SR (A14) The wave speeds SL and SR are approximated as the smallest and largest eigenvalues at the corresponding node. To ensure the stability of this explicit scheme, a Courant–Friedrichs– Lewy (CFL) stability condition per cell is used for all eigenvalues λp at each Unk : 1t
λp Un ≤ 1. (A15) k h k Nonlin. Processes Geophys., 21, 987–1005, 2014

1004

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry

Appendix B: Coupled model The (continuous Galerkin) finite element implementation of LSWE or LVBM uses linear polynomials for solving φ, ψ, and η approximately, while the finite volume implementation for NSWE approximates h and u with a constant value. Since u = ∂x φ, the velocity of the two models is approximated with the same order of polynomials. By coupling both models, in the simulation area, we can rewrite Eq. (52) as Mkl η˙ l − Skl φl − Bkl ψl − δk1 (hu)|x=B − = 0 Mkl φ˙ k + gMkl ηk = 0 ! β˜ Akl ψl + Bkl φl + Gkl ψl − δk1 hu |x=B − = 0. hb

(B1a) (B1b) (B1c)

In the finite volume implementation, the boundary condition is inserted through the numerical flux at x = B by using the coupling condition (16) as follows:     ˘ x ψ˘ hu hb ∂x φ˘ + β∂ = . (B2) h hb + η˘

Acknowledgements. We are pleased to acknowledge funding for Wenny Kristina from the Netherlands Organization of Scientific Research (NWO), Earth Sciences division. The idea to couple potential flow (dispersive) models in a simulation area with the full NSWE in the onshore area came originally from the late Howell Peregrine, and was used by Vijaya Ambati, Onno Bokhove, and Frank Klaver (Klaver, 2009). Edited by: V. Shrira Reviewed by: three anonymous referees

References Antuono, M. and Brocchini, M.: The boundary value problem for the nonlinear shallow water equations, Stud. Appl. Math., 119, 73–93, 2007. Antuono, M. and Brocchini, M.: Solving the nonlinear shallowwater equations in physical space, J. Fluid Mech., 643, 207–232, 2010. Audusse, E., Bouchut, F., Bristeau, M. O., Klein, R., and Perthame, B.: A fast and stable well-balaced scheme with hydrostatic reconstruction for shallow water flows, SIAM J. Sci. Comput., 25, 2050–2065, 2004. Bokhove, O.: Flooding and Drying in Discontinuous Galerkin Finite-Element Discretizations of Shallow-Water Equations. Part 1: One Dimension, J. Sci. Comput., 22, 47–82, 2005. Brocchini, M. and Peregrine, D. H.: Integral flow properties of the swash zone and averaging, J. Fluid Mech., 317, 241–273, 1996. Carrier, G. F. and Greenspan, H. P.: Water waves of finite amplitude on a sloping beach, J. Fluid Mech., 4, 97–109, 1958. Carrier, G. F., Wu, T. T., and Yeh, H.: Tsunami runup and drawdown on a plane beach, J. Fluid Mech., 475, 79–99, 2003.

Nonlin. Processes Geophys., 21, 987–1005, 2014

Choi, B. H., Kaistrenko, V., Kim, K. O., Min, B. I., and Pelinovsky, E.: Rapid forecasting of tsunami runup heights from 2-D numerical simulations, Nat. Hazards Earth Syst. Sci., 11, 707– 714, doi:10.5194/nhess-11-707-2011, 2011. Choi, B. H., Pelinovsky, E., Kim, K. O., and Min, B. I.: Estimation of run-up Heights of the 2011 off the Pacific Coast of Tohoku Earthquake Tsunami Based on Numerical Simulations, The Open Oceanography Journal, 6, 5–13, 2012. Cotter, C. J. and Bokhove, O.: Variational water-wave model with accurate dispersion and vertical vorticity, J. Eng. Math., 67, 33– 54, 2010. Didenkulova, I. and Pelinovsky, E.: Run-up of long waves on a beach: the influence of the incident wave form, Oceanology, 48, 1–6, 2008. Ezersky, A., Tiguercha, D., and Pelinovsky, E.: Resonance phenomena at the long wave run-up on the coast, Nat. Hazards Earth Syst. Sci., 13, 2745–2752, doi:10.5194/nhess-13-2745-2013, 2013. Gagarina, E., van der Vegt, J., and Bokhove, O.: Horizontal circulation and jumps in Hamiltonian wave models, Nonlin. Processes Geophys., 20, 483–500, doi:10.5194/npg-20-483-2013, 2013. Harten, A., Lax, P. D., and Van Leer, B.: On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Rev., 25, 35–61, 1983. Horrillo, J., Kowalik, Z., and Shigihara, Y.: Wave Dispersion Study in the Indian Ocean-Tsunami of December 26, 2004, Mar. Geod., 29, 149–166, 2006. Kâno˘glu, U.: Nonlinear evolution and runup and rundown of long waves over a sloping beach, J. Fluid Mech., 513, 363–372, 2004. Kâno˘glu, U. and Synolakis, C. E.: Long wave runup on a piecewise linear topographies, J. Fluid Mech., 374, 1–28, 1998. Klaver, F.: Coupling of numerical models for deep and shallow water, M.Sc. thesis, University of Twente, Enschede, the Netherlands, 2009. Klopman, G., Van Groesen, E., and Dingemans, M.: A variational approach to Boussinesq modelling of fully non-linear water waves, J. Fluid Mech., 657, 36–63, 2010. Kristina, W., Van Groesen, E., and Bokhove, O.: Effective Coastal Boundary Conditions for Dispersive Tsunami Propagation, Memorandum 1983, Department of Applied Mathematics, University of Twente, Enschede, the Netherlands, 2012. Li, Y. and Raichlen, F.: Solitary wave run-up on plane slopes, J. Waterw. Port C. Div., 127, 33–44, 2001. Liu, Y., Shi, Y., Yuen, D. A., Sevre, E. O. D., Yuan, X., and Xing, H. L.: Comparison of linear and nonlinear shallow wave water equations applied to tsunami waves over the China sea, Acta Geotechnica, 4, 129–137, 2009. Luke, J. C.: A variational principle for fluids with a free surface, J. Fluid Mech., 27, 395–397, 1967. Madsen, P. A. and Schaffer, H. A.: Analytical solutions for tsunami run-up on a plane beach: single waves, N-waves and transient waves, J. Fluid Mech., 645, 27–57, 2010. Madsen, P. A., Murray, R., and Sorensen, O. R.: A new form of the Boussinesq equations with improved linear dispersion characteristics, Coast. Eng., 15, 371–388, 1991. Mei, C. C.: The applied dynamics of ocean surface waves, World Scientific Publishing Singapore, 1989. Miles, J. W.: On Hamilton’s principle for surface waves, J. Fluid Mech., 83, 153–158, 1977.

www.nonlin-processes-geophys.net/21/987/2014/

W. Kristina et al.: EBC for tsunami wave run-up over sloping bathymetry Neetu, S., Suresh, I., Shankar, R., Nagarajan, B., Sharma, R., Shenoi, S. S. C., Unnikrishnan, A. S., and Sundar, D.: Trapped waves of the 27 November 1945 Makran tsunami: observations and numerical modeling, Nat. Hazards, 59, 1609–1618, doi:10.1007/s11069-011-9854-0, 2011. Pelinovsky, E. N. and Mazova, R. Kh.: Exact Analytical Solutions of Nonlinear Problems of Tsunami Wave Run-up on Slopes with Different Profiles, Nat. Hazards, 6, 227–249, 1992. Ryrie, S. C.: Longshore motion generated on beaches by obliquely incident bores, J. Fluid Mech., 129, 193–212, 1983. Stefanakis, T. S., Dias, F., and Dutykh, D.: Local run-up amplification by resonant wave interactions, Phys. Rev. Lett., 107, 124502, doi:10.1103/PhysRevLett.107.124502, 2011. Synolakis, C. E.: The run-up of solitary waves, J. Fluid Mech., 185, 523–545, 1987.

www.nonlin-processes-geophys.net/21/987/2014/

1005

Titov, V. V., Moore, C. W., Greenslade, D. J. M., Pattiaratchi, C., Badal, R., Synolakis, C. E., and Kâno˘glu, U.: A New Tool for Inundation Modeling: Community Modeling Interface for Tsunamis (ComMIT), Pure Appl. Geophys., 168, 2121–2131, 2011. Toro, E. F., Spruce, M., and Speares, W.: Restoration of the contact surface in the HLL-Riemann solver, Shock Waves, 4, 25–34, 1994. van Groesen, E.: Variational Boussinesq Model, part 1: Basic equations in Cartesian coordinates, Technical Report of LabMathIndonesia, 2006. Wei, Y., Bernard, E. N., Tang, L., Weiss, R., Titov, V. V., Moore, C., Spillane, M., Hopkins, M., and Kâno˘glu, U.: Realtime experimental forecast of the Peruvian tsunami of August 2007 for U.S. coastlines, Geophys. Res. Lett., 35, L04609, doi:10.1029/2007GL032250, 2008. Zakharov, V. E.: Stability of periodic waves of finite amplitude on the surface of a deep fluid, J. Appl. Mech. Tech. Phy., 9, 190–194, 1968.

Nonlin. Processes Geophys., 21, 987–1005, 2014