On the Resonance Hypothesis of Tsunami and Storm Surge Runup

. Resonance has recently been proposed as the fundamental underlying mechanism that shapes the ampliﬁcation in coastal runup for both Tsunamis and storm surges. It is without doubt that the resonance plays a rôle in runup phenomena of various kinds, however we think that the extent at which it plays its role has not been completely understood. For incident waves, the best approach to investigate the rôle played by the resonance would be to calculate the normal modes by taking radiation damping into account and then test how those modes are excited by the incident waves. There are a small number of 5 previous works that attempt to calculate the resonant frequencies but they do not relate the amplitudes of the normal modes to those of the incident wave. This is because, by not including radiation damping, they automatically induce a resonance that leads to inﬁnite amplitudes, thus preventing them from predicting the exact contribution of the resonance to coastal runup. In this study we consider two different coastal geometries: an inﬁnitely wide beach with a constant slope connecting to a ﬂat-bottomed deep ocean and a bay with sloping bottom, again, connected to a deep ocean. For the fully 1-D problem we 10 ﬁnd signiﬁcant resonance if the bathymetric discontinuity is large. For the 2-D ocean case the analysis shows that the wave conﬁnement is very effective when the bay is narrow. The bay aspect-ratio is the determining factor for the radiation damping. One reason why we include a bathymetric discontinuity is to mimic some natural settings where bays and gulfs may lead to abrupt depth gradients such as the Bay of Tokyo. The other reason is, as mentioned above, to test the role played by the depth discontinuity for resonance.


Introduction
During the last decades, several analytical and numerical studies of coastal runup were published (see Synolakis (1987), Carrrier et al. (2003), Kânoglu (2005), Kânoglu and Synokalis (2005), Özeren and Postacioglu (2012), and Stefanakis et al. (2015)) most of which made use of Carrier-Greenspan transformations (Carrrier et al. (2003)).Some of these works, for instance Stefanakis et al. (2015), Stefanakis et al. (2011) and Ezersky et al. (2013b), identified the resonance as the fundamental factor for the runup amplification.The bulk of the present study will be dedicated to determine those physical settings in which this might be the case.
In the past, several researchers looked at resonance aspect of the coastal runup.Among those, the ones that are the most relevant to the discussion in the present study are Stefanakis et al. (2015), Stefanakis et al. (2011), Carrier and Noiseux (1983), Ezersky et al. (2013a), Fuentes et al. (2015), Volker et al. (2010) and Yamazaki and Cheung (2011).The last two of these studies report coastal resonance mechanisms leading to amplified runups during the 2009 Samoa and 2010 Chile Tsunamis is the side view of the the incident wave at instant t = 0 .The black dot is the shoreline and (x s , z s ) are its coordinates with z s being equal to runup r .(B) is the geometry of the channel seen from above.In MODEL-1 the sloping channel is connected to a deeper channel of the same width (see broken lines in (B).In MODEL-2 the channel opens to a semi-infinite ocean.The width of the channel is 2a .Non-dimensional quantities are defined as x = x /L , ,y = y /L , z = z /(αL ) , η = η /(αL ) , D = D /(αL) and a = a /L .respectively.The work by Carrier and Noiseux (1983) does not explicitly mention resonance but their formulation clearly shows a constructive interference by multiple reflections of an obliquely incident Tsunami wave for a particular set of incidence angles.Some of the studies looked at the resonant mechanism in experimental settings (Ezersky et al. (2013a), Abcah et al. (2016)).
In our modelling, we will be considering monochromatic incident waves, simply because this makes it easier to shed light on the resonance.However, the mathematical algorithm we will develop is not limited to monochromatic waves but is capable of calculating runup for any kind of offshore source including the earthquakes, submarine landslides or atmospheric pressure perturbations.
This article deals with the transient runup response of a sloping channel (or bay) to an incident wave.Here the term transient is used for the wave evolution before the standing wave regime sets in.As such, the problem presents itself as an initial value problem.Such an initial value problem can be difficult and expensive to handle by purely numerical approaches.The reason for the difficulty is that, on the offshore boundary it is difficult to make a distinction between the incident and reflected waves.One way of surmounting this difficulty numerically is to take the computational domain so large that by the time the reflected waves arrive at the offshore boundary, the standing wave regime would have set in in the coastal region.For instance, Stefanakis et al.
(2011) did one-dimensional numerical simulations to understand the runup amplification by nonleading long waves.They cast the problem as a boundary value problem and imposed an offshore boundary condition, at distance L from the undisturbed shoreline, for the wave height as η = ±2η I 0 sin(ω t ) where t is time.However, this model, as will be explained in detail later, tends to overestimate the runup if the open boundary lies in one of the nodes of the standing wave that will eventually set in.When the incident and reflected wave distinction is made (see, for example, Antuono and Brocchini (2010)), the runup amplification factor (defined as r /(2η I 0 ) where η I 0 is the incident wave amplitude and r is the runup) remains finite as long as the frequency of the incident wave is real.This has independently been shown in Antuono and Brocchini (2010).
Our purpose in the present work is to determine the way in which the free modes near the coast are excited by the incident waves by taking the radiation damping into account.We will examine resonance in two different geometric settings: first a 1D slope which connects to a 1D channel with a flat bottom and then, again a 1D slope that connects to a semi-infinite 2D ocean with a flat bathymetry (see Figure ( 1)).When the wavelength of the incident wave is much shorter than the width of the sloping channel the completely 1D model is a good approximation to the natural case and we can neglect the geometric spreading of the waves at the toe of the sloping bay (this geometry will be referred to as MODEL-1 during the rest of the article).If this is not the case, then a 2D model (MODEL-2) near the mouth of the bay becomes necessary (see Figure 1).The Coriolis acceleration is neglected in both cases as we limit ourselves to small scales.Storm surges and meteotsunamis can excite Kelvin waves in the coastal areas and such waves travelling along the coast can come across bay mouths and force the normal modes of the bay, so the analysis here is not confined to one specific type of wave.
MODEL-1 is actually solvable, through fast Fourier transforms (see Ezersky et al. (2013b)), without necessarily resorting to the coastal free modes.However generalizing this solution approach to 2D in the deep ocean part is computationally very expensive, requiring a solution of an integral equation for each frequency component to calculate the transient response.Hence, the real importance of the technique developed in this manuscript becomes more apparent in MODEL-2 which can be of great engineering importance in modelling coastal amplifications of Tsunamis and storm surges in places like the Tokyo Bay which has a sloping bathymetry (Kataoka et al. (2013)).

Basic equations
MODEL-1 consists of a channel of constant slope α, length L that connects to a another channel of uniform depth αL + D (see Figure 1).Parameter D is the discontinuity in depth at toe of the slope.The maximum depth of the sloping channel is αL .
The governing equations we shall use over the sloping part of the geometry are non-linear shallow water equations (1) where u is depth-averaged velocity in offshore-pointing x direction, t is time, g is the acceleration due to gravity and η is the wave height.For the flat part of the domain, for MODEL-1, the linearized versions of the same equations, are used.Note that, when we shall eventually start to discuss MODEL-2, we will generalize equations ( 3) and (4) into two dimensions.Let us now define non-dimensional quantities as A hodograph transformation introduced by Carrrier et al. (2003), also called Carrier-Greenspan transformation (hereafter referred to as CG), uses "distorted time" λ and a potential ϕ defined as where The non-linear shallow water equations accordingly become (Carrrier et al. (2003)) To treat the incident wave problem, the initial conditions everywhere are and for the wave height η and the fluid velocity u respectively.We assume that both quantities are initially zero over the slope (0 < x < 1).The minus sign in ( 11) is due to the fact that the progressive wave advances in the negative-x direction towards the coast on the left.

Green's function and the free mode expansion
For the flat part of the domain (x > 1) we propose the following solution where ηI (ω) and R(ω) are the temporal Fourier transforms of the incident and the reflected waves (to be calculated) respectively.Here ηI ω is defined as The integrand in the right-hand side of ( 13), for all times, can be inferred from the initial condition as η I (t, x = 1) = η I (0, 1 + t √ D + 1).Note that the choice of the point x = 1 is completely arbitrary.Let us now propose a solution over the slope as where A 0 (ω) is the amplitude of the wave over the slope and J 0 is Bessel function of the first kind of order zero.
The corresponding wave height near x = 1 − (toe of the slope) is then given as Note that we performed a linearization here by taking λ = t and σ = √ x, because we assume that in the deeper part of the domain, waves are small.The unknown coefficients R(ω) and A 0 (ω) are to be determined from the continuity conditions at the toe of the slope (see A1 for the matrix form of these two equations).The zero subscript for A is used because for MODEL-1, having a fully one-dimensional geometry, only one coefficient is needed.Later, other coefficients will also be needed when we will consider two spatial dimensions.
The linearized free surface and flux continuity conditions at the toe (x = 1) are used to obtain A 0 (ω) and R(ω) as Here by linear boundary condition at x = 1 we mean that x = σ 2 , λ = t and depth is equal to x rather than x + η.
Both A 0 (ω) and R(ω) have simple poles at ω = 0.The non-vanishing poles are on the upper complex plane due to the conservation of energy through the radiation damping, because on these poles, the ratio A 0 (ω)/η I (ω) and R(ω)/η I (ω) both diverge and the wave can sustain itself without an incident wave and is fueled solely by the initial conditions (see Longuet-Higgins (1967) for a similar radiation damping problem over a circular submerged sill), these are called the free modes.Note that the conservation of energy requires that the amplitude of the incident wave is equal to that of the reflected wave with |η I (ω)| = |R(ω)| when the frequency is real.However this condition is relaxed when the frequency is no longer real because the energy density averaged over one cycle of oscillation evolves in time.See also Synolakis (1988) for a rigorous proof that the frequencies of the free modes are indeed on the upper complex half-plane for D = 0. Our argument based on the conservation of energy is more general and can be applied to any bathymetric profile.
Recently Stefanakis et al. (2015) attempted to calculate resonant frequencies for a geometry similar to Ezersky et al. (2013b) but imposed an offshore Dirichlet condition of the form η = 2η 0 sin(ωt) at a specified point x = L using a numerical wavemaker.With this boundary condition their results diverged at certain discrete real frequencies (equation ( 46) in Stefanakis et al. (2015)).They attributed this to runup resonance.However one can always adjust the frequencies in such a way that the reflected wave and the incident wave are off-phase by π at x = L, thereof leading to a destructive interference, requiring an infinite runup yielding to such a finite oscillation at x = L.There is no simple relation between the amplitude of oscillation at a specified point and the amplitude of the incident wave unless the incident and reflected waves are in phase at this particular point.Therefore relating the ratio between η 0 and the runup to the resonance is meaningless.This problem can partially be remedied by imposing an artificial relaxation zone (Stefanakis et al. (2015) and Madsen and Schäffer (2010)).In this work, we replace this approach with a physically realistic initial value problem where our definition of resonance is based on the ratio between the runup and the amplitude of the incident wave.Not only that the resonant frequencies we will calculate will be complex, but also that their real parts will be substantially higher than those found by Stefanakis et al. (2015) when the discontinuity, D, is small (see Table 1).Only in the limit of large D, the complex roots of the denominator (i of A 0 (ω) approach those of Dirichlet condition J 0 (2ω) = 0. Note that all of the resonant frequencies that we calculate in this work are independent of any nonlinearity.The reason for this is that we linearize the boundary conditions at the toe of the slope.The modes themselves (J 0 (2ωσ)), on the other hand, are affected by shoaling nonlinearity.Note that by the frequencies of the modes, we mean the frequency of the oscillation at the toe of the slope (x = 1).In section 4.2 we will show that for higher modes the shoaling nonlinearities are dominant over those created by the boundary condition at the toe.
Now consider an incident wave of the following form: where δ is the Dirac's delta function.The zero index of t in (18) relates to the phase of the incident wave.The response, ϕ, to such Dirac-type incident wave will be called Green's function G(λ, t 0 , σ).Remember that A 0 (ω)J 0 (2ωσ) is the temporal Fourier transform of ϕ .In ( 16) if we replace ηI by exp(−iωt 0 )/2π, then A 0 (ω)J 0 (2ωσ) will become the Fourier transform of the Green's function.This Green's function in Fourier domain is given as Any incident wave can be expressed in terms of a linear superposition of Dirac functions, the response, ϕ, will then be Because of the linearisation at the toe of the slope, the upper limit of the integration, t(λ, σ = 1), in (20) can be simply replaced by λ.We will use G(ω, t 0 , σ) to recover G(λ, t 0 , σ).Accordingly the potential will become where the integration over the frequencies will be transformed to a series of residues.This is not a closed integral over the complex plane but a real line integral between −∞ and ∞ except ω = 0 which we circumvent with an infinitesimal semicircle on the lower half-plane.The reason we use the lower half-plane is that we want the Green's function to vanish for negative values of λ.The whole integral can be cast into a closed integral by connecting ∞ to −∞ along a semi-circle on the upper complex plane and can be calculated using a residue summation.Here ω = 0 is not the only pole.As a matter of fact G(ω, t 0 , σ) has many poles in the upper half-plane (see equation 19).These poles are symmetrical with respect to the imaginary axis.Our aim is to understand the excitation of these free modes by the incident waves.
It is important to note that as x → ∞ these free modes diverge.However this is not a problem if one wishes to find a solution in the coastal zone.We will expand the reflected wave in terms of these free modes.For the numerical calculation, the free mode expansion is truncated at a finite term N .This finite series also diverges for x → ∞, thus the truncation error between the real reflected wave and the finite series expansion grows as x increases but this is also not a problem if one wishes to find the wave field near the coast because the discrepancy (or error) propagates towards offshore.
To perform the free mode expansion approach, let us rewrite the convolution (20) as In an effort to calculate solitary wave runup, Synolakis (1987) evaluated a similar integral in his linear approach.However his integral was tasked to compute directly the free surface, η, rather than the potential and consequently did not have a singularity at ω = 0.The Fourier transform of the particular solitary wave Synolakis (1987) considered approached zero so fast along the infinite-radius integration contour that he was able to close his contour on the lower half-plane (upper half-plane in his convention).Thus his complex integration loop did not contain any of the complex frequencies we mentioned above.The only poles that remained within his closed contour were those of the Fourier transform of the solitary wave he considered and were, therefore, independent of the geometry.His technique is limited to this particular incident wave forcing.In his article, Synolakis (1987), did express the solution in terms of a summation over discrete frequencies but these frequencies can not be interpreted as free mode frequencies not only because they are independent of the geometry but also because the result is only valid for times smaller than a critical time t c .To see this, let us write down the integral (given as equation 2.6 in Synolakis (1987)) in our convention for D = 0: where Φ is the Fourier transform of incident wave.In the integrant here, the exponential term exp(iωt) diverges on the lower semi-circle.This divergence is counter-balanced by the Fourier transform of the incident wave for t < t c but not later.
To do our calculations for wave evolution, we will need the non-vanishing poles of A 0 (ω) which we calculate using the Müller scheme (see Press et al. (2007)).For large D an asymptotic approximation is presented n the Appendix A. A comparison between these two approaches is displayed in Table-1.As seen in this table, as D increases, radiation towards offshore becomes less efficient, therefore making the complex parts of the roots smaller.The real parts, also decrease as D increases.This is because the waves that propagate from the shore towards offshore create reflecting waves continuously because of the variable depth of the slope over which they are travelling.However, we have radiation damping in our case, some of the waves that reach the toe of the channel escape from the sloping part of the channel.In the absence of radiation damping, all waves reaching the toe of the channel would reflect back.These waves contribute to standing waves of low frequencies over the sloping channel, because of the long distance they travel.Consequently, the relative weight of this low frequency component would increase if there is no or little radiation.In the case of any non-uniform bathymetry, at the high-frequency limit, ray theory can be used and the reflections will become minimal.This explains why the damping factor of the higher modes (imaginary parts of the eigenfrequencies) become larger.In the real geophysical settings where the discontinuities are less abrupt than the geometry depicted here (such as the edges of coastal shelves), the short waves will cross over the toes of the slopes, essentially without reflection.
The geometry considered by Stefanakis et al. (2015) features two consecutive slopes.If one fixes the slope angles in their work to a single value, the resulting geometry would be the same as ours.The frequencies that they would have come up with in their solutions would have been the real roots of J 0 (2ω).These frequencies are substantially smaller than the real parts of the frequencies we calculate (2.98 versus 2.4 for the fundamental mode) with radiation damping.Now let us return to the integral ( 22).This integral can be written completely in terms of λ and σ as: for λ > 2(1 − σ) because of the causality (ϕ becomes zero otherwise, because the disturbance takes 2(1 − σ) to travel from the toe of the slope to the point σ).A residue summation transforms (24) into: where J 0 and J 1 are the derivatives of J 0 (2ω k ) and J 1 (2ω k ) with respect to 2ω k .The first integral in (25) comes from the residue associated with the pole at ω = 0.It is calculated by using the fact that J 0 (0) = 1 and J 1 (0) = 0. Using the Leibnitz rule, the partial derivatives of the potential are then and The equation given in ( 9) is, in essence, the linear wave equation with cylindrical symmetry.Because of this, the partial derivative of its regular solution with respect to σ must be zero at σ = 0.A quick inspection of ( 27) reveals that the terms in the first and second line are not equal to zero individually, and their collective sum will involve a truncation error.This truncation error is subject to an amplification in the estimation of ∂ σ (ϕ)/σ which is used to calculate both η and u.To remedy this, we use the fact that, near σ = 0, the value of ∂ σ ϕ/σ is approximately equal to 2∂ 2 λλ ϕ due to L'Hôspital's rule.Hence, we use the numerical derivative of (26) to find the non-linear contribution to runup.

Resonance sensitivity for MODEL-1
In this section we consider a monochromatic incident wave of type η I 0 (ω) sin(ω(t + (x − 1)/ √ D + 1 )).We will study the evolution in the large time limit, namely the standing wave regime.The transient regime will be discussed in the next subsection.and 5 for the red dot-dashed curve.The black bullets correspond to the limiting runup given in the expression (37), normalized to ηI 0 , for the infinite slope.Here σ0 is taken to be one.
The resulting wave on the slope in the linear approximation will then be for t → ∞.This equation follows naturally from equation ( 15) where instead of evaluating the integral, we add the two contributions coming from ω and −ω.Taking into account that J 0 (0) = 1 the associated runup, r = η(t, σ = 0), becomes in the linearised theory.Here φ A is the argument of complex number A 0 (ω).In a real situation that would occur in the nature where the monochromatic incident wave gets generated at a particular instant t = 0, the expression The frequency of the incident wave is equal to the real part of the first-mode resonant frequency ω1 (see Table 1).The continuous blue curve is obtained from the series of residues (see equation 26).The dashed blue curve is the runup calculated using the Fast Fourier Transform approach.The green horizontal line on the top part of the figure is the limiting amplitude given by ω I |A0(ω I )|/η I 0 where ω I is the frequency of the incident wave.The red bullets are the runup produced by a wave-maker on an infinite constant slope.The action of the wave-maker is represented by the "tsunamigenic" seafloor motion given by h(t, x) = −2η I 0 δ(x−1) cos( ω1t)θ(t)/ ω1 where θ is the Heaveside function, h is the seafloor uplift and denotes the real part.in (29) provides a limiting value of amplitude of oscillation of the runup.This limiting amplitude, normalised to the amplitude of the incident wave, is displayed in Figure (2) as a function of 2ω where ω is the frequency of the incident wave.In this figure local maxima of the limiting amplitude of the runup can be observed for D = 1 and D = 5 but not for D = 0 where this amplitude steadily increases with the frequency of the incident wave.Therefore there is no resonance for D = 0.
In the limit of large D, the local maxima of the limiting amplitude of the runup occur at the frequencies where J 0 (2ω) = 0 and the value of the limiting amplitude of the oscillation of the runup at those maxima is given as where 2ω k is the k th root of Bessel function J 0 .Hence, the runup sensitivity to ω increases as D increases.
These results indicate that in the real world where D is often much smaller than one, the resonance is an insignificant phenomenon.For two-slope cases the resonance is even less significant because the imaginary parts of eigenfrequencies are decreasing functions of second slope, β.
One last remark relates to the power laws for runup, provided by Didenkulova et al. (2009).The Figure (2) shows that, when D is large, the denominator of ( 16) has roots that are closer to the real axis.This essentially means that for large D explicit relations for power laws might not be possible to derive.

Transient regime
The resonant phenomena we discussed above do not set in immediately upon the entrance of the incident wave into the slope region.It is important to know how fast the limiting amplitude of oscillation of the runup will be reached, because in a real situation the incident wave will have a finite duration, a fact that was not taken into account in the large time limit analysis.
For that purpose the residue series in equation ( 26) must be evaluated.Here, the incident wave is, taken as η I  1 for the values of ω1.All curves are obtained using series of residues.Note that in the bottom figure the wave arrives at the shore at t = 2 and the first reflection from the toe of the slope reaches the shore at t = 6 for which the continuous plot includes a slight glitch.In the main text we elaborate on this glitch.do not "switch on" at a given time, at once.Therefore their initial profiles do not have discontinuous spatial derivatives.As a matter of fact, such discontinuities will trigger short-frequency waves which can not be properly modeled using shallow water approach anyway.The discontinuity mentioned here is seen much more clearly in the shoreline velocity field in Figure ( 5).

Nonlinear effects
We assume the incident wave to be linear.Considering that we are dealing with a monochromatic incident wave, this makes sense.Because any nonlinearity over the flat part of the ocean would have generated higher harmonics during the propagation.
As long as the waves do not break, the nonlinearity arising from the shoaling over the slope is accounted for by the CG approach.This particular nonlinearity, as indicated by Pelinovsky and Mazova (1992), does not affect the maximum shoreline velocity, it does, on the other hand, affect the timing of the maximum.After the CG transformation the equations become linear, the only remaining nonlinearity in the (λ, σ) space for the runup is the u 2 s /2 term in (7) with u s being the fluid velocity on the shoreline.Figure (6) shows the difference between the two values of runup calculated by taking and not taking this u 2 s /2 term into account.Another subtelty in this figure is that, for the nonlinear case the time is given as t = λ + u s .
For t → ∞, the nonlinear runup is given as ω|A 0 | sin(ωλ) − 0.5u 2 s which approximately becomes The reason we performed this expansion with respect to small shoreline velocity is that we eventually want to compare the nonlinear effects on the runup, caused by the nonlinear boundary condition at x = 1 with the contribution by the nonlinear terms in the above equation.
At the near-resonant frequencies nonlinear effects will become important even in the deep part of the slope, rendering the linearized boundary conditions that we apply at x = 1 invalid, compromising the relations for A 0 (ω) and R(ω) (equations ( 16) and ( 17)).In Appendix B we give the results from a first-order perturbative approach to correct the boundary conditions for nonlinearity.We performed a comparison between the nonlinearity brought in by the boundary condition and that caused by the nonlinear part of (32).The details of this perturbative analysis can be found in Appendix B. The results suggest that, the first-order correction associated with the nonlinearity of the boundary condition at x = 1 is given by where the superscripts refer to the first order perturbation.When D becomes very large, at the near-resonant frequencies (J 0 (2ω) ≈ 0), the time-dependent and time-independent parts of the solution for η over the slope, respectively become These indicate that the nonlinearity brought in by the nonlinear part of ( 32) is comparable to that caused by the boundary conditions at x = 1 only for the first mode and when D is infinite.This is to be expected, because shoaling is insignificant for the first mode.The wave must travel a distance of several wavelengths in order to experience significant amplification associated with shoaling.For the first mode (2ω ≈ 2.4) the analysis yields N (1) ≈ −0.41|A 0 | 2 and C (1) ≈ 0.097.The coefficient of the second harmonic in (32) is 1.4|A 0 | 2 .For the higher modes the dominance of the (32), nonlinear effect associated with shoaling, increases rapidly because of the fourth power of ω.

Resonance for infinite slope
In this section we investigate the resonant frequencies of the waves produced by a wave-maker placed on an infinite, constant slope.The reason for this practice is that the work that claims significant resonance (Stefanakis et al. (2011)) considered an infinite slope.In the analysis we will allow the waves to progress in the offshore direction, unrestricted.When the the wavelength produced by the wave-maker matches the distance of the wave-maker to the shoreline one might expect a resonance to occur.To see if this is really the case, let us go into a little further detail about the nature of the wave-maker.When the water is sufficiently deep, in the vicinity of the wave-maker linear shallow water equations will apply.The effect of the wave-maker which starts its action at a given time, will be equivalent to hypothetical volume injections and suctions at a rate 2η I 0 √ depth sin(ωt) = (2η I 0 σ 0 ) sin(ωt) where σ 2 0 is the horizontal distance between the wave-maker and the shoreline.Consequently, two waves of equal amplitude of η I 0 start moving in opposite directions if ω is not smaller than σ 2 0 .For ω → 0 the asymmetry introduced by the slope becomes significant and the wave maker radiates essentially in the offshore direction (see the black bullets Figure 2).Due to the superposition principle, the action of the wave-maker will not prevent the waves reflected from the shore to freely propagate in the offshore direction.The free surface response to such a source on the slope is given in equation (2.13) in Özeren and Postacioglu (2012).To adapt the source in Özeren and Postacioglu (2012) to our case, the expression S(σ, λ) needs to be replaced by 2η I 0 sin(ω λ)δ(σ − σ 0 )σ 0 where σ and λ are integration variables.In Figure (3), the red bullets show the runup associated by such a source over an infinite slope as a function of time.For t → ∞ a simple analytical solution that satisfies the radiation condition (wave progressing in +x direction for x >> x 0 ) can be found for a source of the form mentioned above.The presence of the source at x 0 induces a discontinuity in the fluid velocities such that where x 0 is both the horizontal position of the source and the non-dimensional depth at that position.The corresponding solution as t → ∞ is then given as where Y 0 is the Bessel function of the second kind.For σ >> σ 0 the relation given in ( 36) is proportional to sin(ω(t The amplitude of the resulting runup then becomes and for large frequencies it is reduced to 2 √ πωσ 0 η I 0 .The black-dotted curve in Figure 2 shows this expression normalized by η I 0 as a function of 2ω for σ 0 = 1.As seen in this figure, the resonance for the infinite slope is almost non-existent.It is also worthwhile to mention that the runup tends to zero for ω → 0 for the infinite slope.The reason for this is that the wave generated during one period of the forcing activity spreads over to an infinite distance offshore.It is clear from Figure 2 that the wave generated by the artificial source/sink that we introduce emulates satisfactorily an incident wave if the wavelength of the generated wave is less than the distance to the shore.

Normally incident wave from a 2D ocean into a bay
Analytical studies that consider runup in 2D are rare.There are, however, some studies that combine 1D analytical approaches with 2D numerical simulations such as Choi et al. (2011).In this section we shall study the MODEL-2.We consider an incident wave of the form η I 0 (ω) exp(iω(t + (x − 1)/ √ D + 1 )) coming from the open ocean into a sloping bay of width 2a.This wave will be reflected back by the the shallower part of the bay.This reflected wave will then be subject to geometrical spreading in the open ocean.Consequently, the one-dimensional nature of the waves will be lost in the deeper part of the bay.We will now derive a simple mathematical formulation to model these phenomena.
The governing equations, linearized in the deeper part of the channel and the open sea, for the potential, ϕ, read where h(x) = x is the undisturbed depth in our case.For a constant-slope beach Mei et al. (2004) page 155 and Zhang and Wu (1999) used the confluent hypergeometric function M .We will now adapt their solution to a sloping bay.Let us start by defining a potential inside the channel as: where the unknown coefficients A 0 , A 1 , .. are to be determined from the boundary conditions at the mouth of the channel.Note that the solution given in (39) satisfy the linear shallow water equations in the deeper part of the channel.For 1 − x >> a the y-dependent part of the solution quickly decays.The cosine term with integer values of n is both symmetrical about the x-axis and its derivative with respect to y for y = ±a is zero, satisfying, therefore the no-flux condition on the channel sides.In the the deeper part of the channel where 1 − x is of the order of a, the velocities in x and y directions can be derived from the potential given in equation( 39 where î and ĵ are unit vectors in x and y directions respectively.Note that the Hankel function H (2) 0 satisfies both the linearized wave equation in two dimensions and the radiation condition for ω > 0. For ω < 0, on the other hand, Hankel function will be replaced by its complex conjugate alongside the rest of the terms.The integral in equation ( 40) represents the potential of the waves radiating from the mouth of the channel.Function S(ω, ỹ) is the unknown virtual sources distribution along the mouth of the channel.This source distribution will be determined by matching the potentials given in equations ( 39) and (40).
We want the solution in the open sea to satisfy the no-flux condition at (x = 1, |y| > a).According to Mei et al. (2004) page 194 , the solution given in (40) satisfies this condition as To match the depth-integrated u along the mouth we need to satisfy The second condition to be satisfied is the continuity of η itself, across the mouth.This condition reads Here it is important to note that the condition ( 43) is an integral equation for the source distribution.The condition (42), on the other hand, is not an integral equation because the field inside the channel is governed using summations rather than an integral.In order to solve a system that simultaneously involves an integral equation and a set of algebraic equations we shall expand the source distribution, S(ω, y), in terms of even-order Legendre polynomials as where the even indices of the Legendre polynomial were used to make sure that the solution is symmetric with respect to the symmetry axis of the sloping channel (y = 0).The coefficients Sn (ω) in ( 44) and A 0 , A 1 , ..., A N −1 in (39) (we truncate the series in (39) to N − 1) are solved in such a way to minimize the following penalty integral along the mouth: Thus we have 2N unknowns.The integral ( 45) is evaluated numerically using Gauss quadrature.For precision, the number of quadrature points we use is larger than 2N .The equations ( 42) and ( 43) have to be satisfied at all quadrature points, making the resulting system over-determined.We use weighted least-squares approach to solve this over-determined system.Algebraically speaking, this simply corresponds to writing down equations ( 42) and ( 43) for each quadrature point and multiplying each equation by the square-root of the corresponding quadrature weight and solving the resulting linear, over-determined system by using conventional least-squares method.
Note that in equations ( 39) and ( 44 vanish for n = 0.The conservation of the net flux requires that −2aωA 0 J 1 (2ω) = 2a S0 where J 1 is Bessel function of order 1.When the frequency of the incident wave is equal to that of the free oscillation of the system, the coefficients A 0 and S0 should diverge.In order to determine the frequencies of these free mode oscillations (natural frequencies), the roots of 1/A 0 are sought in the complex plane using the Müller method.
A quick look at Table 2 reveals that for any mode, with the decreasing channel width (2a), imaginary parts of the normal mode frequencies decrease.However, for lower modes this effect is slightly more pronounced.This is because the time necessary for the waves to travel over the open ocean across a distance of half channel width a is equal to a/ √ D + 1 and when this time is much less than the period of a longitudinal oscillation within the sloping channel, the geometrical spreading in the open sea is very efficient.This makes the free surface in the vicinity of the channel mouth almost flat due to the fast escape of the waves, rendering the channel mouth boundary condition effectively a Dirichlet condition with η = 0. Consequently, the waves reaching the mouth reflect very efficiently back towards the shore, limiting the radiation damping.
The rays that the waves follow in the sloping channel are straight lines at the shallower parts and they bend towards the corners as they get near to the mouth of the channel, due to geometrical spreading (see the streamlines in Figure C.1).If the width of the channel increases, this corner effect will penetrate deeper into the channel, making the rays longer.Longer rays will decrease the frequencies of the free oscillations.This can be observed in Table 2. Overall, because of this ray bending the frequencies in MODEL-2 are lower than those in MODEL-1.Having said that, for a fixed value of channel width, if D increases, the discrepancy between MODEL-1 and MODEL-2 decreases because for D → ∞ the free surface in open ocean side of the domain will be flat for all free modes for both MODEL-1 and MODEL-2, making the behaviour the same within the channel.Now let us turn our attention to the transient response for MODEL-2.As before, we shall model the transient response to an incident wave of the form η I (t + x/ √ D + 1) and, as the first step, we shall look at the particular case in which η I is a Dirac's delta function.For the sloping channel, such a response for the potential ϕ corresponds to the integration of equation ( 39) with respect to ω.This integration reduces down into a residue summation which we shall explain below.Since we are not interested in the y-direction dependence of the runup, we are only interested in A 0 .But since they are inter-related through the boundary conditions, we do end up having to calculate A 1 , A 2 , ....A n as well, even though we do not need them individually for our physical interpretations.
Remember that in MODEL-1 we had an analytical expression for A 0 , therefore also for ϕ (see equations ( 16) and ( 14)).For MODEL-2, on the other hand, we do not have a closed-form relation for A 0 except when ω → 0. Obviously, when ω → 0, the free surface over the slope becomes flat and in this regime A 0 (ω) is approximately equal to 2η I (ω)/iω for both MODEL-1 and MODEL-2.For the rest of the frequencies, ω k , that make A 0 singular, we simply calculate circular integrals around each ω k on the complex plane.Consequently, for ϕ we have where c k are integral contours with infinitesimal radii around each ω k .This approach is much faster than calculating the frequency integrals exclusively along the real axis.This is because, on the real axis the integrand becomes oscillatory and the accuracy can only be sustained using small integration steps (in Model 1 we used several tens of thousand integration points to keep the integral stable).This is not feasible because for each integration point, an integral equation (equation ( 45)) needs to be solved to calculate the relevant A 0 .On the other hand, for each ω k , the complex integral A 0 (ω)J 0 (2ωσ)dω can be satisfactorily calculated using just four points because the integrand for ω → ω k becomes proportional to 1/(ω − ω k ) which is not oscillatory.A further practicality of this approach for an operational activity such as predicting a storm surge runup is that these contour integrals are independent of the structure of the incident wave and therefore can be pre-calculated for a particular bay or channel geometry and be injected in during the operational calculations.To relate the runup to 47 we follow the same procedure as in MODEL-1 to compute the partial derivatives of ϕ.In Figure ( 7) the runup, shoreline velocity and shoreline acceleration is given as a function of time.
In Figure C.1 the rays not only bend towards the mouth of the channel but also they coalesce.This coalescence reflects certain features that relate to the energy exchange to and from the channel.
To summarize, for both MODEL-1 and MODEL-2, when the incident wavelength is larger than the channel length, the runup amplitude tends to 2η I , in other words the wave becomes blind to the bathymetric variation.However an interesting case arises in MODEL-2 when the wavelength is much larger than the width of the channel but is still less than the length of the channel.In this case an analytical solution can be found using a conformal mapping approach (Mei et al. (2004) page 218 for flat-bottomed channel case).In Appendix C we give the mathematical details of the analysis in which we generalize the conformal mapping approach for the sloping channel with D = 0.As seen in the Appendix, the corresponding A 0 (ω) becomes where Γ ≈ 1.781 is the exponential of the Euler-Mascheroni constant.The runup amplitude is simply then given as ω|A 0 (ω)|.
The Figure (8) shows this value as a function of twice the frequency of the incident wave, calculated using both the integral equation and conformal mapping approaches.Kajiura (1977) also used conformal mapping to reach the same conclusion (note that his equivalent figure includes the amplification factor rather than the runup, therefore his values are half of those reported here).In this work we essentially extended the solution of Kajiura (1977) to short waves using integral equations.47)) normalized by the incident wave amplitude as a function of time.The channel width, 2a is equal to 0.2.The incident wave is given as η I 0 sin ( ω1(t + (x − 1))) θ(t+(x−1)) where ω1 is the lowest free mode frequency for MODEL-2 (see Table 2 for D = 0), smoothed by multiplying it with the tanh fuction as in Figure ( 5).The continuous red curve is the shoreline velocity.The dashed green curve is the shoreline acceleration.The horizontal blue line is the maximum runup calculated using conformal mapping.

Conclusions
In this work we studied the resonance aspect of the coastal runup as a response to incident waves.The analysis follows a normal mode approach and examines the sensitivity of those normal modes to a given incident wave to produce coastal runup.
In MODEL-1, significant runup sensitivity, in other words resonance, occurs only when D is large.Large values of D are not encountered very often in the nature, not even in the shelf breaks.In the two-dimensional open ocean and finite-width sloping channel case (MODEL-2) resonance occurs when the aspect ratio of the bay (width/length) is small.This kind of bay (or channel) geometry exists in many coastal regions such as the Bay of Tokyo, making the results relevant in engineering practice.The residue method developed here can actually be generalized for more complicated channel geometries (such as piecewise constant slopes with varying width) by performing a "fusion" of this method with the boundary elements technique.This is because boundary elements technique is recently proving very efficient for solving Helmholtz equation in multiple dimensions (Gumerov and Duraiswami (2008), Takahashi and Hamada (2009)).So the fusion should work as follows: for a series of complex frequencies, the boundary elements determines those frequencies for which the response diverges (for a single run, 5 the boundary elements, formulated in the temporal Fourier domain, can only calculate the response for a monochromatic incident wave).These are the normal mode frequencies for the particular geometry chosen.Then for any incident wave train (this can be, for instance, buoy data time series for an incoming Tsunami) we can calculate the contour integrals in (47) to compute the response much faster than inverse Fourier transforming the boundary element results.The operational FFT-based models that use buoy data for Tsunami warning (such as Lin et al. (2014)) can not assimilate the buoy data continuously as it Note that this expression does not give information on the transients, it is only valid as t → ∞.An incoming wave arriving from far-away, considered in this work would probably be highly linear, because otherwise it would have lost its monochromaticity due to nonlinear effects.However, near the resonant frequencies the amplitudes of oscillations over the slope (even the deeper part further off-shore) are generally much larger than the amplitude of the incident wave.Nonlinear effects might therefore be important over the entire slope, including the toe region.On the toe, the upper bound of the free surface oscillations, cannot be larger than twice the incident wave amplitude.But the derivative with respect to x of η can be very large because of the discontinuity of u for large D and a theoretical upper limit on it cannot be found.Therefore, the linearization of the boundary conditions might be a serious compromise.
Our aim is to calculate an additional term, ϕ (1) , such that ϕ = ϕ (0) + ϕ (1) (B2) will satisfy the boundary condition at the toe, to the order |A 0 | 2 .The linear approximation to ϕ (0) is given as by using which one can find the linear expressions of the free surface and the fluid velocity.The term that is responsible for the nonlinear correction of the boundary condition at x = 1 is the difference between the ϕ (0) (t, x) and ϕ (0) lin (t, x) near the toe of the slope, denoted as ∆ϕ (0) .To the first order, this is simply obtained through expanding ϕ (0) (t, x) around x = 1 and t : x cos(ω t) sin(ω t) for x → 1 − (B4) and ∆ϕ (0) ≡ 0 for x > 1.We denote the free surface associated with ∆ϕ (1) as η (1) .Since η (0) associated with ∆ϕ (0) has a discontinuity of the order |A 0 | 2 at the toe, this discontinuity will be remedied by adding η (1) = ∂ t ϕ (1) .The boundary condition at the toe then reads For small |z|, f −1 can be approximated by Matching this with the inner solution A 0 (ω)J 0 (2ω √ x) ≈ A 0 (j 0 (2ω)+(x−1)J 1 (2ω)) at x = 1 − (note that x = 1 − means that the distance to the mouth is much smaller than the channel length but much larger than a, the width of the channel), we obtain The outer solution given by equation ( 40) is, at intermediate range, (that is the distance from the mouth of the channel is much larger than the width of the channel but much smaller than the wavelength) given as In obtaining C8, we made the approximation that a << (x − 1) 2 + y 2 , therefore the Hankel function term in 40 can be taken out of the integral and approximated using the logarithm, again assuming that ω << (x − 1) 2 + y 2 .For |z| → ∞ ,in the upper complex plane f −1 (z) can be approximated by: x

Figure 2 .
Figure 2. The limiting amplitude of runup, rω, normalised to the amplitude of the incident wave is shown as a function of the frequency of the incident wave multiplied by 2 for MODEL-1.The depth discontinuity D is 0 for the blue continuous curve, 1 for the green dashed curve,

Figure 3 .
Figure 3. Continuous and dashed '-.' blue curves display the runup normalised to the amplitude of incident wave for Model I with D = 0.

Figure 4 .
figure it takes longer for the runup to reach its limiting value.This is due to the fact that for larger D the imaginary part of the frequency of the free mode is smaller.This means that if the incident wave excitation were to be cut at a given time, the standing wave oscillations over the slope would last longer (before eventually decaying due to the radiation) when D is larger.Similarly, for large D, it takes also longer for the standing wave regime to reach its limiting amplitude.It is also important to note that in both part of Figure (4) the initial time derivative of the runup is very large, leading to the glitch at t = 6 for the D = 5 case on the bottom.No such glitch exists D = 0 because there is minor reflection from the toe.In the real nature, waves

Figure 5 .
Figure 5.The top figure is the shoreline velocity us normalised to the amplitude of the incident wave as a function of t (Model 1).The frequency of the incident wave is ω1.The continuous green curve is the shoreline velocity corresponding to incident wave η = η I 0 sin( ω1(t + (x − 1)/ √ D + 1 ))θ(t + (x − 1)/ √ D + 1).The red, dashed curve is the same for the incident wave that has been smoothed by multiplying it by tanh(t+(x−1)/ √ D + 1).The bottom figure is the time derivative of the shoreline velocity for the smoothed incident wave.

Figure 6 .
Figure 6.Runups normalized to the amplitude of the incident wave for MODEL-I.The incident wave is in the form η = η I 0 sin(ω(t + (x − 1)/ √ D + 1 ))θ(t + (x − 1)/ √ D + 1) tanh(t + (x − 1)/ √ D + 1).The frequency of the incident wave is ω1(D = 0).The blue continuous curve is the linear runup.The dash-dotted green curve includes nonlinear effects associated with CG transform (r = ∂ λ ϕ − 0.5u 2 s and t = λ + us).The amplitude of the incident wave is 1/15.Nonlinear effects due to the boundary condition at the toe of the slope have been neglected.
) using u = −∂ x ϕ and v = −∂ y ϕ respectively.The bulk of the incident wave will be reflected back by the solid boundary at x = 1, |y| > a. Consequently the wave in the open ocean will be η = 2η I 0 cos(ω(t − (x − 1)/ √ D + 1 )) perturbed by the waves radiating from the mouth of the channel and the scattering from the channel mouth corners.The potential in the open sea then reads ), the only terms that are responsible for the net flux from the channel to the open sea are A 0 and S0 since the following integrals a

Figure 7 .
Figure7.The blue continuous curve is the runup (calculated using equation (47)) normalized by the incident wave amplitude as a function of

Figure 8 .
Figure 8.The maximum runup normalized to the amplitude of the incident wave for the standing wave case for two different channel half-width values.The continuous blue curves are computed using the integral equation and the red broken curves are obtained from the conformal mapping (see equation (48)).The green dot-dashed curve is the maximum runup for the infinitely wide channel.In the top figure the half-width of the sloping channel is 0.1, and in the bottom figure the half-width is 0.2.

=Figure C. 1 .
Figure C.1.The blue curves are x + iy = z = f −1 (z = ζ + ni) for ζ varying between -5 and 5 and n = 0.06/n.The three blue curves correspond to n = 1, 2, 3.The smaller the parameter n becomes, further left these blue curves will reach (the left extremities of these curves are approximately 2aπ ln( n) where a is 0.2).The red curves are the streamlines.

Table 1 .
The complex natural frequencies multiplied by 2 using both Müller method and the asymptotic approach (see (A5)) for MODEL-1 are tabulated.Note that 2ω k tend to the roots of J0 for large D. All nonlinearities are neglected in calculating these frequencies.

Table 2 .
The complex natural frequencies multiplied by 2 are tabulated for Model 2. These are essentially those 2ω k that are the roots of 1/A0(ω).In this table those lines with "(int)" list the frequencies calculated using A0(ω) obtained from minimising the penalty integral (45).The lines with "(conf)", on the other hand, use (48) from the conformal mapping formulation for this quantity.Note that 2ω k tend to the roots of J0 for large D or small a.