Long Wave Run-Up Resonance in a Multi-Reﬂection System

: Wave reﬂection and wave trapping can lead to long wave run-up resonance. After reviewing the theory of run-up resonance in the framework of the linear shallow water equations, we perform numerical simulations of periodic waves incident on a linearly sloping beach in the framework of the nonlinear shallow water equations. Three different types of boundary conditions are tested: fully reﬂective boundary, relaxation zone, and inﬂux transparent boundary. The effect of the boundary condition on wave run-up is investigated. For the fully reﬂective boundary condition, it is found that resonant regimes do exist for certain values of the frequency of the incoming wave, which is consistent with theoretical results. The inﬂux transparent boundary condition does not lead to run-up resonance. Finally, by decomposing the left- and right-going waves into a multi-reﬂection system, we ﬁnd that the relaxation zone can lead to run-up resonance depending on the length of the relaxation zone.


Introduction
Long wave run-up is difficult to observe in nature in real-time due to the catastrophic consequences it usually leads to. An illustration of a long wave is a tsunami, which has a wavelength far greater than the water depth and a small amplitude in the open ocean. When a tsunami reaches the nearshore region, its amplitude increases significantly, while its propagation speeds decreases. A large part of the theoretical results for run-up are based on shallow water theory. Analytical linear solutions for periodic long waves on a slope can be found, for example, in Green and Ferrers [1], Lamb [2], and Keller and Keller [3].
Tsunamis are known for their dramatic run-up height. From the data collection of post-tsunami surveys, extreme run-up values have been observed that exceed 20m and the most severe damages were not always caused by the first wave [4][5][6]. Extreme run-up values cannot be explained by classical analytical solutions. A possible explanation for these extreme run-ups is resonance. Miles [7] and Wilson [8] discussed the conditions for resonance in closed and open basins. The Helmholtz mode was also studied, which is the most important mode in natural basins and is normally observed in bays, inlets, and harbors with narrow entrance. Kajiura [9] described the possibility of resonances in bays. Rabinovich and Leviant [10] introduced the concept of shelf resonance, which is also called quarter-wavelength resonance. Neetu et al. [11] found that reflection and refraction from the nearshore topographic features can lead to the amplification of run-up values of non-leading tsunami waves. Resonance effects in the process of run-up have been studied numerically by Stefanakis et al. [12,13] and experimentally by Abcha et al. [14,15].
Obviously boundary conditions play an important role in numerical simulations: they can easily make the difference between successful and unsuccessful computations, or between fast and slow ones. There is still no clear understanding of how to impose both incoming and outgoing waves at a boundary. One generally uses an open boundary condition seaward, which is a numerical artifact that lies in a limited computational domain to solve a mathematical problem that has no boundary at all. Such open boundary conditions are usually called "radiation", "absorbing", "non-reflecting", or "far-field" boundary conditions. They all allow waves to leave the truncated domain, thus avoiding spurious reflections that may pollute the solution in the interior of the computational domain of interest. They can be classified into two groups: non-reflecting boundary conditions and non-reflecting boundary layers. Non-reflecting boundary conditions absorb impinging waves on the artificial boundary. A classical example is the Sommerfeld boundary condition [16,17]. On the other hand, non-reflecting boundary layers have the property of absorbing waves that are traveling inside the layer. Open boundaries are considered in a widespread literature. Examples are Jensen [18], Carter and Merrifield [19], Ma and Madsen [20].
However, in a laboratory wave flume, not all reflections caused by a wave paddle can be eliminated. The resulting wave system is a multi-reflection system. In the framework of wave energy converters, Wei et al. [21] discussed the re-reflection effects on the wave field. In nature, waves on some topographies are easily trapped, such as in narrow bays. When the water depth changes abruptly in shelves with piece-wise geometry, a multi-reflection system can also occur. Such reflective boundary conditions are seldom discussed.
Run-up on a sloping beach with the use of an open boundary condition at the seaward boundary matches the linear analytical solution that shows no resonance at all. However, no matter if it is theory or experiment, there exist some resonant frequencies [2,13,22]. Wave trapping is one of the phenomena that lead to resonance [23,24]. Using reflective boundary conditions leads to the reflection of outgoing waves and leaves them in the computational domain, thus constructing a multi-reflection system. We discuss run-up in such a system.
In this paper, we provide the analytical solution of run-up on a bounded domain that consists of a sloping beach connected to a flat bottom (linear shallow water equations). Conditions for run-up resonance are given. Then, we investigate numerically the same problem (nonlinear shallow water equations) using three different boundary conditions: influx transparent boundary conditions, reflective boundary conditions, and relaxation zones. The outcome of the numerical simulations is shown. For the fully reflective boundary condition, it is found that resonant regimes do exist for certain values of the frequency of the incoming wave, which is consistent with theoretical results. The influx transparent boundary condition does not lead to run-up resonance. By decomposing the left-and right-going waves into a multi-reflection system, we find that the relaxation zone can lead to run-up resonance depending on the length of the relaxation zone. Finally, the applicability of the results to tsunami science is briefly discussed.

Materials and Methods
We consider periodic waves of angular frequency ω incident on a linearly sloping beach with angle α (see Figure 1). The sloping beach starts at a distance L from the shore.
The one-dimensional (1D) nonlinear shallow water equations (NSWE) read ∂η ∂t where x is the horizontal coordinate, t is time, η(x, t) is the free-surface elevation, h(x, t) = h 0 (x) + η(x, t) is the water depth with h 0 (x) the still water depth, u(x, t) is the depth-averaged horizontal velocity, and g is the acceleration due to gravity. The origin of the coordinate system is at the shoreline. The still water depth h 0 (x) is equal to L tan α for −L t ≤ x ≤ −L and −x tan α for −L ≤ x ≤ 0.

Theoretical Model for Run-Up Resonance Based on the Linear Shallow Water Equations
Resonance phenomena play an important role in run-up amplification. They can explain extreme run-up heights and later arrival of waves with larger amplitude. For the geometry of Figure 1, Billingham and King [25] used linear shallow water theory to show that a resonance occurs when 2ω L/(g tan α) is a zero of the Bessel function J 0 . This resonance can also be found both in the numerical results of Stefanakis et al. [13] and in the experimental results of Ezersky et al. [22], because of the presence of a reflective boundary condition seaward.
Even though the numerical simulations shown in the next section are based on the 1D NSWE (1) and (2), we use the 1D linear shallow water equations (LSWE) to explain the resonance phenomenon analytically. The 1D LSWE read ∂η ∂t The phase velocity of the wave is denoted by c(x) = ω/k(x) = gh 0 (x).
Equations (3) and (4) can be combined into a single equation: For a beach with a constant slope, we look for a harmonic solution of the form [26] η( Inserting Equation (6) into Equation (5) leads to where σ = x 0 k(x) dx with k(x) the wave number. Equation (7) is the Bessel equation of the first kind, whose solution can be represented by where A 1 and A 2 are linear coefficients, and J 0 , resp. Y 0 , are the Bessel functions of the first and second kind. As Y 0 tends to infinity at the shore, the coefficient of Y 0 must be zero in the nearshore region if we assume that the amplitude of the free-surface elevation is finite. Therefore, A is expressed only in terms of the Bessel function of the first kind: where R max denotes the run-up. Billingham and King [25] found the resonance mentioned above by imposing a fixed amplitude at x = L. Indeed, σ(x = L) = 2ω L/(g tan α). For a constant-depth region, the wave field can be described as a sum of incident and reflected waves: where A i and A r are the amplitudes of the incident and reflected waves, respectively. Next, we discuss run-up due to the more realistic bathymetry profile that consists of a constant depth region connected to a sloping beach. Hereafter, it will be referred to as the canonical case (see Figure 1). Following Kânoglu and Synolakis [27] and Stefanakis et al. [13], we use the continuity of the free-surface elevation and of the horizontal fluxes to match surface elevations (Equation (11)) and particle velocities (according to Equation (4), it is equivalent to Equation (12)) at the toe of the sloping beach from two adjacent segments, and match the surface elevations at the left boundary Equation (13), respectively. We obtain where The determinant of the system (11)-(13) is given by When the determinant vanishes, run-up becomes resonant. Therefore, the resonance condition reads det(system) = 0. Clearly, the distance L t from the shore to the seaward boundary that appears in σ 3 plays a crucial role in the resonance condition.

Numerical Integration of the Nonlinear Shallow Water Equations-The Boundary Conditions
We integrate the 1D NSWE (1) and (2) with a finite volume solver based on the Characteristic Flux scheme with UNO2 reconstruction for higher order terms. A third-order Runge-Kutta scheme is used for time discretization. The numerical model is described and validated in Dutykh et al. [28]. What is new here is the fact that we test three types of boundary conditions described below.

Influx Transparent Boundary Conditions
Appropriate boundary conditions are needed to let the wave travel in and out of the domain in an undisturbed way. Influx transparent boundary conditions are boundary conditions that make waves on boundaries propagate into the computational domain and are transparent for wave propagation out of the computational domain. In the region of constant depth, the LSWE reduce to two wave equations [2]: and where Φ is the time-integral of the displacement past the plane x, up to the time t: Φ(x, t) = t u dt.
Equations (18) and (19) have the general solutions where the subscripts lg and rg denote the waves traveling to the left and to the right directions, respectively. Differentiating Equation (21) with respect to x and t yields To build a left boundary at x = x l that is transparent for waves propagating to the left but allows waves to propagate to the right, we eliminate Φ lg : This is called the left Influx Transparent Boundary Condition (ITBC).
Moreover, from (3), Therefore, the velocity at the left boundary can be written as Similarly, the right ITBC (x = x r ) that is transparent for waves propagating to the right but makes waves propagate to the left can be written as Therefore, the velocity at the right boundary is The head-on collision of two single waves on a flat bottom with the use of the left and right influx transparent boundary conditions is shown in Figure 2. Waves of the form η 0 (t) = 0.1 sech 2 (0.17 t − 2) are generated at the boundaries. As the boundaries are transparent to the outgoing waves, no reflection occurs at the boundaries.

The Reflective Boundary Conditions
The reflective boundary conditions are boundary conditions that are prescribed on a bounded interval [x l , x r ] and read For the classical NSWE system, we only need to impose the boundary conditions in one of the two dependent variables [29]. Therefore, it is sufficient to take u(x l , t) = u(x r , t) = 0 as reflective boundary condition.
In this paper, analogous reflective boundary conditions are considered. The incident wave η(x l , t) is coming from the left boundary, where the velocity can be expressed by where h(x, t) is the total depth of water and H(x l , t) = h 0 + η i (t) with η i (t) the free-surface elevation of the incident wave at the left boundary. This boundary condition was used by Stefanakis et al. [12]. A vertical wall is used at the right boundary, which means u(x r , t) = 0. Figure 3 shows an incident single wave generated at the left boundary that propagates to the right and meets a vertical wall at the right boundary, where it is reflected. When the reflected wave propagates back towards the left boundary, it reflects again and the wave polarity is reversed. It makes a multi-reflection system.

The Relaxation Zone
The relaxation zone is a kind of non-reflecting boundary layer. It absorbs all the reflected waves in a region that lies at the inlet boundary. It is usually one wavelength long [26,30]. In practice, it can modify an incident wave field to any desired target function. This also extends the applicability of the sponge layer to inlet boundaries, as the target function could be that of the generated wave train. The relaxation zone works as a blend between a target solution and the computed solution. The modified velocity can be described by where u target is a target velocity, u model is the unaltered computed velocity, and c relax is the relaxation coefficient. This technique has been discussed in many papers, such as Mayer et al. [31][32][33][34][35]. As the relaxation is applied explicitly as a postprocessing step, it is possible to apply it at the inlet. If there is no proper relaxation, the waves may propagate upstream toward the inlet boundary, be reflected back into the computational domain and pollute the flow field.
To obtain a smoothly modified internal flow field, it is required that "c relax " and its spatial derivatives are equal to zero along the interface connecting the relaxation zone and the internal field of the computational domain. There are many options to build a proper function of the relaxation coefficient. In this paper, we use the form where b ∈ [0, 2] is a free parameter, and β = 1 − (x − x start )/(x end − x start ) (x start and x end are coordinates of the relaxation zone boundaries) is the nondimensional distance function. The choice of the relaxation coefficient c relax depends strongly on the problem and the length of the relaxation zone. What is important is that the choice is not unique for a given problem.
In numerical simulations, reflected waves from the beach will interact with the wave generation boundary, which will pollute the final results, especially for breaking waves. Therefore, we can test the applicability of the relaxation zone through the modeling of standing waves. If it can absorb reflected waves in front of the inlet, the incident wave profile is kept in the relaxation zone and the standing wave pattern is simulated in the non-relaxed part of the computational zone.
The standing wave is modeled as sketched in Figure 4. The relaxation zone at the left boundary is one wavelength long. The total horizontal length of the domain is three wavelengths. The incident wave amplitude is η 0 = 1 m, the wave period is T = 40 s, and the water depth is h 0 = 20 m. One can see that the relaxation zone can simulate a standing wave pattern well in the non-relaxed part of the computational domain. The relaxation zone makes the boundary as "open". Compared with the ITBC, which can only be applied to regular waves, it can be wildly used for any type of wave, regular or irregular, long or short.

Results for Waves Incident on a Linearly Sloping Beach
Having described the three types of boundary conditions, we now present their effect on the numerical simulations of waves incident on a linearly sloping beach. Periodic monochromatic incident waves η(x l , t) = η 0 sin(ωt) propagate towards a sloping beach, with sloping angle tan α = 0.02 (unless otherwise specified) and horizontal length L = 5000 m (unless otherwise specified). The value tan α = 0.02 is typical of nearshore conditions as seen in [36].

Results with the Influx Transparent Boundary Conditions
The left ITBC is used on the left boundary. Non-dimensional run-up with respect to nondimensional wavelength is shown in Figure 5. One can see that the numerical nondimensional run-up matches the results given by the analytical linear solution before wave breaking. The ITBC is only approximate when the wavelength is too large as can be seen in Figure 5 for λ 0 /L > 3 [23,37,38]. The left boundary is not entirely transparent to outgoing waves when the wavelength is large. The outgoing waves are partially reflected and trapped in the computational domain. However, the effect on run-up is slight.

Results with the Reflective Boundary Conditions
We construct a multi-reflection system for the canonical case of Figure 1, where the reflective boundary conditions are used. Waves generated at the left seaward boundary propagate forward and are reflected at the right boundary. The reflected waves propagate back to the left boundary and are re-reflected, then the re-reflected waves propagate forward again and the process is repeated. A reflection cycle is defined as a cycle during which the wave travels the total domain back and forth. The time it takes is called the reflection period T re f lect . For the canonical case of Figure 1, it can be expressed by In the simulations below, the total simulation time is always four reflection cycles, and the incident periodic wave at x = −L t is expressed by where T is the wave period and η 0 = 0.1 m. First, we considered three combinations of α, L, and L t , thus fixing T re f lect for each combination. We tested the maximum run-up in a multi-reflection system for three cases by varying the nondimensional wave period within the range 0.25 < T/T re f lect < 0.65 (see Figure 6). Resonance effects are clearly visible. The maximum nondimensional run-ups are sharply amplified for certain wave periods, which lead to three peaks in Figure 6a, i.e., the maximum nondimensional run-up is up to 27 when tan α = 0.02, L t = 8000 m, L = 5000 m, and T = 220.8 s (T/T re f lect = 0.2658). Second, we fixed α, L, and T instead of L t , the horizontal distance from the left boundary to the undisturbed shoreline. This time, we tested the maximum run-up as a function of L t , which in turn made T re f lect and T/T re f lect vary through Equation (35). The results are shown in Figure 7. Waves with the same period propagating on the same bathymetry (tan α = 0.02, L = 2000 m) lead to different maximum run-ups with respect to various flat bottom lengths, which indicates that the wave behavior is affected by the location of the left boundary. This is different from the open boundary. However, resonance occurs when T/T re f lect = 0.2658, 0.3595 in the range of 0.25 < T/T re f lect < 0.4. From these two figures, monochromatic forcing is found to lead to resonant run-up by non-leading waves when the nondimensional wave periods T/T re f lect are the roots of the determinant of system (17), which coincides with the resonance condition we gave above.

Results with the Relaxation Zone
Next, we investigate the wave behavior in the presence of a relaxation zone. Its coefficients are given in the equation for c relax (34). However, we use the parameters b = 2 and β = (x end − x)/λ 0 , where λ 0 is the wavelength of the incoming wave. When the full relaxation zone is applied, the value of "c relax " decreases from 1 to 0 along the relaxation zone. When the relaxation zone is only partial, the value of "c relax " decreases from a value that depends on the length of the zone (value less than 1) to 0.
Results are shown in Figure 8. Two different relaxation zones are considered. To emphasize the effect of the relaxation zone, we also show the free-surface elevation time series with the use of the other boundary conditions described above (reflective and left ITBC). One can observe that the partial relaxation zone only alters the wave amplitude, but maintains the other wave characteristics. In other words, it can be considered as a partially reflective boundary condition. When the length of the relaxation zone is one fifth of the incident wavelength, the results coincide with those of obtained with the left ITBC, which means that there is no reflection caused at the left boundary.

Conditions for Run-Up Resonance in a Multi-Reflection System
To investigate the influence of the boundary conditions on the resonance mechanism, we first consider the reflective boundary conditions. In Appendix A, we describe a technique to decompose a wave into left-and right-going waves in a multi-reflection system. Using that technique, we separated the recorded waves at a point x located 2000 m to the right of the left boundary, into left-and right-going waves with L t = 7000, 8000, 9000 m. The other parameters were fixed at L = 2000 m, tan α = 0.02, T = 268.5 s, and η 0 = 0.1 m (see Figure 9). The amplitudes of left-and right-going waves are constant in each reflection cycle. When L t = 8000 m, the amplitude of right-going waves is the sum of amplitudes of the initial incident wave and of the left-going wave in the preceding reflection cycle, which are not observed in the other two cases (see Figure 9).
The wave motion could develop as a standing wave pattern with the use of open boundary condition [26]. For a multi-reflection system, it can be considered as a superposition of a standing wave and traveling waves. When the phase difference between the right-going and left-going waves is π, resonance can occur. When L t = 8000 m, one can observe that the phase difference between the right-going and left-going waves is close to π (see Figure 9b,c), which leads to resonance and amplifies the free surface elevation significantly (see Figure 9a). The period T = 268.5 s (T/T re f lect = 0.2658) used in Figure 9 is the resonant frequency that leads to the first peak in Figure 6a. Therefore, det(system) = 0 gives an analytical resonance condition, and Figure 9 reveals a wave motion resonance condition. We can conclude that the phase difference between the right-going and left-going waves must be π for the waves of the periods at the run-up peaks and det(system) = 0 in Figure 7. Due where m > 1. The average value is considered as the reflection coefficient of the left boundary. Tables 1 and 2 list reflection coefficients corresponding to cases of resonant and non-resonant regimes for the three types of boundary conditions. The amplitudes of the harmonic components being small, tiny errors may lead to deviations, and the reflection coefficients have small discrepancies even with the use of the same boundary condition. The reflection coefficient of the reflective boundary condition is up to 95%, it can be regarded as "full reflection". With the use of a partial relaxation zone (L relax = 0.2λ 0 ) and the left ITBC, amplitudes of the harmonic components are equal to zero when m > 1, which indicates that there is no reflection at the left. The left boundary is open.  In Figure 10, we tested the maximum run-up as a function of wave period and length of the partial relaxation zone, respectively. It can be observed that the maximum run-ups of the resonant regimes are determined by the length of the partial relaxation zone, which are sharply reduced by 40% when L relax = 0.1λ 0 and hold constant when L relax > 0.18λ 0 . It indicates that the relaxation zone does not require a full wavelength to behave as an open boundary.

Conclusions
Long wave run-up resonance in multi-reflection systems is studied in the paper. A simplified geometry involving long periodic waves incident on a linearly sloping beach is considered. By decomposing the surface wave into right-and left-going waves, we are able to obtain a better understanding of the conditions leading to resonance. It is found that wave trapping and wave-wave interactions can lead to run-up resonance. Resonance occurs when a condition obtained by linearizing the equations is satisfied. This condition is written as a determinant that must vanish. This determinant involves the angle of the slope, the length of the slope and the location of the artificial seaward boundary where boundary conditions are applied. The maximum run-up for the resonant regime is largely determined by the intensity of wave trapping.
Meanwhile, the fully and partially reflective boundary conditions are discussed. The reflection coefficient is determined by the length of the partial relaxation zone. The partial relaxation zone can also behave as an open boundary, even if it is less than one wavelength long.
The present study is mathematical. What it says is that the wrong choice of boundary conditions can lead to a misinterpretation of what happens in the real-world. However, there are several physical phenomena that can lead to wave trapping, and therefore to run-up resonance if the bathymetry has the right properties. In recent tsunamis, areas that were close to each other have experienced different run-ups, sometimes with an order of magnitude difference. The results in the present paper could provide a partial explanation.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript.
where A is the amplitude and κ is the phase difference. The subscripts I and R denote the incident and reflected waves, respectively. The second subscripts B and F denote the bounded and free waves, respectively. The superscript n denotes the nth harmonic, and the superscript m denotes the mth reflection cycle.
The waves are decomposed into the incident and reflected waves in terms of the harmonic bound and free waves. The bound waves are in phase-locked modes that travel at the same phase velocity as the first harmonic waves and are bound to them. The free waves propagate at their respective phase velocities, which can be presented by the dispersion relation (nω) 2 = gk (n) tanh k (n) h 0 . (A2) Let us verify the reconstruction of the wave motion provided by Equation (A1) on an example. Considering the canonical problem shown in Figure 1, together with the use of reflective boundary conditions, we reconstruct the recorded waves on the flat segment at x = −4000 m for a resonant and a non-resonant regimes (see Figure A1). A good agreement is found between the reconstructed and recorded waves. A new component is added in each reflection cycle. A jump appears in the time series of the free-surface elevation when it arrives (see Figure A1b). The jumps can vanish when we smooth the amplitude numerically. As the result shown in Figure A1 proves that Equation (A1) can describe the multi-reflection system, we consider that the smoothing process can be ignored.