Effective Coastal Boundary Conditions for Tsunami Wave Run-Up over Sloping Bathymetry

. An effective boundary condition (EBC) is introduced as a novel technique to predict tsunami wave run-up along the coast and offshore wave reﬂections. Numerical modeling of tsunami propagation at the coastal zone has been a daunting task since high accuracy is needed to capture as-5 pects of wave propagation in the more shallow areas. For example, there are complicated interactions between incoming and reﬂected waves due to the bathymetry and intrinsically nonlinear phenomena of wave propagation. If a ﬁxed wall boundary condition is used at a certain shallow depth con-10 tour, the reﬂection 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 ei-15 ther 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 20 run-up heights at the shore and the reﬂection caused by the slope. The reﬂected wave is then inﬂuxed 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) conser-25


Introduction
Shallow water equations are widely used in the modeling of tsunamis since their wavelengths (typically 200km) are far greater than the depth of the ocean (typically 2 to 3km).Tsunamis also tend to have a small amplitude offshore, which 35 is why they generally are less noticeable at sea.Therefore, linear shallow water equations (LSWE) often suffice as a simple model of tsunami propagation (Choi et al., 2011;Liu et al., 2009;Kânoglu and Synolakis, 1998).On the contrary, it turns out that the lack of dispersion is a shortcoming of 40 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 short amount of 45 computation.However, as the tsunami approaches the shore, shoaling effects cause a decrease of the wavelength and an increase of 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 50 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 resolution to preserve the accuracy of the solution near the coast area (Titov et al., 2011;Wei et 55 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 this shortcoming by an investigation of a so-called effective 60 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 ∈ [x s (t), B] with shoreline position x s (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 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 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 to estimate tsunami run-up heights is also suggested by Choi et al. (2011Choi et al. ( , 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 runup height of tsunami waves at the coast is subsequently calculated in separation by employing a linear approach.It is claimed that the linear and nonlinear theories predict the 105 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) shows that there is a difference in the maximum run-up prediction between linear and nonlinear theory.In addition to calculating only the maxi-110 mum run-up height as in Choi's method, our EBC also includes the calculation of reflected waves.Thus, the pointwise wave height in the whole domain (offshore and onshore area) is predicted accurately.For the inundation prediction, we have verified that the method introduced by Choi et al.

115
(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 120 at later times.Stefanakis et al. (2011) investigate 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 125 why in some cases it is not the first wave that results in the highest run-up.
A 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 130 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 sea-135 ward 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 140 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 of the wavelength in a linear approach, and thereafter estimate the location of the EBC point.

145
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 asymp-150 totic closed-form solution, we also numerically simulate the NSWE in the model area using a finite volume method (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 in-troduce 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 coupling conditions derived 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.

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 offshore coordinate x.The artificial boundary is then placed at x = B, while the real (time-dependent) boundary lies at x = x s (t) with x s (t) < B. For example, at rest land starts at x = 0 where the total water depth h(0, t) = 0.This water line is time dependent as the wave can move up and down the beach.
We will restrict 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=δ (1) The approximation for the velocity potential Φ in Eq. ( 1) can be of various kind, but all are based on the idea to restrict the class of wave motions to a class that contains the wave motions one is interested in (van Groesen, 2006;Cotter and Bokhove, 2010;Gagarina et al., 2013).Following Klopman et al. (2010), we approximate the velocity potential as follows 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 con-205 dition at the bottom boundary z +h b (x) = 0.The latter kinematic condition yields For slowly varying bottom topography, this condition is approximated by 210 Substitution of Eq. ( 2) into Eq.( 1) yields the variational principle for Boussinesq equations, as follows (Klopman et al., 2010) where functions β(x), α(x), and γ(x) are given by (5)

220
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 dis-

225
persive 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 limiting system.We a priori divide the domain into two intervals, We choose a parabolic profile function b , in which the x dependence is considered to be parametric when the total water depth h is sufficiently slowly varying.The coefficients in (5) simplify to their linearized counterparts in the simulation area where the linear Boussinesq equation holds (while these coefficients disappear in the model area where the nonlinear depth-averaged shallow water equations hold).
Since the variations are arbitrary, the linear equations emerging from Eq. (8b) for x ∈ [B, L] are as follows and for x ∈ [x s (t), B], we get the nonlinear equations of motion The last two terms in Eq. (8b) are the boundary terms at x = x s .They can be rewritten as follows Eq. ( 10b) into (11), the boundary condition at the shoreline is i.e., the velocity of the shoreline equals the horizontal ve-295 locity 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 using the variations δ φ and δ ψ.Because the depth-averaged 300 shallow water equations are considered, we have where the last equality arises from approximation (2) for the velocity potential.Thus, the variation of δφ becomes Substituting this into Eq.(8b), we get the coupling condition at x = B for the linear model as follows 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 φ = φ and thus δφ = δ φ.From Eq. (8b), the coupling condition for nonlinear model is given by Note that the coupling conditions ( 15)-( 16) are used to transfer the information between the two domains.The coupling conditions (15) gives the information of φ and ψ in the simulation area, provided the information of φ from the model area.Meanwhile, the coupling condition ( 16) gives the information of φ in the model area, provided the information of φ and ψ from the simulation area.
3 Nonlinear Shallow Water Equations

Characteristic form
We will start with the NSWE in the shore region.Using η = −h b + h and velocity u = ∂ x φ, we may rewrite Eq. ( 10) as follows (starred variables are used here for later convenience) 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): in which h 0 is the still water depth at the seaward boundary and u 0 , l 0 , and t 0 are defined below as where g = 1 and γ = 1 are dimensionless gravity acceleration and beach slope, respectively.The NSWE in dimensionless form are then given by 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, 1957;Synolakis, 1987;Pelinovsky and Mazova, 1992;Carrier et al., 2003;Kânoglu , 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 in which c = √ gh and Variables α and β are the so-called Riemann invariants since they do not change their value along the characteristics curves in Eq. ( 21).Assuming the flow to be subcritical 370 (that is |u| < c), the first characteristics 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 at the shoreline).375

A trivial solution of characteristic curve
In the trivial case of no motion (u = η ≡ 0) as well as the dynamic case presented later, we focus on the incoming characteristic curve.In the rest case, it is given by 380 For x = 0, substituting y = √ gγx results in the general solution for variable y as follows 385 with a constant C 2 .When the curve intersects x = B at time t = τ , with h 0 the depth at x = B such that h 0 = γB and y(B) = √ gγB = c 0 , the particular solution is given by 390 In case of no motion, the boundary data α = α 0 (τ ) and β = β 0 (τ ) are as follows Transforming back to the x variable while using these ex-395 pressions, we get the incoming characteristic curve with ω = c 0 /(gγ).Along this characteristic curve, the Riemann invariant is constant.As in our previous paper (Kristina et al., 2012), the characteristic curve of the LSWE are given by dx/dt = ±c 0 .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 for which the characteristic reaches the shoreline position, here x = 0, at time t.It is given by the condition As displayed in Fig. 2, the incoming characteristic curves that reach the shoreline at time t, intersect x = B = 1 at time τ = t − 2 (ω = 1 in this case).Since u = 0 in the rest case, the boundary condition ( 12) is of course satisfied.

Boundary Value Problem (BVP)
Li and Raichlen (2001) and Synolakis (1987) combine linear and nonlinear theory to reduce the difficulties in the assignment 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 follows 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 (1957).We do, however, clarify the 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 waves traveling back into the deeper waters.If the sea is assumed in the rest state during the initial condition, the data of η(B, t)=u(B, t)=0 for t < 0. In accordance to the previous trivial case, the initial time where a characteristic meets x = B is labeled as τ and we write x = χ(t, τ ), so we have the data α = α 0 ≡ 2c(B, τ ) − u(B, τ ) + gγτ along the incoming 440 characteristic curves and β = β 0 ≡ 2c(B, τ )+u(B, τ )−gγτ along the outgoing characteristic curves.Then we can rewrite Eq. ( 21) as 445 which means that the boundary values are carried by the incoming and outgoing characteristic curves.To be concise, we write χ t = ∂ t χ and χ τ = ∂ τ χ.Our aim is to obtain a closed equation for the dynamics and we focus on the incom-450 ing characteristic by fixing α = α 0 .We can rewrite Eq. (29a) as follows Here β = β(χ, t) since we are moving along an incoming 455 characteristic curve.By taking the total t derivative of β, we obtain in which the last equality comes from Eq. ( 30).In addition, 460 the τ -derivative of Eq. ( 30) gives We still need an explicit expression for β t which can be obtained by rewriting Eq. (21b) in the following way 465 Combining Eqs. ( 31)-( 33), we get the following differential equation for the incoming characteristic curves: with boundary conditions 475 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 con-480 sistency χ τ must be zero at the shoreline since generally χ tt = gγ.

Boundary value assignment 530
In order to calculate the unknown functions A 1 (s) and A 2 (s), we need to assign the boundary conditions (34).In (ν, ξ) space, t = τ corresponds to ν = −2ω, and by imposing the first boundary condition, we get The second boundary condition is given by Eq. (36c) in which is evaluated at τ = τ 0 , i.e., ν = 0 needs to be finite.Evaluating Eq. ( 44) at ν = 0 gives us convergence when the co-545 efficient A 2 (s) is zero, which avoids an unbounded result.Hence, from the first boundary condition (36b), coefficient A 1 (s) is given by 550 Thus, the solution of the incoming characteristic curves at first order is given by 555 The shoreline position must satisfy χ τ | τ =τs = 0, and in the first order approximation it is given by Since τ = τ 0 corresponds to ν = 0 and ξ = t − 2ω, we get

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 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 with ϕ 1 the basis function in element 0 for x ∈ [x s , B] and with φ(B, t) = η(B, t) = 0.For linear polynomials, use of Eq. ( 49) into Eq.( 6) yields where we introduced mass and stiffness matrices M kl , S kl , A kl , B kl , G kl , 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 Provided we let the size of the zeroth element go to zero such that the underline terms in Eq. (50b) vanish, the equations arising from Eq. ( 50) are M kl φk + gM kl η k = 0 (52b) with Kronecker delta symbol δ kl (one when k = l and zero 610 otherwise) and Eq. ( 10) for x ∈ [x s , 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 that uses its internal time step.

615
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 625 Note that for = 0, we indeed find the rest depth h b (x) = γx.The function χ (1) t follows from evaluation of Eq. ( 46) and since t = τ is equivalent to ν = −2ω, we immediately obtain Thus, the solutions of h and u at t = τ are given as follows In order to calculate the solution 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 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).Thus, by imposing the effective boundary condition (EBC) and choosing the location x = B before the nonlinearity arises, we do actually solve for the perturbation expansion in the nonlinear area, but we do not perturb the incoming wave data.
The values ηh and ȗ 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 determined by using the approximation (13) as follows which is the limit from the right at node 1.
The solutions of η = h − h b and u in Eq. ( 55) account for the reflected wave, so we may define for η I and η R are the wave elevations of the incoming and reflected wave respectively at x = B.This superposition is also described in Antuono and Brocchini (2007;2010) and actually in line with our EBC concept since the 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 which is calculated using approximation (57).Thus, we can 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 let it travel over the real bathymetry to the seaward boundary point x = B. From Eq. ( 55), the expres-680 sions for the reflected wave are as follows where the Fourier transform and its inverse for any incoming 685 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 690 In order to verify the EBC implementation, we perform numerical simulations with a code that couples the LSWE in the simulation area with the NSWE in the model area (Bokhove,695 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.
5 Study Case 700 Three test cases are considered.The first one is a synthetic one concerning a solitary wave, such that we can compare with other results.Subsequently, we consider periodic wave influx as the second case to test the robustness of the technique when there is continuous interaction between the in-705 coming and reflected wave.The third case is a more realistic one concerning tsunami propagation and run-up based on simplified bathymetry at the Aceh coastline.The location of the EBC point is determined from the linearity condition = A 0 /h 0 1.From linear theory, the 710 wave amplification over depth is given by the ratio A 0 = A 4 h/h 0 , where A and h are the initial wave amplitude and depth.Hence, the EBC point must be located at depth 715 Since a dispersive model is also used in the simulation area, we will discuss the dispersion effect at this EBC point as well.The non-dispersive condition is given by µ 2 = (k 0 h 0 ) 2 1, where k 0 = 2π/λ 0 is the wavenumber and λ is the wavelength.In linear wave theory, the wavelength decreases with the square root of the depth when running in shallower water, that is λ 0 = λ h 0 /h.Thus, using this relation we can investigate the significance of the dispersion given the information of the initial condition and bathymetry profile.

Solitary wave
The run-up of a solitary wave is studied by means of the wellknown case of Synolakis (1987).A solitary wave centered at x = x 0 at t = 0 has the following surface profile: We benchmark the EBC implementation and the coupling of numerical solutions with experimental data of Synolakis (1987) provided at NOAA Center for Tsunami Research (http://nctr.pmel.noaa.gov/benchmark/).Solitary wave runup over a canonical bathymetry is considered with the scaled amplitude A = 0.0185 and κ = 3A/4 = 0.1178.The initial condition is centered at x 0 = 37.35 over the bathymetry with a constant slope γ = 1/19.85for x < 19.85.The EBC point is located at x = 10 such that the domain is divided into the model area for x ∈ [−5, 10] and the simulation area for x ∈ [10, 80].The spatial grid size is ∆x = 0.25 in the simulation area and ∆x = 0.0125 in the model area for the numerical solution of NSWE.In all cases, several spatial resolutions have been applied to verify numerical convergence.For the time integration, we use the fourth order ode45 solver that uses its own time step in MATLAB.
The simulations with the EBC implementation and the coupling of numerical solutions are only presented for the LSWE model in the simulation area since the initial condition has long wavelength and thus dispersion effects will not appear.Figure 3 shows the time evolution of this profile for scaled time t = 30 − 70 with 10 increments.It can be seen that the EBC implementation and the coupling of numerical solutions agree well with the laboratory data.The comparison of the shoreline movement between the simulation with EBC implementation and the coupling of numerical solutions is shown in Fig. 4. For the simulation till the scaled physical time t = 100, the computational time for the coupled numerical solutions in both domains is 10.9 s.While the computational time of the simulation with the EBC implementation only takes about 18 % of that time.Hence, we notice that the simulation with the EBC reduces the computational time significantly, up to approximately 82 %, compared with the computational time in the whole domain.
In order to show the dispersion effect, we consider a shorter wave with the profile given in Eq. ( 63) for κ = 0.04, x 0 = 150 m, and A = 0.1 m.The bathymetry is given by constant depth 10m for x > 50 m, continued by a constant slope γ = 1/5 towards the shore.A uniform spatial grid ∆x = 1m is used in the simulation area and ∆x = 0.015 m in the model area for the numerical solution of the NSWE.Evaluating Eq. ( 62) for = 0.02 1, the EBC point must be located at h 0 3.3 m.Accordingly, we choose this seaward boundary point at h 0 = 10 m at the toe of the slope, that is     In Fig. 5, we can see the initial profile of the solitary wave.Comparisons between these two simulations at several time steps can be seen in Fig. 6 (left: LSWE, right: LVBM).Comparing the left and right figures, we can see that the wave is slightly dispersed in the LVBM.Because we have flat bathymethy in this case, the dispersion ratio at the simulation area is constant and given by µ 2 = 0.39 < 1.Hence, it is shown that the long waves propagate faster than the shorter ones in LVBM simulations.In Fig. 7, the shoreline movement caused by this solitary wave is shown.The paths of characteristic curves forming the shoreline are also presented in this figure.We can see that the shoreline is formed by the envelope of the characteristic curves.The result with the LVBM shows a lower run-up but higher run-down with some oscillations at later times.
For simulation till the physical time t = 40 s, the computational time for the coupled numerical solutions in both domains is 3.3 times the physical time for the LSWE and  2.2 times the physical time for the LVBM.While the computational time of the simulation using an EBC only takes 0.12 times the physical time for the LSWE and 0.06 times for the LVBM.Hence, we notice that the simulation with 800 the EBC reduces the computational time significantly, up to approximately 97 %, compared with the computational time for the numerical models in the entire domain.The computational time for the LSWE with an EBC is slower than the one with LVBM and an EBC, because the internal time step of the 805 ode45 time step routine in MATLAB required a smaller time step ∆t (compared to the LVBM) to preserve the stability.
The shoreline movement of our result compare well with the one of Choi et al. (2011).We can see the comparison in Fig. 8.The solution of Choi et al. (2011) gives a higher pre-810 diction for the shoreline, but it cannot follow the subsequent small positive wave.It may be caused by neglecting the reflection wave and nonlinear effects in their formulation.We also compare the free-surface profile for several time steps in Fig. 9.The implementation of the hard-wall boundary condi-815 tion at x = B in the method of Choi et al. (2011) causes inaccuracies in the prediction of the point-wise wave height in the entire domain.In this case, the effect of reflected waves for the shoreline movement prediction is small, but it may become important when a compound of waves arrives at the 820 coastline.

Periodic wave
Using the same bathymetry profile in the previous case, we influx a periodic wave at the right boundary (x = L) with the profile: in which A = 0.05 m is the amplitude and period T = 20 s.
A smoothened characteristic function until t = 10 s is used in influxing this periodic wave.We use uniform spatial grid ∆x = 1 m in the simulation area and ∆x = 0.015 m in the model area for the numerical solution of the NSWE.
As in the previous case, we also choose the seaward boundary point at h 0 = 10 m at the toe of the slope, that is at x = B = 50 m.Thus, the simulation area is for x ∈ [50, 250] m and the model area for x ∈ [−5, 50] m.Comparisons between these two simulations at several time steps can be seen in Fig. 10 (left: LSWE, right: LVBM).We can see in the comparison that the wave is slightly dispersed in the LVBM.The dispersion ratio at the simulation area is given by µ 2 = 0.0986 < 1, which is less dispersive than the previous case.In Fig. 11, the shoreline movement caused by the periodic wave as well as the paths of characteristic curves forming the shoreline are shown.Observing the results of this case, we can conclude that the EBC technique can deal  robustly with consecutive interactions between the incoming and the reflected wave.
For simulation till the physical time t = 80 s, the computational time for the coupled numerical solutions in both domains is 1.83 times the physical time for the LSWE and 2.01 times the physical time for the LVBM.While the computational time of the simulation using an EBC only takes 0.07 times the physical time for the LSWE and 0.06 times 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 domain simulation.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 amplification due to resonance effects.We follow this bathymetry profile with tan α = 0.0036, tan β = 0.0414, h 0 = 2500m, and h 1 = 200m.
These choices are roughly characterizing the Indian coast bathymetry (Neetu et al., 2011).The EBC point is located at the edge of the last beach slope.We influx periodic wave (64) 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 runup 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 calculate the amplification of periodic waves.In our result, the NSWE model is employed in the last beach slope region.The period of 875 this wave is approximately 2 hours and it coincides with the observed tsunami at the Makran coast, according to Neetu et al. (2011).In nature, one would not expect a tsunami of monochromatic wave train.The investigation of Stefanakis et al. (2011) for the October 25, 2010 Mentawai Islands tsunami 880 showed that the period of the dominant mode of the incident wave is within the resonant regime, and it explained the fact that the highest run-up is not driven by the leading and highest wave.

Simulation using simplified Aceh bathymetry 885
The bathymetry near Aceh, Indonesia, is displayed in Fig. 14. Figure 14a concerns bathymetry data from GEBCO, with a zero value on land.Figure 14b    The location of the EBC point is also determined from Eq. ( 62).For = 0.02 1, the linear model is valid for h 0 25.1 m.Hence, we choose the EBC point at depth h 0 = 41.4 m, which is located at x = B = 12.4 km.Thus, the simulation area is for x ∈ [12.4, 162.4] km, where we follow the real bathymetry of Aceh to calculate the wave propagation.
It is coupled with the model area for x ∈ [−8.6, 12.4] km, where a uniform slope with gradient γ = 1/300 is used to calculate the reflection and shoreline position.We use an irregular grid according to the depth with ratio h 0 /h as the decrease of the wavelength when traveling from a deep region with depth h to a shallower region with depth h 0 in linear wave theory.The grid size used in the simulation area is ∆x = 305 m at the shallowest point near x = B.This choice of spatial resolution is fairly close to a tsunami numerical simulation (Horrillo et al., 2006 use ∆x = 100 m offshore and ∆x = 10 m onshore in one dimensional simulations).
For the numerical solution of the NSWE in the model area, a uniform grid with ∆x = 3 m is used.In Fig. 15, we show the initial profile.Comparisons between these two simulations at several time steps can be seen in Fig. 16.In this case, the wave elevation measured at B has been deformed from its initial condition due to the reflection from the bathymetry before entering the model area, see Fig. 16a and b.We hardly see any differences between the LSWE and LVBM simulations because the wavelength is much larger than the depth.The dispersion ratio at the initial condition is given by µ 2 = 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 shoreline 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 with a 3.3 m runup 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 reduce the computational times up to 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 x s (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, first breaking is approximately taking place at t = 45 min.The corresponding free-surface profiles for several times before and after breaking are shown in Fig. 19.

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 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 975 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.

980
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 also have validated our approach with the 985 laboratory experiment of Synolakis (1987) for the run-up of 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 990 good prediction of the wave run-up as well as the wave reflec-tion, 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 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.
An extension of this EBC technique to the case where the NSWE model is used both in the simulation and model area follows directly from the variational methodology.The analytical benchmark for this case is provided by Carrier et al. (2003) and Kânoglu (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 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ânoglu and Synolakis (1998) can be obtained by requiring the normal velocity at the shoreline wall boundary to be zero.Another characteristic for the outgoing or reflected waves must be derived (either for the LSWE or NSWE model), but the coupling between the numerical and analytical model remains the same as has been derived in this article.

Finite volume implementation
The conservative form of NSWE is given by and the topographic term s s = −gh db/dx 0 . (A3) The system (A1) is discretized using a Godunov finite volume scheme.First the domain [A, B], with some fixed A < x s (t) is partitioned into N grid cells with grid cell k occupying x k− 1 2 < x < x k+ 1 2 .The Godunov finite volume scheme is derived by defining a space-time mesh with ele-1040 ment x k− 1 2 < x < x k+ 1 2 and t n < t < t n+1 and integrating Eqs.(A1) over this space-time element which is a forward Euler explicit method.
To ensure that the depth is non-negative and that the steady 1060 state of a fluid at rest is preserving, the approach of Audusse ( 2004) is used.The numerical flux F is then defined as where the interface values are given by with the waterdepths h (k+ 1 2 ) − and h (k+ 1 2 ) + chosen as follows, to ensure non-negativity of these depths The discretization of the shallow water equations thus reads The HLL flux (Harten et al., 1983;Toro et al., 1994) is used as the numerical flux.It is given by The wave speed S L and S R are approximated as the smallest and largest eigenvalue 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 Appendix B

Coupled model
The (continuous Galerkin) finite element implementation of LSWE or LVBM uses linear polynomials for approximately solving φ, ψ, and η.While the finite volume implementation for NSWE approximates h and u with a constant value.Since u = ∂ x φ, the velocity of the two models are approximated with the same order of polynomials.By coupling both models, in the simulation area we can rewrite Eq. ( 52) as

Fig. 1 .
Fig. 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.
t), where η = h − h b is the wave elevation and h = h(x, t) the total water depth above the bathymetry b = −h b (x) with h b (x) the rest depth.Time runs from t ∈ [0, T ]; partial derivatives are denoted by ∂ t , et cetera, the gradient in the vertical plane as ∇ = (∂ x , ∂ z ) T and the acceleration of gravity as g.
at x = B = 50 m.Therefore, we divide the domain into the simulation area for x ∈ [50, 250] m and the model area for x ∈ [−5, 50] m.

Fig. 3 .
Fig. 3. Run-up of a solitary waves over a canonical bathymetry at 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 and dotted-dashed lines are the coupling of LSWE and NSWE model respectively, and the symbols are laboratory data of Synolakis (1987).

Fig. 4 .
Fig. 4. The shoreline movement of a solitary wave introduced in Synolakis (1987).The LSWE model coupled to the NSWE is shown by the dashed line, while the solid one is the shoreline movement of the LSWE model with an EBC implementation.

Fig. 5 .
Fig. 5.A solitary wave initial condition for the NSWE (dotteddashed line) coupled to the linear model (dashed line), and the linear model with the EBC implementation (solid line) at x = 50 m.The solid and dashed lines are on top of another.

Fig. 6 .
Fig. 6.Free-surface profiles of solitary wave propagation are shown for the coupled linear model (left: LSWE, right: LVBM) with the NSWE (dashed and dotted-dashed lines), and for the linear model with an EBC implementation (solid line), at times (a) t = 10 s, (b) t = 20 s, (c) t = 30 s, (d) t = 40 s.The solid and dashed lines are on top of another.

Fig. 7 .
Fig. 7.The shoreline movement of a solitary wave for the linear model (a: LSWE, b: LVBM) coupled to the NSWE (dashed line) and the linear model with an EBC implementation (solid line).Paths of the first-order characteristic curves are shown by the thin lines.

Fig. 8 .
Fig. 8.Comparison of the shoreline movement of Choi et al. (2011) (dashed line) and the LSWE with an EBC simulation (solid line) for solitary wave case.

Fig. 9 .
Fig. 9. Free-surface profiles of solitary wave propagation are shown for the coupled LSWE with the NSWE (dashed and dotted-dashed lines), for the LSWE with an EBC implementation (solid line), and for the LSWE with the method of Choi et al. (2011) (solid line with 'o' marker) at times (a) t = 10 s, (b) t = 20 s, (c) t = 30 s, (d) t = 40 s.The solid and dashed lines are on top of another.

Fig. 10 .
Fig. 10.Free-surface profiles of periodic waves are shown for the coupled linear model (left: LSWE, right: LVBM) with the NSWE (dashed and dotted-dashed lines), and for the linear model with an 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.

Fig. 11 .
Fig. 11.The shoreline movement of periodic waves for the linear model (a: LSWE, b: LVBM) coupled to the NSWE (dashed line) and for the linear model with an EBC implementation (solid line).Paths of the first-order characteristic curves are shown by the thin lines.
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 890 an earthquake of magnitude of Mw 9.1 at the epicenter point (95.854• E, 3.316 • N), shown by the symbol in Figure 14a.Presently, we consider the following initial N -wave profile η(x, 0) = Af (x)/S with f (x) = d dx exp −(x − x0) 2 /w0 2 and S = max (f (x)) , (65) 895 and the initial velocity potential is zero.We take A = 0.4 m, the position of the wave profile x 0 = 107.4km, and a width w 0 = 35 km.

Fig. 14 .
Fig. 14.Bathymetry near Aceh (a) and the cross section (b) at (95.0278 • E, 3.2335 • N)-(96.6583• E, 3.6959 • N).The solid line concerns the bathymetry data and the dashed line concerns the approximation used in the simulations.

Fig. 15 .
Fig. 15.The initial condition for the Aceh case is shown for the linear model coupled to the NSWE (dashed and dotted-dashed lines) and for the linear model with an EBC implementation (solid line).The solid and dashed lines are on top of another.

Fig. 16 .Fig. 17 .
Fig. 16.Free-surface profiles of Aceh simulations with the linear model (left: LSWE, right: LVBM) coupled to the NSWE are shown by the dashed and dotted-dashed lines and of simulations for a linear 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.

Fig. 18 .
Fig. 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 approximately at t = 45 min.

Fig. 19 .
Fig. 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 another.
cells, we define the mean cell averageU k = U k (t) as U k (t) :h k = x k+ 1 2 − x k− 1 2 .The function U k is 1050piecewise constant in each cell.A numerical flux F is defined to approximate the flux f