# A theoretical model of tsunami runup in Chile based on a simple bathymetry

M. A. Fuentes, J. A. Ruiz, Armando Cisternas

<figure><img src="image_1.png" /></figure>

## To cite this version:

M. A. Fuentes, J. A. Ruiz, Armando Cisternas. A theoretical model of tsunami runup in Chile based on a simple bathymetry. Geophysical Journal International, 2014, 196 (2), pp.986-995. hal-00947615

HAL Id: hal-00947615

https://hal.science/hal-00947615

Submitted on 5 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A theoretical model of tsunami runup in Chile based on a simple bathymetry

Mauricio A. Fuentes,<sup>1,2</sup> Javier A. Ruiz<sup>1</sup> and Armando Cisternas<sup>1,3</sup>

Accepted 2013 October 10. Received 2013 October 7; in original form 2013 July 11

## SUMMARY

We present a theoretical solution for the water wave elevation along the shoreline and a closed formula is derived for the runup distribution in the 2+1 dimensional (2 spatial + 1 temporal) linear case. An integral transform method is used to solve the linear shallow water wave equations for a sloping beach model extended along the coastline. At large regional scale, this simple bathymetry model gives a good approximation for the seafloor geometry along the Chilean subduction zone. In the bidimensional problem, for special cases, we found an analytic expression for the runup distribution along the shoreline that includes the obliquity of the incident incoming wave approaching the coast. Our solution establishes that obliquity can reduce maximum runup in important ways. It suggests that the distribution along the shoreline is determined mainly by the maximum unidimensional runup modulated, by a function controlling the horizontal extent of the incident wave parallel to the shoreline. This function was introduced to include the physical expected boundedness of the leading wave approaching the shore. Our theoretical solution was applied to model the runup generated by the 2010, $M_w$ 8.8 Maule mega-thrust earthquake. The model fits rather well the runup to first-order. The theoretical solution proposed in this study also provides new insights to the near-field tsunami discriminant formula proposed by Okal & Synolakis.

**Key words:** Fourier analysis; Non-linear differential equations; Tsunamis; Continental margins; convergent.

## 1 INTRODUCTION

The problem of tsunamis is important, particularly in Chile, due to the subduction process between the Nazca and Sudamerican plates. Many large tsunamigenic earthquakes have been generated over the seismogenic contact interface, such as, the 1960, $M_w$ 9.5 Valdivia (e.g. Kanamori & Cipar 1974; Cifuentes 1989) and the 2010, $M_w$ 8.8, Maule mega-thrust events (e.g. Lay et al. 2010; Fritz et al. 2011). Also, the Peruvian and northern Chile subduction zones are under a high tsunami hazard (Okal et al. 2006).

The maximum runup is possibly the most important tsunami parameter to quantify the tsunami hazard at which the population is exposed along a coastal region. Several numerical and theoretical models have been developed to better understand the runup and the generation, propagation and inundation processes (Kajiura 1970; Synolakis 1987; Shuto 1991; Liu et al. 1994; Tadepalli & Synolakis 1994, 1996; Zhang & Wu 1999; Yamazaki et al. 2009).

In particular, Okal & Synolakis (2004) proposed an empirical formula to characterize the runup distribution along the shoreline, but this expression was derived from numerical modelling in order to discriminate the nature of the source of the runup of near-

field tsunamis; by comparing the runup from seismic dislocations sources and offshore underwater landslides. This fitting model is a zero-order approximation and the simple runup formula does not allow the introduction of a different function distribution along the coast.

The 1 spatial + 1 temporal dimensional runup problem (hereinafter referred to as 1+1) in a sloping beach has been studied theoretically in the linear and non-linear cases (Carrier & Greenspan 1958; Synolakis 1987), and it is very well understood. However, in the 2 spatial + 1 temporal dimensional case (hereinafter referred to as 2+1), due to the huge complexity involved on the governing equations there are not many analytical results and there are fewer solutions for the runup distribution, $R(y)$, along the coastline (Okal & Synolakis 2004).

Along the Chilean subduction zone, large tsunamigenic earthquakes are generated, and the runup problem becomes important. At large geographical scale, the seafloor geometry along the Chilean coast can be approximated by a sloping beach model. This simple bathymetry can also be extended along the north-south direction, because the coastline is roughly linear and keeps parallel to the trench axis.

<footnote n="1"><sup>1</sup>Department of Geophysics, Faculty of Physical and Mathematical Sciences, University of Chile, Santiago, Chile. E-mail: mfuentes@dim.uchile.cl</footnote>

<footnote n="2"><sup>2</sup>Department of Mathematical Engineering, Faculty of Physical and Mathematical Sciences, University of Chile, Santiago, Chile</footnote>

<footnote n="3"><sup>3</sup>Institut de Physique du Globe. University of Strasbourg, Strasbourg, France</footnote>
The aim of this study is to derive theoretically an analytical formula for the runup distribution along the coastline in the 2+1 linear case. In this paper we give a theoretical study to solve for the water wave elevation and to obtain the runup distribution along the coast at the shoreline. In the next sections, an overall review of the theory is presented, including the theoretical methodology used to solve the shallow water wave equations for the 2+1 linear problem. We finalize showing a particular solution for the bidimensional runup distribution, that includes the obliquity of the incident wave. As a simple application, we model the runup generated along the Chilean coast by the near-field tsunami due to the 2010, M<sub>w</sub> 8.8, Maule mega-thrust earthquake.

## 2 BASIC THEORY

Considering the water like a very low viscosity fluid, the Euler equations governing the fluid motion are,

$$ \rho \left( \frac{\partial}{\partial t} + \mathbf{u} \cdot \nabla \right) \mathbf{u} + \nabla p + \rho g \hat{z} = 0, \quad (1) $$

where $\mathbf{u} = (u, v, w)$ is the velocity field, $p$ the fluid pressure, $\rho$ the density (assumed constant), $g$ the gravity acceleration, $\hat{z}$ the unit vector corresponding to the z-direction and $\nabla$ is the gradient operator. Assuming that the water is incompressible, the mass conservation equation is,

$$ \mathrm{div}(\mathbf{u}) = 0. \quad (2) $$

Fig. 1 shows the cartesian coordinate system used. Neglecting the dependence of the horizontal components of the velocity field with depth, the kinematic boundary conditions at the sea surface $\eta$ over the datum level ($z = 0$) and at the sea bottom $h$ are,

$$ w(x, y, \eta) = \eta_t + u(x, y, \eta)\eta_x + v(x, y, \eta)\eta_y \quad (3) $$

$$ w(x, y, -h) = -u(x, y, -h)h_x - v(x, y, -h)h_y. \quad (4) $$

The governing eq. (1) combined with (2), (3) and (4) allow the derivation of the classical dimensional non-linear shallow water wave equations system,

$$ \begin{align} u_t + u u_x + v u_y + g \eta_x &= 0 \notag \\ v_t + u v_x + v v_y + g \eta_y &= 0 \notag \\ \eta_t + ((\eta + h)u)_x + ((\eta + h)v)_y &= 0, \tag{5} \end{align} $$

where $u$ and $v$ represent the depth-integrated velocity field components.

The water elevation $\eta$ is estimated by solving the set of eqs (5) when an incoming leading wave approaches the shore (shallow water) and evaluating it at the shoreline, it is possible to obtain the runup.

<figure><img src="image_3.png" /><figcaption>Figure 1. Scheme of variables and coordinate system.</figcaption></figure>

## 2.1 1-D propagation problem

For completeness, we present here the 1+1 case for the runup problem in a sloping beach model. Assuming that the incident waveform is independent of the y-coordinate, the non-linear set (5) is reduced to,

$$ \begin{align} u_t + u u_x + g \eta_x &= 0 \notag \\ \eta_t + ((\eta + h)u)_x &= 0. \tag{6} \end{align} $$

This non-linear set has been studied in exact form by Carrier & Greenspan (1958) and Synolakis (1987), among others. They used variational methods (Riemann invariants) to find an explicit solution. In particular, Synolakis (1991) found that the linear theory can estimate the runup well enough as compared to the non-linear results.

The set (6) can be linearized and transformed into a one single second order partial differential equation (PDE),

$$ \eta_{tt} - g(h\eta_x)_x = 0, \quad (7) $$

where the bathymetry for a sloping beach model is given by,

$$ h(x) = \begin{cases} \tan(\beta)x & \text{if } x \le x_0 \\ d & \text{if } x > x_0 \end{cases} \quad (8) $$

as shown in Fig. 2. The origin of the coordinate system is located at the shoreline, $x$ is positive seaward, $\beta$ represents the angle of the sloping beach and $d$ is the constant depth of the ocean for $x > x_0$.

### 2.1.1 Time-space Fourier transform formalism

We use the time-space Fourier transform that can be easily extended to solve a problem involving more spatial dimensions.

(i) It can be shown that when $x > x_0$, the general solution for the water elevation is,

$$ \eta(x, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \{ A_i(\omega)e^{i\omega t - ik_0(\omega)x} + A_r(\omega)e^{i\omega t + ik_0(\omega)x} \} d\omega, \quad (9) $$

where the dispersion relation is given by $\omega = k_0(\omega)c_0$, being $c_0 = \sqrt{gd}$ and the amplitudes $A_i$, $A_r$ represent the reflected and incident coefficients, respectively.

(ii) When $0 < x < x_0$, a useful change of variable for eq. (7) is the following: $v = \frac{2}{c_0}\sqrt{x}x_0$, $\zeta(v, t) = \eta(x, t)$.

The next step is to apply the time-space Fourier transform to eq. (7) in order to solve a boundary value problem. After the change of variable, to convert the PDE into a single ordinary differential equation (ODE), the problem can be solved in a straightforward manner. The ODE is shown in (A5) and its solution is,

$$ \bar{\bar{\zeta}}(k, \omega) = \frac{C}{\sqrt{w^2 - k^2}}, \quad (10) $$

where $C = C(\omega)$ is a function independent of $k$ and $\bar{\bar{\zeta}}$ denotes the time-space Fourier transform.

Splitting eq. (10) into real and imaginary parts, it is possible to write a general solution in the $(k, \omega)$ Fourier space as,

$$ \bar{\bar{\zeta}}(k, \omega) = \frac{C}{\sqrt{w^2 - k^2}} H(|w| - |k|) - i \frac{C}{\sqrt{k^2 - w^2}} H(|k| - |w|) \quad (11) $$

for a real $C(\omega)$ and where $H(\cdot)$ is the Heaviside step function. Inverting the spatial part of the previous equation and using integral forms of the cylindrical Bessel functions of first and second kind
<figure><img src="image_4.png" /><figcaption>Figure 2. Sloping beach model.</figcaption></figure>

(Abramowitz & Stegun 1964), the two general solutions for the PDE in terms of the (v, t) variables are,

$$\zeta_1(v, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} B(\omega)e^{i\omega t} J_0(\omega v) \, d\omega$$

$$\zeta_2(v, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} D(\omega)e^{i\omega t} Y_0(\omega v) \, d\omega,$$

where B(ω) and D(ω) are functions independent of v.

Finally, because the wave height must remains finite at the shoreline, the bounded solution over the sloping beach, in terms of the original space variable x, can be written as,

$$\eta(x, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} B(\omega)e^{i\omega t} J_0(2k_0\sqrt{x\bar{x}_0}) \, d\omega. \quad (12)$$

After imposing at $x = x_0$ the continuity conditions of the solution and its slope at the intersection between the two bathymetry domains, one obtains exactly the same known kernel $B(\omega)$ (Synolakis 1987),

$$B(\omega) = \frac{2}{J_0(2k_0x_0) + iJ_1(2k_0x_0)} e^{ik_0x_0} A_i(\omega),$$

named also the amplification factor for the transmitted solution (see Appendix A, eq. A11).

### 2.1.2 Runup for the solitary wave solution

Taking the most regular solution, which means fast decay of the solution and their derivatives at infinity, from the Boussinesq equation (Boussinesq 1871), one obtains the classic solitary wave profile,

$$\eta(x, t) = H\mathrm{sech}^2(\gamma[x - x_1 \pm ct]), \quad (13)$$

where $H$ is the initial maximum amplitude, $x_1$ is the centred initial location, $\gamma = \frac{1}{d}\sqrt{\frac{3H}{4d}}$ is a factor related to the wavelength and $c = c_0\sqrt{1 + \frac{H}{d}}$ is the solitary wave speed, which differs from the linear tsunami wave speed $c_0 = \sqrt{gd}$, because the non-linear and dispersive effects are compensated in the non-linear Boussinesq equation (Russel 1844; Boussinesq 1871). The Fourier spectrum for this solitary wave solution is given by,

$$\Phi(\omega) = \frac{4d^3 \pi}{3c} e^{ik_c(x-x_1)} k_c \mathrm{cosech}(\alpha k_c) \quad (14)$$

with $k_c = \frac{\omega}{c}$ (Synolakis 1987).

Synolakis (1987) has shown, using linear superposition and complex contour integration, that the runup ($\mathcal{R}$) formula for the initial value problem, for an incident solitary wave can be estimated by,

$$\mathcal{R} = 2.831H \left(\frac{H}{d}\right)^{\frac{1}{4}} \sqrt{\cot(\beta)}, \quad (15)$$

which he termed the runup law. The computation in the linear theory is done by maximizing over the time the water elevation at the shore (vertical wall). Synolakis (1987) showed that this approach to first-order is quite similar to the runup predicted by the non-linear theory, what actually takes into account effects overland. Ward (2010) has a different known formula to estimate the runup: $\mathcal{R} = (H)^{\frac{1}{4}}(d)^{\frac{3}{4}}$, but the theory he used is different from the one we are using in this work.

However, by applying the boundary value condition at $x = x_0$, the runup found in this study using the time-space Fourier transform is,

$$\mathcal{R} = 2.831H \left(\frac{H}{d}\right)^{\frac{1}{4}} \left(1 + \frac{H}{d}\right)^{\frac{1}{4}} \sqrt{\cot(\beta)}. \quad (16)$$

One can note that (16) is slightly different from (15), due to the extra factor $(1 + H/d)^{1/4}$ that appears when solving a boundary value problem, because the solitary wave speed $c$ is greater than the linear tsunami wave speed $c_0$. However, for a wide range of values of $H$ and $d$, keeping $H/d$ small ($H/d$ is usually less than 0.1 for non-breaking waves), eqs (15) and (16) produce about the same estimation of $\mathcal{R}$.

## 2.2 2-D tsunami propagation problem

In this section, we study the 2+1 tsunami wave propagation. We take the same bathymetry model (Fig. 2), but now extended in the y-direction. The seafloor depth can be written as,

$$h(x, y) = \begin{cases} d & \text{if } x > x_0, y \in \mathbb{R} \\ x \tan(\beta) & \text{if } x \le x_0, y \in \mathbb{R}. \end{cases} \quad (17)$$

It is straightforward to show that the linearized form of set (5) is,

$$\eta_{tt} - g[(h\eta_x)_x + (h\eta_y)_y] = 0. \quad (18)$$

Starting with PDE (7), we can do a similar analysis using the time-space Fourier transform applied to the (y, t) variables.

(i) For $x > x_0$, the solution for the water elevation is,

$$\eta(x, y, t) = \frac{1}{(2\pi)^2} \iint_{\mathbb{R}^2} (A_i e^{ik_x x} + A_r e^{-ik_x x}) e^{i\omega t - i\kappa_y y} dk_y d\omega, \quad (19)$$
with $k_x = \sqrt{k_0^2 - k_y^2}$, $k_0 = \frac{\omega}{c}$, and where $k_x$, $k_y$ are the wavenumber vector components in the x- and y-directions.

(ii) For $0 < x < x_0$, the bounded solution is,

$$ \eta(x, y, t) = \frac{1}{(2\pi)^2} \iint_{\mathbb{R}^2} B(k_y, \omega) X(x) e^{i\omega t - ik_y y} dk_y d\omega, \quad (20) $$

where the function $X(x)$, is $X(x, k_y, \omega) = e^{-k_y x} M(\mu, 1, 2k_y x)$ (see eq. B5), $\mu = \frac{1}{2}(1 - \frac{k_0^2 x_0}{k_y^2})$ and $M(a, b, z)$ is the Kummer confluent hypergeometric function. Similar to the 1+1 case, $B$ is the amplification transmitted factor, but is now dependent on two variables. See Appendix B for more details.

Following the same steps as in the 1+1 case, we apply the matching conditions at $x = x_0$ to retrieve the $B$ factor in terms of the incident coefficient, $A_i$. Then, the water elevation solution on the sloping beach, evaluated at the shoreline, is

$$ \eta(0, y, t) = \frac{1}{(2\pi)^2} \iint_{\mathbb{R}^2} \frac{2\Phi(k_y, \omega)e^{i\omega t - ik_y y}}{X(x_0) - \frac{i}{k_x} \frac{\partial X}{\partial x}(x_0)} dk_y d\omega, \quad (21) $$

where the boundary condition at $x_0$ is given by $\bar{\eta}(x_0, k_y, \omega) = \Phi(k_y, \omega) = e^{ik_0 x_0} A_i(\omega)$, and $\Phi$ represents the incident wave Fourier spectrum.

Eq. (21) is the exact theoretical solution for the water wave elevation along the shoreline for any given incident wave.

## 3 RESULTS FOR 2+1 LINEAR THEORY

One issue when computing expression (21) to retrieve a closed analytical solution for $\eta(0, y, t)$, for a given $\Phi(k_y, \omega)$, is that $k_x$ can be real or complex. It is now more difficult to use the Cauchy's residue theorem when analysing the poles on the complex plane of the kernel. To overcome this difficulty, one technique has been proposed and used by Carrier & Noiseux (1983) and Zhang & Wu (1999), and it will be used later.

In the following sections, we obtain an explicit analytical solution for the runup in the 2+1 case, where an incident solitary wave is analysed.

### 3.1 The runup of a 2+1 incident solitary wave of infinite horizontal extent

Taking a 2+1 incident solitary wave with an incoming angle $\theta$ with respect to the y-direction,

$$ \eta(x, y, t) = H \operatorname{sech}^2 \{\gamma [\sin(\theta)(x - x_1) - \cos(\theta)y + ct]\}. \quad (22) $$

Note that to evaluate the water elevation, one needs to compute the time-space Fourier transform of eq. (22). It is possible to compute explicitly the Fourier spectrum $\Phi(k_y, \omega)$, because it remains the same structure that the 1+1 solitary wave, its expression is shown in Appendix B (see eq. B7).

We verified, under these hypothesis, that the runup along the shoreline is independent of $y$, which is not surprising because the 2+1 incident solitary wave has the same height $H$ all along the transverse direction. The water elevation at the shoreline as a function of $y$ is,

$$ \eta(0, y, t) = \frac{4d^3}{3} \int_{-\infty}^{\infty} \frac{k_c \operatorname{cosech}(\alpha k_c)}{X(x_0) - \frac{i}{ak_c} X'(x_0)} \times e^{ik_c(\sin(\theta)(x_0 - x_1) + \cos(\theta)y + ct)} dk_c. \quad (23) $$

From eq. (23), we find that $\eta(0, y, t)$ can be written as $\eta(0, 0, t^*)$, where $t^* = t + \frac{\cos(\theta)y}{c}$. In order to obtain the runup $\mathcal{R}(y)$, we have to maximize $\eta$ over $t^*$. This transformation of $t$ to $t^*$ allows to see that the runup is independent of $y$, but the time of the maximum is not. The factor $e^{ik_c t^* c}$ acts as a phase change and the kernel of the integral does not depend on $y$.

Additionally, a simple way to verify the solution (21) is to check, if the runup in the 2+1 case contains correctly the result found in 1+1, when an incident 2+1 solitary wave is analysed. Here, we take the limit case when $\theta \to \frac{\pi}{2}$, namely normal incidence to the beach. This produces exactly the same runup as the 1+1 linear theory,

$$ \mathcal{R}(y) = 2.831H \left(\frac{H}{d}\right)^{\frac{1}{4}} \left(1 + \frac{H}{d}\right)^{\frac{1}{4}} \sqrt{\cot(\beta)}, $$

and it is constant along the shoreline. Details of the derivation can be found in Appendix B.1.

### 3.2 The particular solution

It is not straightforward to obtain a closed analytical general form for the runup from eq. (21). In order to retrieve a particular solution, we restart from the 2+1 linear theory (Appendix B). But this time, we consider the linear shallow water equations with the following transformation,

$$ s = y - \sec(\theta)c_0 t, $$

used by Carrier & Noiseux (1983) to study reflected waves due to oblique incident tsunamis. Here, $\theta \in [0, \pi]$ represents the incident angle of the wave front with respect to the shoreline. This transformation allows to simplify the partial derivates with respect to time and $y$ variables, reducing the problem to only two variables, $x$ and $s$.

The transformed equations in the Fourier space for the two different bathymetry domains are,

(i) $x > x_0$.

$$ \hat{\eta}_{xx} + \tan^2(\theta) k_s^2 \hat{\eta} = 0, $$

which is a simple ODE, and $k_s$ is the wavenumber associated to the $s$ variable. $\hat{\eta}$ represents the 1-D space Fourier transform.

(ii) $0 < x < x_0$.

$$ x\hat{\eta}_{xx} + \hat{\eta}_x + (k_s^2 \sec^2(\theta)x_0 - k_s^2 x)\hat{\eta} = 0, $$

which is a type of a hypergeometric differential equation.

We apply the matching condition at $x = x_0$ to express the transmitted amplification coefficient in terms of the incident one. As in the 1+1 case, it is possible to compute $\hat{\eta}(x, k_s)$ on the sloping beach region. Then, applying the inverse Fourier transform, one obtains the solution $\eta(x, y, t)$ evaluated at the shoreline,

$$ \eta(0, y, t) = \frac{1}{\pi} \int_{-\infty}^{+\infty} \frac{\Phi(k_s) e^{-ik_s s}}{M_0(k_s) - iN_0(k_s)} dk_s, \quad (24) $$

The functions $M_0$, $N_0$ in this kernel are given in terms of the hypergeometric function M (Appendix B, eq. B4),

$$ M_0(k_s) = e^{-k_s x_0} M(\mu, 1, 2k_s x_0) \quad \text{and} \\ N_0(k_s) = \cot(\theta) e^{-k_s x_0} [M(\mu, 1, 2k_s x_0) - 2\mu(-k_s)M(\mu, 2, 2k_s x_0)]. $$

In this case, $\mu(k_s) = \frac{1}{2}[1 - k_s \sec^2(\theta)x_0]$, and $\Phi$ is given by the boundary condition at $x_0$ as before, and corresponds to $\Phi(k_s) = \hat{\eta}(x_0, k_s)$.

It is possible to show that the functions $M_0$ and $N_0$, being depending on $k_s$, are even and odd functions, respectively. This implies that actually the integral (24) is real valued.
### 3.2.1 The Particular solution: damped incident wave

The results shown in Section 3.1 show that the runup due to an incident solitary wave is constant along the shoreline. It is because the incoming wave is infinitely extent and has the same shape all along the transverse direction, it means that keeps the same height and wavelength when it arrives at the distance $x_0$ from the coast, no matter how far seaward is the incoming wave from $x_0$.

In a real tsunami, the travelling waves loss energy, there are also directivity effects from the tsunami source (Kajiura 1970), and dispersive effects. Then, it is expected that the incoming wave to the shore is bounded along the transverse direction and the wave height is not constant. One type of finite crest wave was explored by Kânoğlu et al. (2013) in an analytic treatment.

The aim of this section is to look for a runup formula that includes physical properties observed not only on theoretical and numerical models, but also on real tsunami propagation, when modelling the incident wave.

In order to extend the 1-D problem into a 2+1, we propose an initial damped solitary wave given by,

$$ \eta(x, y, t) = H \operatorname{sech}^2(y[\sin(\theta)(x - x_1) - a \cos(\theta)s]) f(p) \quad (25) $$

with $p = a \sin(\theta)y - \cos(\theta)(x_1 - x)$ the distance along the line of the wavefront. Note there is a change of scale between $x$ and $y$ due to the factor $a = \sqrt{1 + \frac{H}{d}}$, but it is imperceptible because $a \approx 1$. $f: \mathbb{R} \to [0, 1]$ is such that $\lim_{|p| \to \infty} f(p) = 0$. This assumption about the function $f$ is taken because an infinitely wide incoming wave cannot be realistic, hence the function $f$ introduces the natural decay of the wave at its ends due to the intrinsic finite characteristic of the tsunamigenic seismic source.

Fig. 3 illustrates an example of an initial incoming wave (centred at $x = x_1$ and $y = 0$) for this model in a 2-D sloping beach bathymetry. The $f(p)$ function sketches the decay that one can impose to the initial wave along the transverse direction.

Following the same formalism than Synolakis (1987), we perform contour integration, so that zeros in the denominator (poles) on the integrand of (24) in the complex plane are important. To do it, the asymptotic expansion of the functions $M_0$ and $N_0$ (Abramowitz & Stegun 1964; Zhang & Wu 1999) are used to approximate,

$$ M_0(k) - iN_0(k) = [kx_0\pi \tan(\theta)]^{-\frac{1}{2}} e^{ikx_0v(\theta) - i\frac{\pi}{4}} + O(|k|^{-1}), \quad (26) $$
with $v(\theta) = (\frac{\pi}{2} - \theta) \sec^2(\theta) + \tan(\theta)$.

Except for some constant factors, the kernel of integral (24) has the same asymptotic expansion than the one studied in the 1+1 dimensional case. This allows us to try a similar mathematical proof to the one given by Synolakis (1988). Both kernels remain close enough to first-order. By virtue of the Rouché's theorem, also both kernels have the same number of roots inside a particular region in the complex plane.

Performing contour integration and the series calculation by the Cauchy's residue theorem, we obtain the runup formula, as

$$ \mathcal{R}(y) = f[a \sin(\theta)y - (x_1 - x_0) \cos(\theta)] \sqrt{\sin(\theta)} \mathcal{R}_0, \quad (27) $$

where $a = \sqrt{1 + \frac{H}{d}}$ and $\mathcal{R}_0$ is the 1-D runup generated by an incident solitary wave given by expression (16).

Formula (27) suggests that the runup distribution along the coastline is now given by the transmitted initial waveform to the sloping beach, but is modulated by the 1-D runup and the factor $\sqrt{\sin\theta}$, being maximum in the case of normal incidence, as one might expect.

### 3.2.2 Runup comparison for different model parameters

In the following, we present the results obtained for the runup distribution along the $y$-direction for different model parameters. Let's point out that, Synolakis (1987) studied the effect of $\beta$, $d$ and $H$ on the 1-D runup $\mathcal{R}_0$. Then, in this study, we will focus on the analysis of the initial waveform $f(p)$ and the incident angle.

<figure><img src="image_5.png" /><figcaption>Figure 3. Sketch of the initial incoming wave for the particular solution analysed in the 2+1 case, showing the decay function, $f(p)$, along the transverse direction.</figcaption></figure>
We set-up a sloping beach model for a typical Chilean subduction zone, by approximating the seafloor bathymetry assuming a mean ocean depth of $d = 4$ km, and with a characteristic distance from the trench axis to the coastline of about 115 km, which gives a sloping beach angle of, $\beta = 2^\circ$. The initial height of the incident wave was fixed to 10 m at the constant depth domain.

We choose three different $f(p)$ functions: Gaussian, Lorentzian and one bimodal function, which is a combination of two Lorentzian functions. All these functions are built using a characteristic length $L$, which is associated with the length of the initial wave. We set $L = 450$ km, that attempts to capture the fault length of a large earthquake (let's say, $M_w \approx 8.8$) that could be generated along the Chilean subduction zone.

The functions used are,

(i) $f_1(p) = e^{-(\frac{2p}{L})^2}$

(ii) $f_2(p) = \frac{1}{1+(\frac{2p}{L})^2}$

(iii) $f_3(p) = 0.9713\{\frac{1}{1+(\frac{8p}{L})^2} + \frac{1}{2[1+(\frac{8p}{L}-4)^2]}\}$

that by construction their maximum value is one.

The initial locations of the incident wave analysed are,

(i) $x_1 - x_0 = 20$ km, and,

(ii) $x_1 - x_0 = 200$ km.

Note that the closest distance to $x_0$ corresponds to an initial wave generated close to the trench.

Fig. 4(a) shows the runup distribution along the coastline for different incident angles using a Gaussian function. The maximum runup is reached for normal incidence ($\theta = 90^\circ$), and it decreases when $\theta$ diminishes. Because of the symmetry with respect to $\theta$, we considered angles bounded between 0° and 90°. The runup shape for each $\theta$ is preserved, but the peak value shifts to the right and the shape gets wider when $\theta$ decreases. This behaviour is related to

the obliquity of the incident wave to the sloping beach and because the runup maximum time as a function of the y-direction becomes steeper (Fig. 4d). The maximum runup time has been cut to the zero threshold because negative times have no physical sense, and it occurs when the incident wave is nearby the toe beach and the characteristic length is comparable to the distance from $x_0$ to the shoreline.

Fig. 4(b) shows the results obtained in the case of a Lorentz function. One can observe a similar behaviour as in the Gaussian function case. The main difference is that the runup distribution along the coastline attenuates slower at the ends in comparison to a Gaussian function. In the case of a bimodal function [$f_3(p)$], one obtains similar conclusions as in the two previous cases (Fig. 4c).

Fig. 5 presents the runup distribution along the y-direction for an incident wave located at $x_1 - x_0 = 200$ km and for a Gaussian, Lorentzian and Bimodal function respectively. We observe a common pattern in the distribution of $R(y)$, in the three cases and for different incident angles, as the maximum runup value decreases, the peaks shift to the right and the shape becomes wider when $\theta$ decreases.

If we compare the runup distribution for a same $f(p)$ function, but generated for an incident wave located at 20 km and 200 km distance from $x_0$ (Figs 4 and 5), as expected, the maximum runup is the same, the peak shifts to the right and the $R(y)$ shape becomes broader when increasing the distance from $x_0$ for a given incident angle. A similar effect is observed for the other functions $f(p)$, Lorentzian (Figs 4b and 5b) and Bimodal (Figs 4c and 5c). The maximum runup time is the same for the three functions considered (Fig. 5d), but for a given y distance and $\theta$ angle, the maximum runup time is longer in the case $x_1 - x_0 = 200$ km than the $x_1 - x_0 = 20$ km case, because of the distance that the incoming wave has to travel.

<figure><img src="image_6.png" /><figcaption>Figure 4. Runup distribution modelled along the y-direction for different incident angles ($\theta = 20^\circ, 45^\circ, 60^\circ$ and $90^\circ$). The initial incident wave is located at $x_1 - x_0 = 20$ km and the damping functions used are (a) Gaussian, (b) Lorentzian and (c) Bimodal Lorentzian. (d) Maximum runup time as a function of y.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same caption as Fig. 4, but using $x_1 - x_0 = 200$ km.</figcaption></figure>

### 3.2.3 Application to real data: Runup and tsunami generated by the 2010, $M_w$ 8.8 Maule Earthquake

In this study we have derived a theoretical formula for the runup distribution, that varies along the coastline. Our results showed that one can describe the runup with either a simple or bimodal function distribution.

Based on a numerical simulation, Okal & Synolakis (2004) proposed a simple formula that they used to empirically fit and to globally adjust the measured runup of past tsunamis along straight coastlines.

In the following analysis, we attempt to adjust to first-order the runup generated by the 2010, $M_w$ 8.8, Maule earthquake. For near-field tsunamis we considered the co-seismic vertical static displacement at the surface generated by the Maule earthquake. In this application, we used the imaged slip provided by the USGS (G.Hayes, http://earthquake.usgs.gov/earthquakes/eqinthenews/2010/us2010t fan/finite_fault.php) and we computed the vertical static displacement at the free surface (Fig. 6a) using the Okada's formulas (Okada 1985), that allows to compute the static displacement field in an elastic half-space. The uplift is basically characterized by two zones having large amplitudes.

The two distinctive zones with large vertical static displacement suggest the use of a bimodal $f(p)$ function. We used a $f_3$ function and, by trial-and-error, we tuned the function parameters to adjust the observed runup. The $f(p)$ function was constructed as a combination of two Lorentzian functions with different amplitudes and central locations: $f(p) = \epsilon[f_2(bp) + af_2(bp + c)]$, where a, b, c are shape-parameters calculated by trial and error and $\epsilon$ is an amplitude factor to ensure the maximum value of $f$ to one. Fig. 6(b) shows the observed and modelled runup distributions along the south-north direction.

The sloping beach model was approximated by, $\beta = 2^\circ$ and $d = 4$ km, which gives $x_0 \approx 115$ km, a reasonable approximate distance from the Chilean trench to the shore at the latitude of Maule

<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Vertical coseismic static displacement computed for the 2010, $M_w$ 8.8, Maule earthquake, grey star corresponds to the epicentre. (b) Comparison of the observed (grey squares) and modelled (solid black line with $H = 9$ m and dashed black line with $H = 4$ m) runup as a function of the longitude. The adjusted function is, $f(p) = \epsilon[f_2(bp) + af_2(bp + c)]$, with $\epsilon = 0.9372$, $a = 0.6$, $b = 3.5$ and $c = 1000$ km. The observed runup data were downloaded from NOAA's website, http://www.noaa.gov.</figcaption></figure>

(Moscoso et al. 2011). The incident solitary wave was located at $x_1 - x_0 = 20$ km distance with an initial wave height of, $H = 9$ m.

Note that $H = 9$ m adjust best the maximum runup observed (~29 m), however the maximum uplift is in the order of 4 m. If we
use this latter value, the maximum runup modelled is reduced by a factor of 0.36. Note that keeping *d* constant, the scaling factor is given by $(\frac{H_1}{H_2})^4$, where $H_1$ and $H_2$ are two different initial wave heights. The two largest peak-runups can be attributed to several hydrodynamic and local bathymetric effects, such as, non-linearity, constructive interference, local resonance (Yamazaki & Cheung 2011) and focusing of water waves, that are not taken into account in this theoretical work.

## 4 CONCLUSIONS

In this work, we have applied an Integral Transform Method based on a time-space Fourier Transform to compute analytically the water elevation and runup, by solving the linear shallow water wave equations in both, 1+1 and 2+1 dimensions, for a sloping beach bathymetry model.

In the 1+1 case, we have shown that the runup formula is quite similar to the expression derived by Synolakis (1987), but differing by an extra factor which is close to one under the validity range of the shallow water wave equation $(\frac{H}{d} \ll 1)$.

In the 2+1 propagation problem, we obtained a new solution for the water wave elevation for a sloping beach model extended along the coastline. In particular, a theoretical solution is proposed for the water elevation at the shoreline for any incident wave approaching the shore. We have shown that in the case of an infinitely wide incident solitary wave, the runup does not vary along the coastline, which have not physical sense because in a real tsunami the leading wave is bounded.

In order to derive a closed runup formula, we used the Carrier & Noiseux's (1983) transformation to simplify the partial differential equation and a damping function $f(p)$ was introduced to bound the incident wave along the transverse direction. This function allows the inclusion of the boundedness physical condition of the leading wave. Under these assumptions, we obtained a runup formula that varies along the shoreline,

$$ \mathcal{R}(y) = f[a \sin(\theta)y - (x_1 - x_0) \cos(\theta)] \sqrt{\sin(\theta)} \mathcal{R}_0, \quad (28) $$

where $\mathcal{R}_0$ is the 1+1 runup and $f(p)$ is a damping function, referring to the characteristics of the initial profile.

The previous formula describes the runup distribution along the coastline, and it shows that the shape depends mainly on the transverse form of the incident wave approaching the shore. This result has several implications to model, for instance, runup of near-field tsunamis, because it allows to model a more realistic initial incoming waves. The runup distribution modelled for the tsunami generated by the 2010, $M_w$ 8.8, Maule earthquake shows that one can describe rather well to first-order the observed runup amplitudes. Similarly to the near-field discriminant formula proposed by Okal & Synolakis (2004), the runup closed form proposed in this study can improve first order runup predictions for future tsunamis.

## ACKNOWLEDGEMENTS

Special thanks to Jaime Campos who helped with contacts with the French Embassy and LIA Montessu de Ballore, a collaboration program between CNRS-France and University of Chile, for a short visit at Strasbourg, France. Financial support from University of Chile scolarship for short visit to Honolulu, Hawaii is deeply grateful. This work was partially funded by CONICYT and Fonde-

cyt iniciación, Grant number 11100134. We thank Steven Ward and one anonymous reviewer for their constructive suggestions that helped us to improve the manuscript.

## REFERENCES

Abramowitz, M. & Stegun, I.A. eds, 1964. *Handbook of Mathematical Functions*, National Bureau of Standards.

Boussinesq, J., 1871. Théorie de l'intumenscense liquide appelée onde solitaire ou de translattion se propageant dans un canal rectangulaire, Institute de France, Académie des Sciences, *Comptes Rendus*, **72**, 755-759.

Carrier, G.F. & Greenspan, H.P., 1958. Water waves of finite amplitude on a sloping beach, *J. Fluid Mech.*, **17**, 97-110.

Carrier, G.F. & Noiseux, C.F., 1983. The reflection of obliquely incident tsunamis, *J. Fluid Mech.*, **133**, 147-160.

Cifuentes, I.L., 1989. The 1960 Chilean earthquakes, *J. geophys. Res.*, **94**(B1), 665-680.

Fritz, H.M. et al., 2011. Field survey of the 27 February 2010 Chile tsunami, *Pure appl. Geophys.*, **168**, 1989-2010.

Kajiura, K., 1970. Tsunami source, energy and the directivity of wave radiation, *Bull. Earthq. Res. Inst.*, **48**, 835-869.

Kanamori, H. & Cipar, J.J., 1974. Focal process of the great Chilean earthquake May 22, 1960, *Phys. Earth planet. Inter.*, **9**, 128-136.

Kánoğlu, U., Titov, V.V., Aydin, B., Moore, C., Stefanakis, T.S., Zhou, H., Spillane, M. & Synolakis, C.E., 2013. Focusing of long waves with finite crest over constant depth, *Proc. R. Soc. A*, **469**, doi:10.1098/rspa.2013.0015.

Lay, T., Ammon, C.J., Kanamori, H., Koper, K.D., Sufri, O. & Hutko, A.R., 2010. Teleseismic inversion for rupture process of the 27 February 2010 Chile (Mw 8.8) earthquake, *Geophys. Res. Lett.*, **37**, L13301, doi:10.1029/2010GL043379.

Liu, P.L.-F., Cho, Y.-S., Yoon, S.B. & Seo, S.N., 1994. Numerical simulations of the 1960 Chilean tsunami propagation and inundation at Hilo, Hawaii. Recent Developments in Tsunami Research, pp. 99-115.

Mosocoso, E., Grevemeyer, I., Contreras, E., Flueh, E., Dzierma, Y., Rabbel, W. & Thorwart, M., 2011. Revealing the deep structure and rupture plane of the 2010 Maule, Chile earthquake ($M_w$ = 8.8) using wide angle seismic data, *Earth planet. Sci. Lett.*, **307**, 147-155.

Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space, *Bull. seism. Soc. Am.*, **75**, 1135-1154.

Okal, E.A. & Synolakis, C.E., 2004. Source discrimination for nearfield tsunamis, *Geophys. J. Int.*, **158**, 899-912.

Okal, E.A., Borrero, J.C. & Synolakis, C.E., 2006. Evaluation of Tsunami Risk from Regional Earthquakes at Pisco, Peru, *Bull. seism. Soc. Am.*, **96**, 1634-1648.

Russel, S., 1844. Report on waves, in *Proceedings of the 14th Meeting of the British Association for the Advancement of Science*, John Murray, pp. 311.

Shuto, N., 1991. Numerical simulation of tsunami—its present and near future, *Nat. Hazards*, **4**, 171-191.

Synolakis, C.E., 1987. The runup of solitary waves, *J. Fluid Mech.*, **185**, 523-545.

Synolakis, C.E., 1988. On the roots of $J_0(z) - iJ_1(z)$, *Q. Appl. Maths., XLVI*(1), 105-108.

Synolakis, C.E., 1991. Tsunami run-up on steep slopes: how good linear theory really is, *Nat. Hazards*, **4**, 221-234.

Tadepalli, S. & Synolakis, C.E., 1994. The run-up on N-waves on sloping beaches, *Proc. R. Soc. Lond. A*, **455**, 99-112.

Tadepalli, S. & Synolakis, C.E., 1996. Model for the leading waves of tsunamis, *Phys. Rev. Lett.*, **77**, 2141-2145.

Ward, S., 2010. *Tsunami. Encyclopedia of Solid Earth Geophysics*, Springer Press.

Yamazaki, Y. & Cheung, K.F., 2011. Shelf resonance and impact of near-field tsunami generated by the 2010 Chile earthquake, *Geophys. Res. Lett.*, **38**, L12605, doi:10.1029/2011GL047508.
Yamazaki, Y., Kowalik, Z. & Fai Cheung, K., 2009. Depth-integrated, non-hydrostatic model for wave breaking and run-up, *Int. J. Numer. Methods Fluids*, **61**(5), 473–497.

Zhang, J.E. & Wu, Y.T., 1999. Oblique long waves on beach and induced longshore current, *J. Eng. Mech.*, **125**(7), 812–826.

## APPENDIX A: GENERAL TREATMENT FOR 1+1 LINEAR THEORY

We use through out the manuscript the following Fourier transform definition. If $f(\vec{x}, t) : \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}$, we define the direct and inverse $(n+1)$-dimensional time-space Fourier transform respectively as,

$$ \bar{\bar{f}}(\vec{k}, \omega) = \int_{-\infty}^{\infty} \int_{\mathbb{R}^n} f(\vec{x}, t) e^{i\vec{k}\cdot\vec{x} - i\omega t} d\vec{x} dt $$

$$ f(\vec{x}, t) = \frac{1}{(2\pi)^{n+1}} \int_{-\infty}^{\infty} \int_{\mathbb{R}^n} \bar{\bar{f}}(\vec{k}, \omega) e^{-i\vec{k}\cdot\vec{x} + i\omega t} d\vec{k} d\omega, $$

where $\bar{\bar{}}$ denotes the time-space Fourier transform, $\vec{x} \in \mathbb{R}^n$ is the location vector and $\vec{k} \in \mathbb{R}^n$ is the wavenumber vector. We adopt the standard notation $\hat{\{\}}$ when only one variable is transformed.

Starting with the linearized PDE (7), there are two cases, depending on the domain we seek the solution:

(i) $x > x_0$.

At the open sea, the bathymetry is constant: $h(x) = d$, then eq. (7) is reduced to,

$$ \eta_{tt} - g d\eta_{xx} = 0 \qquad (A1) $$

which is the classical homogeneous wave equation with velocity propagation given by $c_0^2 = g d$. A general solution in the frequency domain for eq. (A1) is,

$$ \hat{\eta}(x, \omega) = A_i(\omega)e^{-ik_0(\omega)x} + A_r(\omega)e^{ik_0(\omega)x}, \qquad (A2) $$

where $A_i$ and $A_r$ represent the incident and the reflected coefficients, respectively. Applying the inverse Fourier transform for the water elevation, the solution is

$$ \eta(x, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \{ A_i(\omega)e^{i\omega t - ik_0(\omega)x} + A_r(\omega)e^{i\omega t + ik_0(\omega)x} \} d\omega. \quad (A3) $$

(ii) $0 < x < x_0$.

This domain corresponds to the sloping beach region, where the bathymetry is $h(x) = x \tan(\beta)$, so eq. (7) reduces to

$$ \eta_{tt} - g \tan(\beta)(x\eta_{xx} + \eta_x) = 0. \qquad (A4) $$

Defining a change of variables (from $x$ to a $v$) as follow,

$$ v := \frac{2}{c_0} \sqrt{x_0 x} $$

$$ \eta(x, t) := \zeta(v, t). $$

PDE (A4) can be transformed into,

$$ v\zeta_{vv} + \zeta_v - v\zeta_{tt} = 0. \qquad (A5) $$

Applying the time-space Fourier transform, eq. (A5) becomes an ODE,

$$ \frac{\partial}{\partial k} \left[ (k^2 - \omega^2) \bar{\bar{\zeta}} \right] = k \bar{\bar{\zeta}}. \qquad (A6) $$

Solving the previous ODE, one obtains

$$ \bar{\bar{\zeta}} = \frac{C}{\sqrt{w^2 - k^2}} H(|w| - |k|) - i \frac{C}{\sqrt{k^2 - w^2}} H(|k| - |w|) \quad (A7) $$

with $C = C(\omega)$ a function independent of $k$ and $H(\cdot)$ is the Heaviside step function.

Applying the inverse time-space Fourier transform, it is possible to write a general solution,

$$ \zeta(v, t) = \frac{1}{(2\pi)^2} \int_{-\infty}^{\infty} C e^{i\omega t} \\ \times \left( \int_{-|\omega|}^{|\omega|} \frac{e^{-ikv} dk}{\sqrt{w^2 - k^2}} - i \int_{|k|>|\omega|} \frac{e^{-ikv} dk}{\sqrt{k^2 - w^2}} \right) d\omega. \quad (A8) $$

In virtue of integral forms of the cylindrical Bessel functions (Abramowitz & Stegun 1964):

$$ J_0(\xi) = \frac{2}{\pi} \int_0^1 \frac{\cos(u\xi)}{\sqrt{1 - u^2}} du $$

$$ Y_0(\xi) = \frac{2}{\pi} \int_1^\infty \frac{\cos(u\xi)}{\sqrt{u^2 - 1}} du $$

and just taking the regular solutions at $x = 0$, the general solution in this case, over the sloping beach domain, can be written finally in terms of the original variables as,

$$ \eta(x, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} B(\omega)e^{i\omega t} J_0(2k_0\sqrt{xx_0}) d\omega. \quad (A9) $$

In the frequency domain, the solution is simply

$$ \hat{\eta}(x, \omega) = B(\omega) J_0(2k_0 \sqrt{xx_0}). \qquad (A10) $$

## A1 Matching continuity conditions of the solutions

The next step is to set some conditions in order to relate the $B(\omega)$ transmitted coefficient with $A_i$ and $A_r$.

We match eqs (A2) and (A10) by imposing continuity of $\hat{\eta}$ at $x = x_0$ and also, continuity on their first derivatives. Then,

$$ A_i(\omega)e^{ik_0 x_0} + A_r(\omega)e^{-ik_0 x_0} = B(\omega)J_0(2k_0 x_0) $$

$$ A_i(\omega)e^{ik_0 x_0} - A_r(\omega)e^{-ik_0 x_0} = -iB(\omega)J_1(2k_0 x_0) $$

and retrieving $B(\omega)$ and $A_r(\omega)$ in terms of the incident coefficient, $A_i(\omega)$, we obtain

$$ B(\omega) = \frac{2A_i(\omega)e^{ik_0 x_0}}{J_0(2kx_0) + iJ_1(2k_0x_0)} $$

$$ A_r(\omega) = A_i(\omega)e^{2i\left[k_0 x_0 - \arctan\left(\frac{J_1(2k_0 x_0)}{J_0(2k_0 x_0)}\right)\right]} $$

If we give the boundary value at $x = x_0$, that is to say, $\hat{\eta}(x_0, \omega) = \Phi(\omega)$, this means that $A_i(\omega) = e^{-ik_0 x_0} \Phi(\omega)$.

Finally, the transmitted solution to the sloping beach domain is,

$$ \eta(x, t) = \frac{1}{\pi} \int_{-\infty}^{+\infty} \Phi(\omega) \frac{J_0(2k_0 \sqrt{x_0 x}) e^{i\omega t}}{J_0(2k_0 x_0) + iJ_1(2k_0 x_0)} d\omega, \quad (A11) $$

which is the solution proposed by Synolakis (1987).

## APPENDIX B: GENERAL TREATMENT FOR 2+1 LINEAR THEORY

Similar to the solving approach shown in Appendix A, we start with PDE (18) and the bathymetry described by eq. (17), there are two cases:

(i) $x > x_0$
At the open sea, we obtain again the 2+1 wave equation with
$c_0^2 = gd$

$$
\eta_{xx} + \eta_{yy} = \frac{1}{c_0^2} \eta_{tt} \qquad (B1)
$$

then, the general solution can be written as

$$
\eta(x, y, t) = \frac{1}{(2\pi)^2} \iint_{\mathbb{R}^2} (A_i e^{ik_x x} + A_r e^{-ik_x x}) e^{i\omega t - ik_y y} dk_y d\omega, \quad (B2)
$$

where the relation $k_x^2 + k_y^2 = k_0^2$ holds with $k_0 = \frac{\omega}{c_0}$.

(ii) $0 < x < x_0$.

Applying the time-space Fourier transform to eq. (18) over the
$(y, t)$ variables, the PDE transforms into an ODE,

$$
x \bar{\eta}_{xx} + \bar{\eta}_x + (k_0^2 x_0 - k_y^2 x) \bar{\eta} = 0 \quad (B3)
$$

which is a hypergeometric differential equation.

With standard techniques for solving ODEs, the general solution
for (B3) is

$$
\bar{\eta} = A e^{-k_y x} U(\mu, 1, 2k_y x) + B e^{-k_y x} M(\mu, 1, 2k_y x), \quad (B4)
$$

where $\mu = \frac{1}{2}(1 - \frac{k_y^2 x_0}{k_y})$, and $M, U$ are the confluent hypergeometric functions of first and second kind, respectively. $U$ is divergent at $x=0$, so we set $A=0$.

Applying the same matching continuity conditions, as in Section
A.1, one obtains the following coefficients,

$$
\begin{align*}
B(k_y, \omega) &= \frac{2A_i(k_y, \omega)e^{ik_x x_0}}{X(x_0) - \frac{i}{k_x}X'(x_0)} \\
A_r(k_y, \omega) &= \frac{X(x_0) + \frac{i}{k_x}X'(x_0)}{X(x_0) - \frac{i}{k_x}X'(x_0)} A_i(k_y, \omega)e^{2ik_x x_0},
\end{align*}
$$

where

$$
X(x) = e^{-k_y x} M(\mu, 1, 2k_y x). \tag{B5}
$$

The solution in the sloping beach domain can be written as

$$
\bar{\eta}(x, k_y, \omega) = \frac{2A_i(k_y, \omega)e^{ik_x x_0}}{X(x_0) - \frac{i}{k_x} \frac{\partial X}{\partial x}} X(x).
$$

Now, applying the inverse Fourier transform to the previous expres-
sion, the solution is

$$
\eta(x, y, t) = \frac{1}{(2\pi)^2} \iint_{\mathbb{R}^2} \frac{2\Phi(k_y, \omega)X(x)e^{i\omega t - ik_y y}}{X(x_0) - \frac{i}{k_x} \frac{\partial X}{\partial x}(x_0)} dk_y d\omega. \quad (B6)
$$

B1 Incident solitary wave and limit case

Considering the incoming solitary wave (22), its Fourier spectrum
is,

$$
\Phi(k_y, \omega) = \frac{4d^3}{3} \frac{\pi}{c} e^{ik_c \sin(\theta)(x-x_1)} k_c \cosech(\alpha k_c) \cdot 2\pi \delta[k_y - \cos(\theta)k_c], \quad (B7)
$$

where $\delta(\cdot)$ is the Dirac delta distribution.

Inserting (B7) in (B6) and using the fact $\delta$ is the Dirac delta (that
works like an evaluation functional), we get

$$
\eta(x, y, t) = \frac{4d^3}{3} \int_{-\infty}^{\infty} \frac{X(x)k_c \cosech(\alpha k_c)}{X(x_0) - \frac{i}{ak_c}X'(x_0)} \times e^{ik_c[\sin(\theta)(x_0-x_1) + \cos(\theta)y + ct]} dk_c. \quad (B8)
$$

Here, X has been evaluated with $k_y = \cos(\theta)k_c$. So when $\theta \to \frac{\pi}{2}$, it
can be proven, using the series expansion of the Kummer function
(Abramowitz & Stegun 1964), that

$$
\lim_{\theta \to \frac{\pi}{2}} X(x) = J_0(2ak_c \sqrt{x}x_0),
$$

and this limit was also used by Zhang & Wu (1999). This allows us
to write that

$$
\lim_{\theta \to \frac{\pi}{2}} \eta(x, y, t) = \frac{1}{\pi} \int_{-\infty}^{+\infty} \Phi(\omega) \frac{J_0(2k_0 \sqrt{x_0 x}) e^{i\omega t}}{J_0(2k_0 x_0) + i J_1(2k_0 x_0)} d\omega. \quad (B9)
$$

Which is exactly the same solution obtained in (A11). Evidently,
this will produce the same 1+1 runup (16) independent of y.