Tsunami propagation kernel and its applications

Tsunamis rarely occur in a specific area, and their occurrence is highly uncertain. Generated from their sources in deep water, they occasionally undergo tremendous amplification over decreasing water depth : 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 5 onshore waveforms of water surface elevation and flow velocity from observed/simulated wavedata apart from the shore. Kernel convolution involves isolating an incident-wave component from the offshore wavedata 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. 10 The prediction capability of the kernel method was demonstrated through application to real-world tsunami cases.

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 propaga-40 tion over a sloping coast. 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 45 readily inverted into the time-domain kernel function. Madsen and Schäffer (2010) derived an approximate kernel for wave profile on the shore using asymptotic expressions of Bessel functions. More recently, the author derived the general kernel for 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 50 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 55 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 Figure 1. Schematic illustration of the boundary value problem.
wavedata apart 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, Section 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, Section 5 provides brief conclusions of this study.

Kernel represenation
Let us consider one-dimensional propagation of a tsunami from 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 70 its seaside 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).

80
Therefore, the main physics :::::: process :: of :::::::: practical :::::: interest : can be described by the linear equations as long as the ratio of wave amplitude to water depth is small at the boundary (Antuono and Brocchini, 2007). This condition is not restrictive when we place the boundary in deep water :::::::::::::::::::::::::: (Antuono and Brocchini, 2007). Another important fact, which has been often 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 85 uniform slope can be described by the following equations: Here, the asterisk superscript represents a dimensional variable. Therefore, t * represents time, η * (x * , t * ) is the water surface 90 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 * ∝ |u * | u * ), but this introduces non-linearity into the problem, which makes the analytical work infeasible. Instead, we employ the linear damping term as follows;

110
Here, the spatial and time scales are chosen to be the horizontal slope length, l and half of the wave travel time over the entire slope l/gβ, 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 to the one-way travelling time from/to the boundary to/from the shore. Of note, the shallow water approximation requires l L 0 /β, where L 0 is the representative wavelength.
Eliminating η from (5) and (6) yields a single-variable wave equation of u as 120 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 (7). Then, we have whereû is the Laplace transform of u, which is given bŷ 125 with a complex number frequency parameter s. Equation (8) 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.

130
Accordingly, the water surface elevation in the Laplace domain,η(x, s), can be obtained from (5) and (10) such that To determine G(s), we use the boundary condition,η(1, s) =η 0 (s), whereη 0 (s) 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 u(x, s) = −sΨ 1 (x, s)η 0 (s), whereΨ 0 (x, s) andΨ 1 (x, s) are the transfer functions forη andû in the Laplace domain, respectively, which are given aŝ 140 Here,η andû are expressed as the product of each transfer function and sη 0 in the Laplace domain. The factor of s was isolated in (12) and (13) where the asterisk denotes the convolution operation. Ψ n (x, t) results from the inverse Laplace transform ofΨ n (x, s), which is expressed using the Bromwich integral such that Here, γ is a positive real number greater than every singularity ofΨ n (s). The common form of solutions in (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 Ψ 1 ::: Ψ 1 is not well defined. The two transfer functions in (14) can be related at the shoreline asΨ 1 (0, s) = sΨ 0 (0, s), which turns into their relation in time domain as Ψ 1 (0, s) = ∂Ψ 0 (0, s)/∂t :::::::::: Ψ 1 (0, t) = With this relation, the kinematic condition at the shoreline results from (15) and (16); where u s (t) ≡ u(0, t) represents the shoreline velocity, and ξ(t) ≡ η(0, t) expresses the vertical shoreline displacement. This suggests that the choice of the bounded solution at the shoreline in (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 (15) and (18).
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 Ψ 0 . A similar merge of singularities also happens to Ψ 1 at the offshore boundary (x = 1). The infinite series in (23) is convergent except at t − m , where it turns into a harmonic series for large k; the convergence properties of the 205 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 (22) into (15) and (16) yields where t > t 0 ≡ 2(1 − x 1 2 ) and η = u = 0 otherwise. These representations suggest that η(x, t) ≈ η 0 (t − t 0 ) and u(x, t) ≈ 0 hold when the incoming wave is much longer than the slope length. This can be viewed as the seabed slope works like a vertical wall against a relatively long wave. The convolution describes wave deformation over the slope under the observedwave boundary condition. Despite cyclic singularities, the kernel convolution is well defined, and the convolution can be 215 numerically performed, which will be presented in the next section.

Kernel convolution
Kernel representations derived in the previous sections allow us to predict coastal waveforms from observed wavedata apart from the shore. The prediction can be made instantly by kernel convolution of the offshore wavedata. This is advantageous over conventional ways of numerically solving the shallow-water equations, which involve both time and spatial integration 220 of the equations over a physical domain. Especially, 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.

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; where j is 0 or 1 (t − m = t 0m and t + m = t 1m ). To avoid midpoint singularities, the integration interval is divided into multiple These changes of variables will concentrate integration points of time near the singular points using the double-exponential function. Consequenty, 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 of t 0m = t 1m . Because the integrand decays with a doubleexponential rate, actual numerical integration is performed within a finite interval of p j . Therefore, the convolutions in (25)   240 and (26) are rewritten as dτ jm dp j dp j , dτ jm dp j dp j , where M is an integer ceiling of (t − t 0 )/4 when α = 0 and can be a smaller number when α > 0 because the kernels decay 245 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 ψ n dτ jm /dp j . We read the kernel table and then perform computation of the infinite integrals in several segments to predict a coastal waveform.

Validation
The convolution method is demonstrated and validated by comparing numerical results with an analytical solution for a simple 250 problem. Here, we consider a simple problem in which a monochromatic wave is observed apart from the shore that contains both incident and reflected waves. The observed wave is now given as 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.

255
To derive the equilibrium solution, we first take the Laplace transform of (32) and substitute it into (12). Then, we havê .
This can be inversely transformed into time domain by the residue theorem. However, the residues associated with the zeros of the modified Bessel function, s k ands k in (19), yield a transient solution because it accommodates the decaying exponential function of time with the damping factor α. Therefore, we need to consider only the residues of s = ±iλ for the equilibrium 260 solution :: in ::: case ::: of ::::: α > 0.
Consequently, the :::::: damped : equilibrium solution for water surface elevation,η(x, t), is given by: where the wave amplitude a(x) and the phase φ(x) are given by 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.

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 wavedata are available at different locations along the affected coastline (Kawai et al., 2013).
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 1,500 m and the depth contour is significantly curved. To apply the kernel method, 320 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 Figures 7b and 7c.
The numerical error was confirmed to be negligibly small.
Figures 8b-c compare the computed waveforms at the shoreline, ξ * , and the observed wavedata at P2 with the different values of damping parameters, respectively. The first term in (25), η * 0 (t − t 0 ), is also plotted to highlight the slope effect described 335 by kernel convolution. The comparison of ξ * and η * 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 wavedata 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 340 predict the coastal waveform from the observed data a few kilometres off the coast when the α value is appropriately chosen.

Case 2: Kamaishi and Ryoishi
Case 2 is located on the coast ranging from Kamaishi to Ryoishi, which are highly devastated areas 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 Figure 7, TM1 was located 76 km off the coast at a depth 345 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 is no data for onshore waveforms. Nevertheless, intensive post-event 350 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 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 355 GPS803 via the kernel convolution of the wavedata at TM1, assuming the water depth to be linearly decreasing from TM1 to the shore (β = 0.02). Figure 9b and 9c 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 occurs probably because the observed wave at TM1 is located far from the coast and contains wave components that do not propagate 360 shoreward. The resulting run-up height on the shore is sensitive to the choice of the damping parameter, as shown in Figure   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 wavedata at TM2 graphed in Figure 9e. The computed 365 results at GPS802 show better agreement with the observation than those predicted from TM1, as shown in Figure 9f. The results in 9c and 9f 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 ( Figure 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 370 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 wave profile at the shoreline from the wavedata at GPS802 that was available for a much longer period than 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 375 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).
The kernel representation of coastal tsunami evolution provides an efficient method for predicting onshore tsunami waveforms based on a single wave observation apart 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 380 time scale. In this case, the lead times for the real-time prediction of shoreline motions are 16 min 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 385 is well-predicted when the damping parameter is set in the range of 0.5-1.0.

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 390 from observed wavedata 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 apart 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 in the short past. The 395 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 400 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 405 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 semi-empirical kernel and to improve kernel formulation for better representation of tsunami decay in shallow water.
Data availability. The observation data of the 2011 Tohoku tsunami are available from the website of the Nationwide Ocean Wave Information Network for Ports and Harbour (https://nowphas.mlit.go.jp/pastdatars/data/NOWPHAS_Tsunami_data.zip).

410
Appendix A: Kernel convergence properties This appendix presents the convergence properties of the infinite series in (23). Firstly, we define the sequence of the series as Then, the infinite series can be rewritten as We look into the behavior of the sequence for large k. When k is large, we can use the approximation of c k in (20) and the following approximations as well; Furthermore, the Bessel functions can be approximated base on their asymptopic expansion for a large augument

420
which is valid only for large z. Since the augment of the Bessel function in the numerator of (A1) contains x 1 2 , we discuss the property of the sequence separately in two cases; (a) x = 0 and (b) x > 0.
(a) x = 0 Here we discuss the behavior of only f 0k (0, t) because f 1k (0, t) is not well defined. Since J n (0) = 1, we can approximate 425 (A1) for large k as 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 = 4m−2 such that 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. Substituing t = t ± m into (A7) and separting it into two parts, we have The orginal sequence is now expressed as sum of the two sequences. When t = t + m , either sequence vanishes and the other one forms a convergent series. Therefore, the infinite series is convergent at t = t + m . On the other hand, either sequence forms a convergent series and the other one generates a harmonic series when t = t − m . This confirms that the infinite series in (A2) slowly diverges to positive and negative infinity alternatively at t − m as shown in Figure 3b-c and 4.

445
Author contributions. The author confirms sole responsibility for study conception and design, analysis and manuscript preparation.
Competing interests. I declare I have no competing interests.