Tsunami propagation kernel and its applications

Tsunamis rarely occur in a specific area, and their occurrence is highly uncertain. Suddenly generated from their sources in deep water, they occasionally undergo tremendous amplification in shallow water to devastate low-lying coastal areas. Despite the advancement of computational power and simulation algorithms, there is a need for novel and rigorous approaches to efficiently predict coastal amplification of tsunamis during different disaster management phases, such as tsunami risk assessment and real-time forecast. This study presents convolution kernels that can instantly predict onshore waveforms of water surface elevation and flow velocity from observed/simulated wave data away from the shore. Kernel convolution involves isolating an incident-wave component from the offshore wave data and transforming it into the onshore waveform. Moreover, unlike previously derived ones, the present kernels are based on shallow-water equations with a damping term and can account for tsunami attenuation on its path to the shore with a damping parameter. Kernel convolution can be implemented at a low computational cost compared to conventional numerical models that discretise the spatial domain. The prediction capability of the kernel method was demonstrated through application to real-world tsunami cases.

1 Introduction

Tsunamis pose a major threat to low-lying coastal areas worldwide. They are initiated by the rapid displacement of seawater, which is often triggered by earthquakes and/or submarine landslides. Given the initial water surface displacement, we are currently able to simulate tsunami propagation from the source to coastal areas and then assess tsunami hazards in coastal communities with a practical level of accuracy. However, the occurrence of a tsunami is highly uncertain; thus, we need to prepare for potential tsunami hazards by considering many different source scenarios that are often based on scarce historical data. This involves performing numerous tsunami simulations that are often implemented by a two-step approach, i.e. deep-water tsunami modelling and coastal tsunami modelling. The first step predicts tsunami propagation from the source to coastal shelves; the linear shallow-water equations are solved on a relatively coarse grid. Then, the second step simulates coastal tsunami evolution using non-linear shallow-water equations discretised on a finer grid that can resolve coastal bathymetry. Although the two steps can be combined by grid refinement algorithms, the coastal tsunami simulation still requires high computational costs; this may limit the range of uncertainty covered by the scenarios. There is a demand for novel approaches to rapidly predict coastal tsunami evolution despite growing computational power. Such methods will also contribute to the real-time prediction of coastal tsunami impacts from a tsunami waveform that is observed or predicted in deep water.

The coastal evolution of tsunamis has been analytically studied by many researchers using shallow-water equations because the vertical motion of seawater is small compared to the horizontal motion during the propagation phase. The main physics of wave transformation can be described by linear shallow-water equations, while the non-linear effect is responsible for the wave distortion in shallow water (Carrier and Greenspan, 1958). When the wave height is large, relative to water depth, short-wavelength components exhibit wave breaking owing to non-linear wave distortion, which considerably dissipates the wave energy. However, most past tsunamis in world oceans are known to have propagated as non-breaking waves. Analytical solutions for non-breaking wave evolution over a uniform slope have been derived for incident transient waves in different forms, which allow us to predict tsunami run-up height on the shore by prototype incident waves (Synolakis, 1987; Tadepalli and Synolakis, 1994; Pelinovsky and Mazova, 1992). This approach can be further extended for non-uniform seabed profiles (Didenkulova et al., 2009) and narrow bays (Zahibo et al., 2006; Rybkin et al., 2014; Shimozono, 2016). The resulting run-up formulas have provided valuable insights into the coastal evolution processes of tsunamis in different forms. However, actual tsunami waveforms are diverse and differ from these prototype waves. Moreover, the highest impact is not necessarily produced by a leading part of tsunamis, especially in cases of far-field sources. Therefore, the parametric formulas for specific wave types do not work directly for real-world tsunami problems.

A possible way to rapidly predict a real-world tsunami transformation is to use kernel representation of long wave propagation over a sloping coast. A coastal tsunami profile on a uniform slope can be expressed in the frequency domain as a product of a Bessel transfer function and an integral transform of the incident-wave profile (Synolakis, 1991). This, in turn, means that the tsunami waveform in the time domain results from a convolution of the incident-wave profile with a kernel function, which is the inverse transform of the Bessel transfer function. Kernel convolution enables us to predict a tsunami waveform given an incident-wave profile of any integral form. A bottleneck for this approach was that the Bessel transfer function could not be readily inverted into the time-domain kernel function. Madsen and Schäffer (2010) derived an approximate kernel for a wave profile on the shore using asymptotic expressions of Bessel functions. More recently, the author of this paper derived a general kernel for a wave profile at an arbitrary location on the slope, which yields an exact solution of shallow-water equations via the convolution with the incident-wave time series (Shimozono, 2020). There have been few attempts to apply kernel representation for the rapid prediction of tsunamis (Choi et al., 2011; Chan and Liu, 2012) because it still has some limitations for real-world applications. First, the real-world problem cannot be represented as the incident-wave boundary value problem. The observed or simulated tsunamis in deep water more or less contain reflected waves from the coast. In addition, real-world tsunamis are attenuated by different factors as they propagate to the shore. The kernel representation based on frictionless shallow-water equations does not account for the wave decay that becomes significant in shallow water.

As an extension of the author's earlier work (Shimozono, 2020), this paper presents a new kernel representation that has higher applicability to real-world tsunami problems. Unlike the previous ones derived under the incident-wave boundary condition, the presented kernel functions are formulated with the observed-wave boundary condition that may include a contribution of reflected waves from the coast. Moreover, the kernel is derived from shallow-water equations with a damping term. Therefore, a damping factor is incorporated into the kernel, which can be used to represent the tsunami decay on its path to the shore. The practical kernel works when an onshore tsunami profile needs to be predicted directly from observed or simulated wave data away from the shore. The next section presents the derivation of kernel functions for water surface elevation and flow velocity from the shallow-water equations under the observed-wave boundary condition. Subsequently, Sect. 3 presents an efficient numerical method of implementing kernel convolution. Section 4 demonstrates the predictive capability of the kernel method through its applications to 2011 Tohoku tsunami cases. Finally, Sect. 5 provides brief conclusions of this study.

2 Kernel functions

2.1 Kernel representation

Let us consider one-dimensional propagation of a tsunami from the deep ocean to coastal water, which is initially at the stationary state, as illustrated in Fig. 1. We define the x^* axis as the positive seaward horizontal coordinate with an origin at the still-water shoreline. There is a wave observation point at $x^{*}=l$, where water surface displacements driven by the incoming tsunami are continuously measured. The seabed slope landward of the observation point is assumed to be uniform, while the seabed on its sea side may have an arbitrary form as long as it smoothly connects to the landward uniform slope. The measured wave profile contains both landward-propagating waves from the source and seaward-propagating waves reflected from the shore. Our problem is to determine a coastal wave profile under the observed-wave boundary condition. This is different from the previous formalism in which coastal wave evolution is predicted under the incident-wave boundary condition by assuming a horizontal bed seaward of the boundary (Keller and Keller, 1964; Carrier, 1971; Synolakis, 1991).

Figure 1 Schematic illustration of the boundary value problem.

Here, we use linear shallow-water equations to formulate the problem. Wave non-linearity becomes significant in shallow water, and the non-linear effect distorts the time and spatial axes. However, it is well known that the wave amplitude is not significantly affected by non-linearity unless the non-linear wave distortion leads to wave breaking (Carrier and Greenspan, 1958; Tuck and Hwang, 1972; Synolakis, 1991). The non-linear shoreline motion can be readily derived from the linear solution via the hodograph transform, and the run-up height is unchanged from the linear case if the boundary value assignment is linearised (e.g. Pelinovsky and Mazova, 1992). Furthermore, the occurrence of wave breaking, which limits the applicable range of the present approach, can be predicted as a breakdown point of the hodograph transform under the same condition. While the linearised boundary value assignment potentially affects the run-up height and the wave breaking condition, non-linear modifications are minor as long as the ratio of wave amplitude to water depth is small at the boundary (Antuono and Brocchini, 2007, 2008). Therefore, the main process of practical interest can be described by linear equations when we place the boundary in deep water. Another important fact, which has often been neglected in analytical models, is that a tsunami decays on its path to the shore by different factors. Therefore, we incorporate a damping term into the momentum equation to account for the tsunami decay. Accordingly, the damped propagation of tsunamis over a uniform slope can be described by the following equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\eta }}^{*}}{\partial t^{*}}}+\mathit{\beta }{u}^{*}+\mathit{\beta }{x}^{*}{\displaystyle \frac{\partial u^{*}}{\partial x^{*}}}=\mathrm{0},\text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial u^{*}}{\partial t^{*}}}+g{\displaystyle \frac{\partial {\mathit{\eta }}^{*}}{\partial x^{*}}}=d^{*}.\end{array}$$

Here, the asterisk superscript represents a dimensional variable. Therefore, t^* represents time, ${\mathit{\eta }}^{*}(x^{*},t^{*})$ is the water surface elevation above the still-water level, $u^{*}(x^{*},t^{*})$ is the positive seaward horizontal velocity, β is the seabed slope, d^* is the damping term and g is the gravitational acceleration.

The quadratic law is often used to model tsunami damping ($d^{*}\propto \left|u^{*}\right|u^{*}$), but this introduces non-linearity into the problem, which makes the analytical work infeasible. Instead, we employ the linear damping term as follows:

$$\begin{array}{}\text{(3)}& d^{*}=-{\mathit{\alpha }}^{*}{u}^{*},\end{array}$$

where α^* (≥0) is the damping parameter, which has a dimension of the reciprocal of time. The linear damping term introduces an exponential wave decay with the characteristic time corresponding to ${\mathit{\alpha }}^{*-\mathrm{1}}$ into the non-dissipative solution (Davies et al., 2020). In virtue of its simple formulation, the linear model has been used to represent damped propagation of tsunamis in deep water. Many researchers have suggested the damping parameter of 𝒪(10⁻⁵) s⁻¹ for deep-water tsunami propagation through comparisons of the model and observations (e.g. Fine et al., 2013; Kulikov et al., 2014). However, the damping parameter could be much higher for nearshore tsunami dissipation. Mazova et al. (1990) analytically derived a dissipative run-up formula of monochromatic long waves over a uniform slope using the linear damping term. They suggested the α^* value of 𝒪(10⁻²) s⁻¹ in a typical situation of calculating tsunami run-up.

Tsunamis decay in shallow water by different factors, such as skin and form drag over the seabed as well as local topography and coastal structures, the scales of which are much smaller than the tsunami wavelength. Here, we attempt to represent the tsunami decay due to different factors collectively by the linear damping term; thus, the parameter α^* works as an empirical parameter rather than as a physical parameter. This treatment makes the resulting kernel representation a semi-empirical one. We will discuss the variability in this parameter through real-world tsunami applications in a later section.

The governing equations are non-dimensionalised by the following dimensionless variables:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{x}^{*}=xl,t^{*}=t\sqrt{l/g\mathit{\beta }},{\mathit{\eta }}^{*}=\mathit{\eta }\mathit{\beta }l,u^{*}=u\sqrt{g\mathit{\beta }l},\\ \text{(4)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}^{*}=\mathit{\alpha }\sqrt{g\mathit{\beta }/l}.\end{array}$$

Here, the spatial scales and timescales are chosen to be the horizontal slope length l and half of the wave travel time over the entire slope $\sqrt{l/g\mathit{\beta }}$, respectively. Namely, the horizontal dimensionless coordinate takes the value of x=0 at the shoreline and x=1 at the offshore boundary. The dimensionless time of t=2 equals the one-way travelling time from/to the boundary to/from the shore. Of note is that the shallow-water approximation requires $l\ll L_{\mathrm{0}}/\mathit{\beta }$, where L₀ is the representative wavelength.

Substituting Eq. (4) into Eqs. (1) and (2) yields the dimensionless forms of the governing equations:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+u+x{\displaystyle \frac{\partial u}{\partial x}}=\mathrm{0},\text{(6)}& {\displaystyle \frac{\partial u}{\partial t}}+{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}+\mathit{\alpha }u=\mathrm{0}.\end{array}$$

Eliminating η from Eqs. (5) and (6) yields a single-variable wave equation of u as

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial t^{\mathrm{2}}}}+\mathit{\alpha }{\displaystyle \frac{\partial u}{\partial t}}-\mathrm{2}{\displaystyle \frac{\partial u}{\partial x}}-x{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}=\mathrm{0}.\end{array}$$

We solve this equation for wave evolution over initially stationary water on a uniform slope given an offshore wave profile of arbitrary form at x=1. For this purpose, we apply Laplace transform to Eq. (7). Then, we have

$$\begin{array}{}\text{(8)}& (s^{\mathrm{2}}+\mathit{\alpha }s)\widehat{u}-\mathrm{2}{\displaystyle \frac{\partial \widehat{u}}{\partial x}}-x{\displaystyle \frac{{\partial }^{\mathrm{2}}\widehat{u}}{\partial x^{\mathrm{2}}}}=\mathrm{0},\end{array}$$

where $\widehat{u}$ is the Laplace transform of u, which is given by

$$\begin{array}{}\text{(9)}& \widehat{u}(x,s)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}u(x,t)e^{-st}\mathrm{d}t\end{array}$$

with a complex number frequency parameter s. Equation (8) is the modified Bessel's equation and has a bounded solution at the shoreline (x=0) as

$$\begin{array}{}\text{(10)}& \widehat{u}(x,s)=G\left(s\right)x^{-\frac{\mathrm{1}}{\mathrm{2}}}{I}_{\mathrm{1}}\left(\mathrm{2}\sqrt{s\left(s+\mathit{\alpha }\right)}{x}^{\frac{\mathrm{1}}{\mathrm{2}}}\right),\end{array}$$

where I_n is the modified Bessel function of the first kind and the nth order and G(s) is an arbitrary function that will be determined by the offshore boundary condition.

Accordingly, the water surface elevation in the Laplace domain, $\widehat{\mathit{\eta }}(x,s)$, can be obtained from Eqs. (5) and (10) such that

$$\begin{array}{ll}{\displaystyle }\widehat{\mathit{\eta }}(x,s)& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{s}}\left(\widehat{u}+x{\displaystyle \frac{\partial \widehat{u}}{\partial x}}\right)\\ \text{(11)}& {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\sqrt{s\left(s+\mathit{\alpha }\right)}}{s}}G\left(s\right)I_{\mathrm{0}}\left(\mathrm{2}\sqrt{s\left(s+\mathit{\alpha }\right)}{x}^{\frac{\mathrm{1}}{\mathrm{2}}}\right).\end{array}$$

To determine G(s), we use the boundary condition, $\widehat{\mathit{\eta }}(\mathrm{1},s)={\widehat{\mathit{\eta }}}_{\mathrm{0}}\left(s\right)$, where ${\widehat{\mathit{\eta }}}_{\mathrm{0}}\left(s\right)$ is the Laplace transform of the observed time series of water surface elevation at x=1. Then, we can express the solutions for water surface elevation and horizontal velocity as

$$\begin{array}{}\text{(12)}& {\displaystyle }\widehat{\mathit{\eta }}(x,s)=s{\widehat{\mathrm{\Psi }}}_{\mathrm{0}}(x,s){\widehat{\mathit{\eta }}}_{\mathrm{0}}\left(s\right),\text{(13)}& {\displaystyle }\widehat{u}(x,s)=-s{\widehat{\mathrm{\Psi }}}_{\mathrm{1}}(x,s){\widehat{\mathit{\eta }}}_{\mathrm{0}}\left(s\right),\end{array}$$

where ${\widehat{\mathrm{\Psi }}}_{\mathrm{0}}(x,s)$ and ${\widehat{\mathrm{\Psi }}}_{\mathrm{1}}(x,s)$ are the transfer functions for $\widehat{\mathit{\eta }}$ and $\widehat{u}$ in the Laplace domain, respectively, which are given as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\widehat{\mathrm{\Psi }}}_{\mathrm{0}}(x,s)={\displaystyle \frac{I_{\mathrm{0}}\left(\mathrm{2}\sqrt{s\left(s+\mathit{\alpha }\right)}{x}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)}{s{I}_{\mathrm{0}}\left(\mathrm{2}\sqrt{s\left(s+\mathit{\alpha }\right)}\right)}},\\ \text{(14)}& {\displaystyle }& {\displaystyle }{\widehat{\mathrm{\Psi }}}_{\mathrm{1}}(x,s)={\displaystyle \frac{x^{-\frac{\mathrm{1}}{\mathrm{2}}}}{\sqrt{s\left(s+\mathit{\alpha }\right)}}}{\displaystyle \frac{I_{\mathrm{1}}\left(\mathrm{2}\sqrt{s\left(s+\mathit{\alpha }\right)}{x}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)}{I_{\mathrm{0}}\left(\mathrm{2}\sqrt{s\left(s+\mathit{\alpha }\right)}\right)}}.\end{array}$$

Here, $\widehat{\mathit{\eta }}$ and $\widehat{u}$ are expressed as the product of each transfer function and $s{\widehat{\mathit{\eta }}}_{\mathrm{0}}$ in the Laplace domain. The factor of s was isolated in Eqs. (12) and (13) such that the transfer functions are invertible to the time domain. Consequently, η and u in the time domain are commonly expressed by the convolution of the rate of water surface displacement at the offshore boundary with the inverse Laplace transforms of the transfer functions as follows:

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }\mathit{\eta }(x,t)={\mathrm{\Psi }}_{\mathrm{0}}(x,t)*{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{d}t}},\text{(16)}& {\displaystyle }& {\displaystyle }u(x,t)=-{\mathrm{\Psi }}_{\mathrm{1}}(x,t)*{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{d}t}},\end{array}$$

where the asterisk denotes the convolution operation. Ψ_n(x,t) results from the inverse Laplace transform of ${\widehat{\mathrm{\Psi }}}_n(x,s)$, which is expressed using the Bromwich integral such that

$$\begin{array}{}\text{(17)}& {\mathrm{\Psi }}_n(x,t)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }i}}\underset{\mathit{\gamma }-i\mathrm{\infty }}{\overset{\mathit{\gamma }+i\mathrm{\infty }}{\int }}{\widehat{\mathrm{\Psi }}}_n\left(s\right)e^{st}\mathrm{d}s.\end{array}$$

Here, γ is a positive real number greater than every singularity of ${\widehat{\mathrm{\Psi }}}_n\left(s\right)$. The common form of solutions in Eqs. (15) and (16) suggests that the convolution kernels accommodate the general processes of wave transformation under the observed-wave boundary condition. An exception arises at the shoreline (x=0), where ${\widehat{\mathrm{\Psi }}}_{\mathrm{1}}$ is not well defined. The two transfer functions in Eq. (14) can be related at the shoreline as ${\widehat{\mathrm{\Psi }}}_{\mathrm{1}}(\mathrm{0},s)=s{\widehat{\mathrm{\Psi }}}_{\mathrm{0}}(\mathrm{0},s)$, which turns into their relation in the time domain as ${\mathrm{\Psi }}_{\mathrm{1}}(\mathrm{0},t)=\partial {\mathrm{\Psi }}_{\mathrm{0}}(\mathrm{0},t)/\partial t$. With this relation, the kinematic condition at the shoreline results from Eqs. (15) and (16):

$$\begin{array}{}\text{(18)}& u_s=-{\displaystyle \frac{\partial \mathit{\xi }}{\partial t}},\end{array}$$

where $u_s\left(t\right)\equiv u(\mathrm{0},t)$ represents the shoreline velocity and $\mathit{\xi }\left(t\right)\equiv \mathit{\eta }(\mathrm{0},t)$ expresses the vertical shoreline displacement. This suggests that the choice of the bounded solution at the shoreline in Eq. (10) is equivalent to imposing the linearised kinematic condition. It should be emphasised that the present solutions are based on the moving boundary condition at the shoreline. The shoreline velocity can be obtained from Eqs. (15) and (18).

2.2 Kernel derivation

The next step is to derive the kernel functions by evaluating the Bromwich integrals in Eq. (17). The residue theorem can be used to determine the inverse Laplace transform of the transfer function. Both ${\widehat{\mathrm{\Psi }}}_{\mathrm{0}}(x,s)$ and ${\widehat{\mathrm{\Psi }}}_{\mathrm{1}}(x,s)$ have an infinite number of poles at s=s_k such that $I_{\mathrm{0}}\left(\mathrm{2}\sqrt{s_k(s_k+\mathit{\alpha })}\right)=\mathrm{0}$. Additionally, ${\widehat{\mathrm{\Psi }}}_{\mathrm{0}}(x,s)$ has a pole at s=0, while the singularity of ${\widehat{\mathrm{\Psi }}}_{\mathrm{1}}(x,s)$ at s=0 is removable. The complex zeros of the modified Bessel function are distributed on the imaginary axis; thus, they can be associated with positive real zeros of the Bessel function of the first kind and zeroth order, c_k, that satisfies J₀(c_k)=0. All poles of the two functions can be summarised as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{s}_{\mathrm{0}}=\mathrm{0},s_k={\displaystyle \frac{-\mathit{\alpha }+i\sqrt{c_k^{\mathrm{2}}-{\mathit{\alpha }}^{\mathrm{2}}}}{\mathrm{2}}}\text{and}\\ \text{(19)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{s}}_k={\displaystyle \frac{-\mathit{\alpha }-i\sqrt{c_k^{\mathrm{2}}-{\mathit{\alpha }}^{\mathrm{2}}}}{\mathrm{2}}},\end{array}$$

where k is a natural number ($k=\mathrm{1},\mathrm{2},\mathrm{3}\mathrm{\dots }$) and c_k is the positive zeros of J₀ in ascending order of magnitude, as graphed in Fig. 2. Equation (19) suggests that the damping factor has an upper limit such that $\mathit{\alpha }<c_{\mathrm{1}}\approx \mathrm{2.405}$. When k is large, the positive zeros of the Bessel function can be approximated as

$$\begin{array}{}\text{(20)}& c_k\approx \left(k-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right)\mathit{\pi }.\end{array}$$

This approximation suggests that c_k linearly increases with k, as shown in Fig. 2.

Figure 2 Positive zeros of the Bessel function of the first kind and zeroth order; J₀(c_k)=0.

Because both functions are holomorphic elsewhere, the Bromwich integrals can be evaluated using the theory of residue such that

$$\begin{array}{ll}{\displaystyle }{\mathrm{\Psi }}_{\mathrm{0}}(x,t)& {\displaystyle }=\text{Res}\left[{\widehat{\mathrm{\Psi }}}_{\mathrm{0}}\right(s_{\mathrm{0}}\left)e^{s_{\mathrm{0}}t}\right]\\ {\displaystyle }& {\displaystyle }+\sum _{k=\mathrm{1}}^{\mathrm{\infty }}\left\{\text{Res}\left[{\widehat{\mathrm{\Psi }}}_{\mathrm{0}}\right(s_k\left)e^{s_{k}t}\right]+\text{Res}\left[{\widehat{\mathrm{\Psi }}}_{\mathrm{0}}\right({\overline{s}}_k\left)e^{{\overline{s}}_{k}t}\right]\right\},\\ \text{(21)}& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{1}}(x,t)& {\displaystyle }=\sum _{k=\mathrm{1}}^{\mathrm{\infty }}\left\{\text{Res}\left[{\widehat{\mathrm{\Psi }}}_{\mathrm{1}}\right(s_k\left)e^{s_{k}t}\right]+\text{Res}\left[{\widehat{\mathrm{\Psi }}}_{\mathrm{1}}\right({\stackrel{\mathrm{‾}}{s}}_k\left)e^{{\stackrel{\mathrm{‾}}{s}}_{k}t}\right]\right\},\end{array}$$

where $\text{Res}\left[\widehat{\mathrm{\Psi }}\right(s\left)e^{st}\right]$ is the residue of $\widehat{\mathrm{\Psi }}\left(s\right)e^{st}$ at the point s. Substituting Eqs. (14) and (19) into Eq. (21) yields

$$\begin{array}{}\text{(22)}& {\mathrm{\Psi }}_{\mathrm{0}}(x,t)=\mathrm{1}+{\mathit{\psi }}_{\mathrm{0}}(x,t),{\mathrm{\Psi }}_{\mathrm{1}}(x,t)=x^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\psi }}_{\mathrm{1}}(x,t)\end{array}$$

with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_n(x,t)=-\mathrm{2}{e}^{-\frac{\mathit{\alpha }t}{\mathrm{2}}}\sum _{k=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{J_n\left(c_k{x}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)}{{\mathit{\lambda }}_k{J}_{\mathrm{1}}\left(c_k\right)}}\sin \left({\displaystyle \frac{{\mathit{\lambda }}_k}{\mathrm{2}}}t+{\displaystyle \frac{n\mathit{\pi }}{\mathrm{2}}}+{\mathit{\theta }}_k\right),\\ \text{(23)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_k=\sqrt{c_k^{\mathrm{2}}-{\mathit{\alpha }}^{\mathrm{2}}},{\mathit{\theta }}_k=\text{atan}\mathrm{2}\left({\mathit{\lambda }}_k,\mathit{\alpha }\right).\end{array}$$

Both kernels consist of infinite series of the combination of the Bessel function of x^½ and the sinusoidal function of t. In addition, the kernels contain the negative exponential function of time and decay with time at a higher rate for larger α. For α=0 without the damping effect, the kernels exhibit oscillatory behaviour in the infinite time length; thus, the kernel convolution involves a full wave history from the initial time at the offshore boundary. Introduction of the damping effect confines the causal relation to the near past. Kernel convolution deals with the separation of incident and reflected waves because it is formulated with the observed-wave profile containing two-way components. Separation requires a longer observed-wave history for a smaller damping factor.

Figure 3 Function Ψ₀(x,t) at (a) x=0, (b) $x=\mathrm{1}/\mathrm{8}$ and (c) $x=\mathrm{1}/\mathrm{2}$. Each panel shows graphs of the function for α=0, α=0.2 and α=1.0. Vertical dotted lines represent the timing of singularities, $t_m^{\pm }$.

Figure 4 Function Ψ₁(x,t) at (a) $x=\mathrm{1}/\mathrm{8}$, (b) $x=\mathrm{1}/\mathrm{2}$ and (c) x=1. Each panel shows the graphs of the function for α=0, α=0.2 and α=1.0. Vertical dotted lines represent the timing of singularities, $t_m^{\pm }$.

Figures 3 and 4 show graphs of Ψ₀(x,t) and Ψ₁(x,t), respectively, on the time axis for different values of the damping factor; namely, α=0, α=0.2 and α=1.0. Ψ₀(x,t) is shown at x=0, $x=\mathrm{1}/\mathrm{8}$ and $x=\mathrm{1}/\mathrm{2}$, while Ψ₁(x,t) is graphed at $x=\mathrm{1}/\mathrm{8}$, $x=\mathrm{1}/\mathrm{2}$ and x=1. Of note is that ${\mathrm{\Psi }}_{\mathrm{0}}(\mathrm{1},t)=\mathrm{1}$ and Ψ₁(0,t) are not well defined. The initial forms of the two kernels agree with those derived under the incident-wave boundary condition (Shimozono, 2020) because the initial part of observed-wave data represents a shoreward-propagating wave. They deviate from the incident-wave kernels in the large time domain exhibiting cyclic singularities. Cyclic singularities occur at the following timing:

$$\begin{array}{}\text{(24)}& t_m^{\pm }=\mathrm{2}\left\{\mathrm{1}\pm (-\mathrm{1}{)}^m{x}^{\frac{\mathrm{1}}{\mathrm{2}}}\right\}+\mathrm{4}(m-\mathrm{1}),\end{array}$$

where m is the natural number. Both Ψ₀ and Ψ₁ diverge to positive or negative infinity at $t_m^{-}$ and form slope discontinuities at $t_m^{+}$. Of note is that 2(1−x^½) corresponds to the propagation time from the offshore boundary to the reference point, whereas 2(1+x^½) represents the propagation time between the two locations via reflection on the shore. Hence, the two kernels have cyclic singularities representing arrivals of past seaward- and shoreward-propagating waves with a period of t=4, which equals the wave round-trip time over the entire slope. Equation (24) can be alternatively described as a series of characteristic curves in the (x,t) plane, $x={\left\{t-\mathrm{2}-\mathrm{4}(m-\mathrm{1})\right\}}^{\mathrm{2}}/\mathrm{4}$, which represent forth and back propagation of a wave signal on the slope. Kernel convolution generates a coastal waveform by processing the two-way propagating waves at the offshore boundary. At the shoreline (x=0) where $t_m^{+}=t_m^{-}$, the divergence and slope discontinuity merge with each other and the bimodal kernel structure disappears from Ψ₀. A similar merging of singularities also happens to Ψ₁ at the offshore boundary (x=1). The infinite series in Eq. (23) is convergent except at $t_m^{-}$, where it turns into a harmonic series for large k; the convergence properties of the infinite series are more rigorously discussed in Appendix A. Therefore, the infinite series can be numerically calculated until the point of convergence except at the singular points.

Substituting Eq. (22) into Eqs. (15) and (16) yields

$$\begin{array}{}\text{(25)}& {\displaystyle }\mathit{\eta }(x,t)={\mathit{\eta }}_{\mathrm{0}}(t-t_{\mathrm{0}})+\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{\mathit{\psi }}_{\mathrm{0}}(x,\mathit{\tau }){\left.{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{d}t}}\right|}_{t-\mathit{\tau }}\mathrm{d}\mathit{\tau },\text{(26)}& {\displaystyle }u(x,t)=-x^{-\frac{\mathrm{1}}{\mathrm{2}}}\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{\mathit{\psi }}_{\mathrm{1}}(x,\mathit{\tau }){\left.{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{d}t}}\right|}_{t-\mathit{\tau }}\mathrm{d}\mathit{\tau },\end{array}$$

where $t>t_{\mathrm{0}}\equiv \mathrm{2}(\mathrm{1}-x^{\frac{\mathrm{1}}{\mathrm{2}}})$ and $\mathit{\eta }=u=\mathrm{0}$ otherwise. These representations suggest that $\mathit{\eta }(x,t)\approx {\mathit{\eta }}_{\mathrm{0}}(t-t_{\mathrm{0}})$ and $u(x,t)\approx \mathrm{0}$ hold when the incoming wave is much longer than the slope length. This can be viewed as the seabed slope working like a vertical wall against a relatively long wave. The convolution describes wave deformation over the slope under the observed-wave boundary condition. Despite cyclic singularities, the kernel convolution is well defined, and the convolution can be numerically performed, which will be presented in the next section.

3 Kernel convolution

Kernel representations derived in the previous sections allow us to predict coastal waveforms from observed-wave data away from the shore. The prediction can be made instantly by kernel convolution of the offshore wave data. This is advantageous over conventional ways of numerically solving the shallow-water equations, which involve both time and spatial integration of the equations over a physical domain. Particularly, an accurate simulation of shoreline response requires high-resolution discretisation. Additionally, the numerical model approaches require iterative calculations for the separation of the incident and reflected waves. Despite the efficiency benefits, the numerical convolution requires special treatment for cyclic kernel singularities. Here, an efficient algorithm for the singular kernel convolution is presented and then validated in a simple case of monochromatic-wave propagation.

3.1 Numerical method

The kernels have cyclic singularities at $t_m^{+}$ and $t_m^{-}$. For the convenience of presenting the convolution method, the singular time points are redefined as follows:

$$\begin{array}{}\text{(27)}& t_{jm}=\mathrm{2}\left\{\mathrm{1}-(-\mathrm{1}{)}^j{x}^{\frac{\mathrm{1}}{\mathrm{2}}}\right\}+\mathrm{4}(m-\mathrm{1}),\end{array}$$

where j is 0 or 1 ($t_m^{-}=t_{\mathrm{0}m}$ and $t_m^{+}=t_{\mathrm{1}m}$). To avoid midpoint singularities, the integration interval is divided into multiple time segments: $[t_{\mathrm{01}},t_{\mathrm{11}}],[t_{\mathrm{11}},t_{\mathrm{02}}],[t_{\mathrm{02}},t_{\mathrm{12}}],\mathrm{\dots },[t_{\mathrm{0}m},t_{\mathrm{1}m}],[t_{\mathrm{1}m},t_{\mathrm{0}m+\mathrm{1}}],\mathrm{\dots }$ Then, based on the kernel form, different variable changes are applied in [t_0m,t_1m] and $[t_{\mathrm{1}m},t_{\mathrm{0}m+\mathrm{1}}]$ based on the double exponential integration formula by Takahasi and Mori (1974) via

$$\begin{array}{ll}{\displaystyle }{\mathit{\tau }}_{\mathrm{0}m}\left(p_{\mathrm{0}}\right)=& {\displaystyle }\mathrm{2}\left\{\mathrm{2}m-\mathrm{1}+x^{\frac{\mathrm{1}}{\mathrm{2}}}\tanh \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\sinh p_{\mathrm{0}}\right)\right\}\\ \text{(28)}& {\displaystyle }& {\displaystyle }\text{for}[t_{\mathrm{0}m},t_{\mathrm{1}m}],{\displaystyle }{\mathit{\tau }}_{\mathrm{1}m}\left(p_{\mathrm{1}}\right)=& {\displaystyle }\mathrm{2}\left\{\mathrm{2}m+(\mathrm{1}-x^{\frac{\mathrm{1}}{\mathrm{2}}})\tanh \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\sinh p_{\mathrm{1}}\right)\right\}\\ \text{(29)}& {\displaystyle }& {\displaystyle }\text{for}[t_{\mathrm{1}m},t_{\mathrm{0}m+\mathrm{1}}].\end{array}$$

These changes of variables will concentrate integration points of time near the singular points using the double exponential function. Consequently, each segment is mapped to an infinite interval of p_j with vanishing endpoint singularities. For shoreline motions (x=0), the segments of j=0 disappear because t_0m=t_1m. Because the integrand decays with a double exponential rate, actual numerical integration is performed within a finite interval of p_j. Therefore, the convolutions in Eqs. (25) and (26) are rewritten as

$$\begin{array}{}\text{(30)}& {\displaystyle }\begin{array}{rl}\mathit{\eta }& (x,t)={\mathit{\eta }}_{\mathrm{0}}(t-t_{\mathrm{0}})+\sum _{m=\mathrm{1}}^M\sum _{j=\mathrm{0}}^{\mathrm{1}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathit{\psi }}_{\mathrm{0}}(x,{\mathit{\tau }}_{jm}(p_j\left)\right)\\ & {\left.{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{d}t}}\right|}_{t-{\mathit{\tau }}_{jm}\left(p_j\right)}{\displaystyle \frac{\mathrm{d}{\mathit{\tau }}_{jm}}{\mathrm{d}{p}_j}}\mathrm{d}{p}_j,\end{array}\text{(31)}& {\displaystyle }\begin{array}{rl}u& (x,t)=-x^{-\frac{\mathrm{1}}{\mathrm{2}}}\sum _{m=\mathrm{1}}^M\sum _{j=\mathrm{0}}^{\mathrm{1}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathit{\psi }}_{\mathrm{1}}(x,{\mathit{\tau }}_{jm}(p_j\left)\right)\\ & {\left.{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{d}t}}\right|}_{t-{\mathit{\tau }}_{jm}\left(p_j\right)}{\displaystyle \frac{\mathrm{d}{\mathit{\tau }}_{jm}}{\mathrm{d}{p}_j}}\mathrm{d}{p}_j,\end{array}\end{array}$$

where M is an integer ceiling of $(t-t_{\mathrm{0}})/\mathrm{4}$ when α=0 and can be a smaller number when α>0 because the kernels decay on the time axis. The convolution is numerically performed by applying the segment-by-segment integration following the trapezoidal rule. The convolution is efficiently implemented by the advance tabulation of ${\mathit{\psi }}_n\mathrm{d}{\mathit{\tau }}_{jm}/\mathrm{d}{p}_j$. We read the kernel table and then perform computation of the infinite integrals in several segments to predict a coastal waveform.

3.2 Validation

The convolution method is demonstrated and validated by comparing numerical results with an analytical solution for a simple problem. Here, we consider a simple problem in which a monochromatic wave is observed away from the shore that contains both incident and reflected waves. The observed wave is now given as

$$\begin{array}{}\text{(32)}& {\mathit{\eta }}_{\mathrm{0}}\left(t\right)=\sin \left(\mathit{\lambda }t\right),\end{array}$$

where λ is the dimensionless angular frequency. It is assumed that the coastal wave is in an equilibrium state, and a partial standing wave is formed owing to wave decay over the uniform slope.

To derive the equilibrium solution, we first take the Laplace transform of Eq. (32) and substitute it into Eq. (12). Then, we have

$$\begin{array}{}\text{(33)}& \widehat{\mathit{\eta }}(x,s)={\displaystyle \frac{\mathit{\lambda }{I}_{\mathrm{0}}\left(\mathrm{2}\sqrt{s(s+\mathit{\alpha })}{x}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)}{I_{\mathrm{0}}\left(\mathrm{2}\sqrt{s(s+\mathit{\alpha })}\right)(s+i\mathit{\lambda })(s-i\mathit{\lambda })}}.\end{array}$$

This can be inversely transformed into the time domain by the residue theorem. However, the residues associated with the zeros of the modified Bessel function, s_k and ${\stackrel{\mathrm{‾}}{s}}_k$ in Eq. (19), yield a transient solution because they accommodate the decaying exponential function of time with the damping factor α. Therefore, we need to consider only the residues of $s=\pm i\mathit{\lambda }$ for the equilibrium solution in the case of α>0.

Consequently, the damped equilibrium solution for water surface elevation, $\stackrel{\mathrm{̃}}{\mathit{\eta }}(x,t)$, is given by

$$\begin{array}{}\text{(34)}& \stackrel{\mathrm{̃}}{\mathit{\eta }}(x,t)=a\left(x\right)\sin \left\{\mathit{\lambda }t+\mathit{\varphi }\left(x\right)\right\},\end{array}$$

where the wave amplitude a(x) and the phase ϕ(x) are given by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }a\left(x\right)={\displaystyle \frac{\sqrt{{\left\{AC\left(x\right)+BD\left(x\right)\right\}}^{\mathrm{2}}+{\left\{AD\left(x\right)-BC\left(x\right)\right\}}^{\mathrm{2}}}}{A^{\mathrm{2}}+B^{\mathrm{2}}}},\\ \text{(35)}& {\displaystyle }& {\displaystyle }\mathit{\varphi }\left(x\right)=\text{atan}\mathrm{2}\left\{AD\left(x\right)-BC\left(x\right),AC\left(x\right)+BD\left(x\right)\right\}\end{array}$$

with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }A=\text{Re}\left(I_{\mathrm{0}}(\mathrm{2}\sqrt{i\mathit{\lambda }(i\mathit{\lambda }+\mathit{\alpha })}\right),\\ {\displaystyle }& {\displaystyle }B=\text{Im}\left(I_{\mathrm{0}}(\mathrm{2}\sqrt{i\mathit{\lambda }(i\mathit{\lambda }+\mathit{\alpha })}\right),\\ {\displaystyle }& {\displaystyle }C\left(x\right)=\text{Re}\left(I_{\mathrm{0}}(\mathrm{2}\sqrt{i\mathit{\lambda }(i\mathit{\lambda }+\mathit{\alpha })x}\right),\\ {\displaystyle }& {\displaystyle }D\left(x\right)=\text{Im}\left(I_{\mathrm{0}}(\mathrm{2}\sqrt{i\mathit{\lambda }(i\mathit{\lambda }+\mathit{\alpha })x}\right).\end{array}$$

The wave amplitude on the slope cannot be represented by simple mathematical functions and exhibits drastic changes with varying λ because the observation point moves between the node and anti-node of the partial standing wave.

To validate the numerical convolution method, the wave amplitude is computed for non-zero α values using Eq. (30). Because the damping effect confines the causal relation in a finite time, we set the total time length of the convolution such that the kernel value sufficiently decays ($t=\mathrm{12}/\mathit{\alpha }$). The numerical convolution was performed over a sufficiently large interval ($-\mathrm{3}<p_j<\mathrm{3}$) for both j=0 and 1 with a common subinterval Δp_j. To investigate the sensitivity of the result to Δp_j, the numerical solutions were obtained with three different subintervals, i.e. Δp_j=0.05, 0.1, and 0.25.

Figure 5 Wave amplitude a as a function of dimensionless angular frequency λ at different locations (x=0, $x=\mathrm{1}/\mathrm{8}$ and $x=\mathrm{1}/\mathrm{2}$) for damping coefficients (a) α=0.1 and (b) α=1.0. Each panel shows exact and numerical solutions with different integration subintervals Δp_j=0.05, Δp_j=0.25 and Δp_j=0.1.

Figure 5 shows the comparison of wave amplitude from numerical convolution and the exact solution at three different locations (x=0, $\mathrm{1}/\mathrm{2}$ and $\mathrm{1}/\mathrm{8}$). The two cases with different values of the damping parameter, α=0.1 and α=1.0, are shown in the left and right columns, respectively. In each figure, exact and numerical solutions of wave amplitude are compared on the λ axis. The dimensionless angular frequency can be expressed as $\mathit{\lambda }=\mathrm{2}\mathit{\pi }l/L_{\mathrm{0}}$ with the representative wavelength L₀; thus, λ=2π represents the situation in which the wavelength is equal to the slope length. In the case of α=0.1, wave amplitude shows large fluctuations on the λ axis owing to considerable reflection. The coastal wave amplitude becomes very large when a partial node is formed at x=1 where wave amplitude is fixed to be unity. The fluctuations are smeared out in the case of α=1.0 owing to the stronger damping effect. Numerical results agree well with the exact values in this range of λ when Δp_j<0.1. Numerical error appears from the high-frequency side when Δp_j=0.25. These results confirm that numerical convolution produces an exact solution with a sufficiently small subinterval relative to the tsunami wave period.

Figure 6 Kernel validation for the formation of a full standing wave ($\mathit{\lambda }=c_k/\mathrm{2}$): (a) incident wave, η_i(t); (b) observed wave at x=1, η₀(t); and (c) shoreline elevation data computed from η_i(t) and η₀(t).

In the case of α=0, the monochromatic wave creates a complete node at x=1 when $\mathit{\lambda }=c_k/\mathrm{2}$. Therefore, we have no wave signal at x=1 in the equilibrium state, and the wave amplitude given by Eq. (35) goes to infinity because the boundary value given by Eq. (32) is not appropriate. The kernel method works even in such a case as long as the full standing wave is formed from an initially stationary state. To demonstrate this, we consider a transient incident wave given by

$$\begin{array}{}\text{(36)}& {\mathit{\eta }}_i\left(t\right)=\tanh \left({\displaystyle \frac{\mathit{\lambda }t}{\mathrm{10}}}\right)\sin \left(\mathit{\lambda }t\right),\end{array}$$

where the tanh function was introduced for a smooth transition from the stationary water surface to a monochromatic wave. As illustrated in Fig. 6a, the incident wave grows to a monochromatic wave over several wave periods. When we set $\mathit{\lambda }=c_{\mathrm{1}}/\mathrm{2}\approx \mathrm{1.2024}$, a node will be formed at x=1 under the incident wave. Using the incident-wave kernel previously proposed by the author (Shimozono, 2020), we can obtain the water surface elevation at x=1, η₀, as shown in Fig. 6b. The null elevation datum occurs after the initial transient phase due to the formation of a full standing wave. Figure 6c compares shoreline elevation profiles, ξ, computed from η_i using the incident-wave kernel and that computed from η₀ using the present kernel. The two results show a perfect agreement, and the wave amplitude converges to the analytical solution of the monochromatic-wave run-up height by Keller and Keller (1964): $\mathrm{2}/\sqrt{J_{\mathrm{0}}(c_{\mathrm{1}}{)}^{\mathrm{2}}+J_{\mathrm{1}}(c_{\mathrm{1}}{)}^{\mathrm{2}}}\approx \mathrm{3.85}$. Even if the null datum occurs at x=1 as a result of the superposition of incident and reflected waves, the kernel convolution can predict a waveform at any location on the slope from the initial transient part of the wave history at x=1. It is worth emphasising that the present kernel works for such a problem because it is constructed for the initial-boundary value problem. Without the initial condition, we could not derive an incident-wave signal from offshore wave data when a full node is formed at the boundary.

4 Kernel applications

We can instantly predict a coastal waveform over the slope from the offshore wave profile via the kernel convolution. To demonstrate its capability in real-world tsunami problems, the kernel method is applied to cases of the 2011 Tohoku tsunami, for which instrumentally recorded wave data are available at different locations along the affected coastline (Kawai et al., 2013). Based on the availability of multiple wave data in the same area, two cases were selected on the northern side of the epicentre as shown in Fig. 7a. Case 1 is located far north of the tsunami source, but tsunami run-up on the coast still reached several metres above the mean sea level. The offshore wave data are available relatively close to the coast, and the coastal tide station also recorded the water surface fluctuation due to the tsunami. In contrast, Case 2 is located in one of the most devastated parts of the coastline, where three offshore stations were aligned perpendicularly to the coast. Located close to the tsunami source, the initial wave was dominant over successive waves, and the tsunami run-up reached up to 20 m in this area (Shimozono et al., 2012). Despite the lack of a coastal wave record, this case is employed owing to the availability of multiple observations in the cross-shore direction. The leading part of the tsunami at each observed location is assumed to have been refracted in the deep water and approached the coastline nearly perpendicularly. This assumption may not be valid for the most offshore point of Case 2 (TM1), where water depth exceeds 1500 m and the depth contour is significantly curved. To apply the kernel method, the water depth is assumed to linearly decrease over the distance between the observed location and the shore in each case, even though the actual seabed slope is not strictly uniform. The cross-shore bathymetry and locations of wave observation in each case are shown in Fig. 7b and c.

Figure 7 Map of the case study sites and wave observation stations for the 2011 Tohoku tsunami: (a) overview of the two case study sites at the northeastern Pacific coast of Japan and (b) and (c) seabed profiles along the cross-shore transects for Case 1 and Case 2, respectively. The markers indicate the locations of wave observation points.

4.1 Case 1 – Mutsu-Ogawara

Case 1 is located at the coast of Mutsu-Ogawara, where relatively small waves of similar amplitude repeatedly attack the coast in a period of approximately 30 min. The offshore wave data were recorded by an ultrasonic wave gauge installed on the seabed at a depth of 43.8 m, approximately 3.4 km off the coast (P1 in Fig. 7b). In addition, an onshore wave record is available from a tide station inside Mutsu-Ogawara Port (P2 in Fig. 7b). Because the offshore wave data contained high-frequency fluctuations that may not be real waves propagating shoreward, the following analyses were performed with the low-pass-filtered time series (<0.003 Hz) as shown in Fig. 8a. The numerical convolution was performed with the filtered wave data with three different values of the damping parameter, α=0.1, α=0.5 and α=1.0, on the assumption of a uniform slope ($\mathit{\beta }=\mathrm{43.8}/\mathrm{3400}\approx \mathrm{0.013}$). The numerical integration in each segment was implemented for $-\mathrm{3}<p_j<\mathrm{3}$ with Δp_j=0.05. The numerical error was confirmed to be negligibly small.

Figure 8 Prediction of the shoreline motions in Case 1: (a) observed and low-pass-filtered wave data at P1 and (b–d) waveforms at the shoreline computed with α=0.1, α=0.5 and α=1.0, respectively, in comparison to the coastal wave record at P2.

Figure 8b and c compare the computed waveforms at the shoreline, ξ^*, and the observed-wave data at P2 with the different values of damping parameters, respectively. The first term in Eq. (25), ${\mathit{\eta }}_{\mathrm{0}}^{*}(t-t_{\mathrm{0}})$, is also plotted to highlight the slope effect described by kernel convolution. The comparison of ξ^* and ${\mathit{\eta }}_{\mathrm{0}}^{*}$ shows that relatively short wave components are significantly amplified over the slope. The computed results with α=0.1 overestimate the short waves on the shore. The observed-wave data at P2 are in good agreement with the computed result of α=0.5, but the initial waveforms are reproduced slightly better with α=1.0. This result suggests that the damping factor varies with wave amplitude, and the non-linear formulation of the damping term may be required to achieve higher accuracy. Nonetheless, despite the crude assumptions, the kernel method works well to predict the coastal waveform from the observed data a few kilometres off the coast when the α value is appropriately chosen.

4.2 Case 2 – Kamaishi and Ryoishi

Case 2 is located on the coast ranging from Kamaishi to Ryoishi, areas which were highly devastated during the 2011 tsunami event. Two bottom pressure sensors (TM1 and TM2) captured the initial waveform of the deep-water tsunami propagating towards the coast (Maeda et al., 2011; ERI, 2011). As shown in Fig. 7, TM1 was located 76 km off the coast at a depth of approximately 1600 m, while TM2 was 47 km off the coast at a depth of approximately 1000 m. In addition, a GPS buoy (GPS802) deployed 15 km off the coast recorded a full profile of the tsunami propagating at a depth of approximately 200 m (Kawai et al., 2013). Because three observation points are aligned on a cross-shore line, the dataset can be used to validate the kernel method for the long-distance propagation of the real-world tsunami with a large amplitude. Coastal wave gauges in shallow water were destroyed by the tsunami; thus, there are no data for onshore waveforms. Nevertheless, intensive post-event surveys provide the range of coastal run-up heights in this area (Shimozono et al., 2012). The ria coast exhibits an intricate coastline, but the dimension of coastline variations is smaller than the tsunami wavelength. Therefore, tsunami propagation is expected to be described with the one-dimensional kernel approach to a large extent.

Figure 9 Prediction of coastal waveforms in Case 2: (a) observed-wave data at TM1, (b) predicted waveforms based on TM1 and observed wave at TM2, (c) predicted waveforms based on TM1 and observed wave at GPS802, (d) predicted shoreline motion based on TM1, (e) observed-wave data at TM2, (f) predicted wave based on TM2 and observed wave at GPS802, and (g) predicted shoreline motion based on TM2.

Figure 9a shows the observed waveform at TM1. The leading part of the tsunami in this region was characterised by a short, impulsive wave of large amplitude riding on a relatively long wave. We predict the waveforms at the locations of TM2 and GPS803 via the kernel convolution of the wave data at TM1, assuming the water depth to be linearly decreasing from TM1 to the shore (β=0.02). Figure 9b and c compare the computed waveform for α=0.1, α=0.5 and α=1.0 with the observed data at TM2 and GPS802, respectively. The computed profile agrees with the observed one at both TM1 and GPS802 when α=0.5. However, the water surface elevation is predicted to be higher after the peak at both locations, and this discrepancy probably occurs because the observed wave at TM1 is located far from the coast and contains wave components that do not propagate shoreward. The resulting run-up height on the shore is sensitive to the choice of the damping parameter, as shown in Fig. 9d. The measured run-up heights after the event significantly varied owing to the intricate coastline; however, the maximum run-up height was up to 20 m in the coastal areas. Therefore, the α value of 0.5–1.0 should be employed for the purpose of predicting the run-up of the tsunami in this area.

Next, we predict the waveform at the location of GPS803 using the wave data at TM2 graphed in Fig. 9e. The computed results at GPS802 show better agreement with the observation than those predicted from TM1, as shown in Fig. 9f. The results in Fig. 9c and f confirm that the observed wave at TM1 did not fully propagate shoreward and suggest that TM1 is too far from the shore for the kernel application. The lower value of α realises better agreement at GPS802, but the maximum run-up height is overestimated with the same value (Fig. 9c). This poses the limitation of the linear formulation of the damping effect that cannot account for relatively high tsunami attenuation in shallow water owing to the quadratic dependence on flow velocity. In particular, the local topographic elements along the intricate coastline caused additional damping effects in the coastal area. Therefore, the empirical damping parameter should be optimised for the target location of the prediction. Figure 10 shows another prediction of the wave profile at the shoreline from the wave data at GPS802 that were available for a much longer period than at TM2. The shoreline displacement shows a similar form to the previous predictions, and the low damping parameter of α=0.1 produces successive high waves that are not observed in this area. The onshore waveform of the leading wave in the case of α=0.5–1.0 agrees well with the previous model result based on the two-dimensional non-linear shallow-water equations (Shimozono et al., 2012).

Figure 10 Prediction of shoreline motions in Case 2: (a) observed-wave data at GP802 and (b) predicted shoreline motion based on GPS802.

The kernel representation of coastal tsunami evolution provides an efficient method for predicting onshore tsunami waveforms based on a single wave observation away from the shore. Coastal waveforms are almost simultaneously obtained with the observed profile because the computation time for the numerical integration is negligibly small compared to the tsunami timescale. In this case, the lead times for the real-time prediction of shoreline motions are 16 and 8 min when predicted from TM2 and GPS802, respectively. The prediction accuracy is limited by the crude assumptions underlying kernel representation that makes the fast prediction possible: one-dimensional propagation over a uniform slope with the linear damping effect. Nonetheless, the damping parameter can be optimised for the prediction of tsunamis at a specific coastal location from an offshore observation point; this can be done with the help of a numerical model. The two cases suggest that coastal run-up is well predicted when the damping parameter is set in the range of 0.5–1.0.

5 Conclusions

This study presents the tsunami propagation kernel that compactly accommodates damped propagation processes of tsunamis over the coastal slope. The kernel representation was derived under the observed-wave boundary condition unlike the previous kernels based on an incident-wave boundary condition. Therefore, it can be directly applied to predict a coastal waveform from observed-wave data that contain both shoreward- and seaward-propagating components. Furthermore, the damping factor was incorporated into the kernel to collectively represent tsunami attenuation that occurs in various ways. Both the water surface elevation and horizontal flow velocity over the slope can be represented as the kernel convolution of the rate of water surface displacement away from the shore. Kernel functions have cyclic singularities to separate the two-way propagating wave in the offshore wave history. The introduction of the damping factor confines this causal relation to the near past. The kernel convolution can be efficiently performed via the change of variables using the double exponential function, and the numerical convolution method worked satisfactorily if the integration subinterval was chosen to be sufficiently small relative to the tsunami wave period.

Kernel representation has potential applicability to real-world problems because it can instantly predict the response of coastal water to incoming wave trains in deep water. It can be applied for the real-time prediction of coastal waveforms from a single offshore observation by bottom-pressure-type or buoy-type wavemeters that are widely deployed off the coasts under an imminent threat of tsunamis. In addition, it can be incorporated into deep-water tsunami simulations to efficiently predict onshore waveforms without resolving coastal bathymetry. Prediction accuracy depends on the choice of the damping parameter because the linear damping term cannot fully represent actual tsunami decay that occurs in diverse ways. Therefore, the damping parameter should be treated as an empirical parameter and can be optimised at each target site through pre-calibration based on numerical simulations. Even though the number of cases is limited, the damping parameter of 0.5–1.0 works for different cases to reasonably reproduce onshore waveforms. Additional studies will be needed to confirm the applicability of the semi-empirical kernel and to improve kernel formulation for better representation of tsunami decay in shallow water.

Appendix A: Kernel convergence properties

This appendix presents the convergence properties of the infinite series in Eq. (23). Firstly, we define the sequence of the series as

$$\begin{array}{}\text{(A1)}& f_{nk}(x,t)={\displaystyle \frac{J_n\left(c_k{x}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)}{{\mathit{\lambda }}_k{J}_{\mathrm{1}}\left(c_k\right)}}\sin \left({\displaystyle \frac{{\mathit{\lambda }}_k}{\mathrm{2}}}t+{\displaystyle \frac{n\mathit{\pi }}{\mathrm{2}}}+{\mathit{\theta }}_k\right).\end{array}$$

Then, the infinite series can be rewritten as

$$\begin{array}{}\text{(A2)}& {\mathit{\psi }}_n(x,t)=-\mathrm{2}{e}^{-\frac{\mathit{\alpha }t}{\mathrm{2}}}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{f}_{\text{nk}}(x,t).\end{array}$$

We look into the behaviour of the sequence for large k. When k is large, we can use the approximation of c_k in Eq. (20) and the following approximations as well:

$$\begin{array}{}\text{(A3)}& {\mathit{\lambda }}_k\approx c_k\approx \left(k-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right)\mathit{\pi },{\mathit{\theta }}_k\approx {\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}.\end{array}$$

Furthermore, the Bessel functions can be approximated based on their asymptotic expansion for a large augment:

$$\begin{array}{}\text{(A4)}& J_n\left(z\right)\approx \sqrt{{\displaystyle \frac{\mathrm{2}}{\mathit{\pi }z}}}\cos \left(z-{\displaystyle \frac{n\mathit{\pi }}{\mathrm{2}}}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\right),\end{array}$$

which is valid only for large z. Since the augment of the Bessel function in the numerator of Eq. (A1) contains x^½, we discuss the property of the sequence separately in two cases: (a) x=0 and (b) x>0.

A1 Case a, x=0

Here we discuss the behaviour of only f_0k(0,t) because f_1k(0,t) is not well defined. Since J_n(0)=1, we can approximate Eq. (A1) for large k as

$$\begin{array}{}\text{(A5)}& f_{\mathrm{0}k}(\mathrm{0},t)\approx {\displaystyle \frac{(-\mathrm{1}{)}^{k-\mathrm{1}}}{\sqrt{\mathrm{2}k-\frac{\mathrm{1}}{\mathrm{2}}}}}\cos \left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(k-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right)\mathit{\pi }t\right\}.\end{array}$$

The sequence decays with increasing k while oscillating due to the combination of the alternating sign and cosine function of k. Therefore, the infinite series of the sequence is generally convergent. However, the oscillations disappear when $t=t_m^{\pm }=\mathrm{4}m-\mathrm{2}$ such that

$$\begin{array}{}\text{(A6)}& f_{\mathrm{0}k}(\mathrm{0},t_m^{\pm })\approx {\displaystyle \frac{-\mathrm{1}}{\sqrt{\mathrm{2}k-\frac{\mathrm{1}}{\mathrm{2}}}}}\cos \left\{{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}-{\displaystyle \frac{m\mathit{\pi }}{\mathrm{2}}}\right\}.\end{array}$$

This sequence forms a general harmonic series, and thus, the kernel slowly diverges to positive or negative infinity depending on the m value at $t=t_m^{\pm }$ as shown in Fig. 3a.

A2 Case b, x>0

When x>0, the sequence can be approximated for large k as

$$\begin{array}{ll}{\displaystyle }{f}_{\text{nk}}(x,t)\approx & {\displaystyle }{\displaystyle \frac{(-\mathrm{1}{)}^k}{\mathit{\pi }{x}^{\frac{\mathrm{1}}{\mathrm{4}}}\left(k-\frac{\mathrm{1}}{\mathrm{4}}\right)}}\cos \left\{\mathit{\pi }\left(k-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right)x^{\frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{n\mathit{\pi }}{\mathrm{2}}}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\right\}\\ \text{(A7)}& {\displaystyle }& {\displaystyle }\cdot \cos \left\{{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\left(k-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right)t+{\displaystyle \frac{n\mathit{\pi }}{\mathrm{2}}}\right\}.\end{array}$$

This sequence also exhibits both oscillation and decay with increasing k, and thus, the infinite series is convergent except at singular points where the oscillations cease. Substituting $t=t_m^{\pm }$ into (A7) and separating it into two parts, we have

$$\begin{array}{ll}{\displaystyle }{f}_{nk}& {\displaystyle }(x,t_m^{\pm })\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{x}^{\frac{\mathrm{1}}{\mathrm{4}}}\left(k-\frac{\mathrm{1}}{\mathrm{4}}\right)}}\\ {\displaystyle }& {\displaystyle }\cdot \cos \left[\mathit{\pi }\left(k-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right)x^{\frac{\mathrm{1}}{\mathrm{2}}}\left\{\mathrm{1}\pm (-\mathrm{1}{)}^m\right\}-{\displaystyle \frac{m\mathit{\pi }}{\mathrm{2}}}\right]\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{(-\mathrm{1}{)}^n}{\mathrm{2}\mathit{\pi }{x}^{\frac{\mathrm{1}}{\mathrm{4}}}\left(k-\frac{\mathrm{1}}{\mathrm{4}}\right)}}\\ \text{(A8)}& {\displaystyle }& {\displaystyle }\cdot \sin \left[\mathit{\pi }\left(k-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right)x^{\frac{\mathrm{1}}{\mathrm{2}}}\left\{\mathrm{1}\mp (-\mathrm{1}{)}^m\right\}+{\displaystyle \frac{m\mathit{\pi }}{\mathrm{2}}}\right].\end{array}$$

The original sequence is now expressed as a sum of the two sequences. When $t=t_m^{+}$, one of the sequences vanishes and the other one forms a convergent series. Therefore, the infinite series is convergent at $t=t_m^{+}$. On the other hand, one sequence forms a convergent series and the other one generates a harmonic series when $t=t_m^{-}$. This confirms that the infinite series in Eq. (A2) slowly diverges to positive and negative infinity alternatively at $t_m^{-}$ as shown in Figs. 3b and c and 4.