Ecole Normale

Superieure De Cachan

Themistoklis S Stefanakis

Denys Dutykh

Francesco Carbone

Claudio Viotti

Ashkan Rafiee

Emile Contal

Emile Richard

Laura O'brien

Emiliano Renzi

David Buffoni

Roxana Tiron

Sarah Gallagher

Savvas Tzanis

Elena Siapani

Elli Androulaki

John Karlos

Gerasimos Skaltsas

Alexandros Toulias

Alexandros Gkiokas

Filippos Stavrinidis

Mara Egglezaki

Adda Koufopoulou

Fotis Bouzanis

Konstantina Iosif

Christos Stavrakis

Katerina Xviii

"She moves in mysterious ways" Bono

, who further distinguished the resonant run-up frequencies from the natural frequencies of the wavemaker. We also show that this resonant run-up amplification can be observed for several wave profiles such as bichromatic, polychromatic and cnoidal. However, resonant run-up amplification is not restricted to infinitely sloping beaches. We varied the bathymetric profile, and we saw that resonance is present in the case of piecewise linear and real bathymetries.

. The solution, which is based on separation of variables and a double Fourier transform in space, is exact, easy to implement and allows the study of realistic waveforms such as N-waves. We show the existence of focusing points for N-wave type initial displacements, where unexpectedly large wave amplitudes may be observed. N-wave focusing is found analytically by both linear dispersive and linear non-dispersive theories and is a possible explanation for run-up amplification beyond what is assumed by widely used scaling relationships.

I could not but start with my two PhD advisors, Professors Frédéric Dias and Nicolas Vayatis. I thank them for the trust they showed towards me, the opportunities they gave me and for the new horizons they opened not only in research but also more generally through our numerous discussions. They are both great mentors and enlightening personalities, each one in his own way. I know that I can always count on their guidance and support.

The PhD is just the end of a long path and the beginning of a new one. Throughout this path I had the continuous intellectual, emotional and financial support which was generously provided by my parents Stefanos and Agape. Together with my sister Eirini

In the third part I will explore the effect of localized bathymetric features on long wave generation. The generation of waves in a fluid domain of uniform depth by uplift or subsidence of a portion of the flat bottom boundary has been elegantly studied by [START_REF] Hammack | Tsunamis -A Model of Their Generation and Propagation[END_REF] for idealized motions. However, even when the final displacement is known from seismic analysis, the deforming seafloor includes relief features such as mounts and trenches. We investigate analytically the effect of bathymetry on the surface wave generation, by solving the forced linear shallow water equation. Our model for bathymetry consists of a cylindrical sill on a flat bottom, to help understand the effect of seamounts on tsunami generation. We derive the same solution by applying both the Laplace and the Fourier transforms in time. We find that as the sill height increases, partial wave trapping reduces the wave height in the far field, while amplifying it above the sill.

Finally, in the last part I will try to explore whether small islands can protect nearby coasts from tsunamis as it is widely believed by local communities. Recent findings for the 2010 Mentawai Islands tsunami show amplified run-up on coastal areas behind small islands [START_REF] Hill | The 2010 Mw 7.8 Mentawai earthquake: Very shallow source of a rare tsunami earthquake determined from tsunami field survey and near-field GPS data[END_REF], compared with the run-up on adjacent locations, not influenced by the presence of the islands. We will investigate the conditions for this run-up amplification by numerically solving the nonlinear shallow water equations. Our bathymetric setup consists of a conical island sitting on a flat bed in front of a plane beach and we send normally incident single waves.

The experimental setup is governed by five physical parameters. The objective is twofold: Find the maximum run-up amplification with the least number of simulations. Given that our input space is five-dimensional and a normal grid approach would be prohibitively computationally expensive, we present a recently developed active experimental design strategy, based on Gaussian Processes, which significantly reduces the computational cost. After running two hundred simulations, we find that in none of the cases considered the island did offer protection to the coastal area behind it. On the contrary, we have measured run-up amplification on the beach behind it compared to a lateral location on the beach, not directly affected by the presence of the island, which reached a maximum factor of 1.7.

Thus, small islands in the vicinity of the mainland will act as amplifiers of long wave severity at the region directly behind them and not as natural barriers as it was commonly believed so far.

Français

Cette thèse est divisée en quatre parties principales. Dans la première partie je vais présenter notre travail sur le run-up des vagues longues et sur les phénomènes d'amplification par résonance. Grâce à des simulations numériques basées sur les équations en eau peu profonde non-linéaires, nous montrons que dans le cas des vagues monochromatiques d'incidence normale sur une plage inclinée, une amplification résonante du run-up se produit lorsque la longueur de la vague d'entrée est 5.2 fois plus grande que la longueur de la plage [START_REF] Stefanakis | Local run-up amplification by resonant wave interactions[END_REF]. Ceci est cohérent avec les expériences en laboratoire de [START_REF] Ezersky | Physical simulation of resonant wave run-up on a beach[END_REF], qui distinguent bien les fréquences résonantes de run-up des fréquences naturelles du batteur.

Nous montrons également que cette amplification résonante de run-up peut être observée à partir de plusieurs profils de vagues (bichromatiques, irrégulières et cnoidales). Cependant, l'amplification résonante du run-up n'est pas limitée aux plages inclinées infinies. En faisant varier le profil bathymétrique, nous avons vu que la résonance est également présente dans le cas de bathymétries linéaires par morceaux et pour des bathymétries réalistes.

Dans la deuxième partie, je vais présenter une nouvelle solution analytique pour étudier la propagation des tsunamis générés par une source non ponctuelle sur une profondeur constante en utilisant la théorie des vagues en eau peu profonde linéaires [START_REF] Kânoğlu | Focusing of long waves with finite crest over constant depth[END_REF]. La solution, qui repose sur la séparation des variables et sur une double transformée de Fourier dans l'espace, est exacte, facile à mettre en oeuvre et permet l'étude d'ondes de formes réalistes comme les ondes en forme de N (N-waves). Nous montrons l'existence de points focalisants pour des déplacements initiaux de type N-wave, où des grandes amplitudes d'onde peuvent être observées de façon inattendue. La focalisation de vagues en forme de N est trouvée analytiquement à la fois par une théorie linéaire dispersive et par une théorie non-linéaire non-dispersive. Cette focalisation fournit une explication possible pour l'amplification du run-up au-delà de ce qui est communément admis par les relations d'échelle classiques.

Dans la troisième partie, je vais étudier l'effet de protubérances localisées sur la génération de vagues longues. La génération des vagues dans un domaine fluide de profondeur uniforme par soulèvement ou affaissement d'une partie du fond plat a été élégamment étudiée par [START_REF] Hammack | Tsunamis -A Model of Their Generation and Propagation[END_REF] pour des mouvements idéalisés.

Cependant, même lorsque le déplacement final est connu grâce à l'analyse sismique, le plancher océanique qui se déforme peut avoir du relief comme des montagnes et des failles. On étudie analytiquement l'effet de la bathymétrie sur la génération des vagues de surface, en résolvant les équations en eau peu profonde linéaires avec for. Notre modèle de bathymétrie se compose d'un rebord cylindrique sur un fond plat, afin de mieux comprendre l'effet des montagnes sous-marines sur la génération des tsunamis. Nous obtenons la même solution en appliquant à la fois les transformées de Laplace et de Fourier dans le temps. Nous constatons que quand la hauteur du rebord augmente, le piégeage partiel de la vague permet de réduire la hauteur des vagues dans le champ lointain, tout en l'amplifiant au-dessus du rebord. Je vais aussi présenter brièvement une solution de la même équation forcée au-dessus d'un cône.

Enfin, dans la dernière partie, je veux voir si les petites îles peuvent protéger les côtes proches de tsunamis comme il est largement admis par les communautés locales. Des découvertes récentes sur le tsunami des îles Mentawai en 2010 montrent un run-up amplifié sur les zones côtières derrière de petites îles [START_REF] Hill | The 2010 Mw 7.8 Mentawai earthquake: Very shallow source of a rare tsunami earthquake determined from tsunami field survey and near-field GPS data[END_REF], par rapport au run-up sur les lieux adjacents, qui ne sont pas influencés par la présence des îles. Nous allons étudier les conditions de cette amplification du run-up en résolvant numériquement les équations en eau peu profonde non-linaires. Notre configuration bathymétrique se compose d'une île conique sur un fond plat en face d'une plage inclinée. Nous envoyons des vagues solitaires avec incidence normale. Le dispositif expérimental est régi par cinq paramètres physiques. L'objectif est double: Trouver l'amplification maximale du run-up avec un nombre minimum de simulations. Etant donné que notre espace d'input est de dimension cinq et qu'une approche classique serait trop coûteuse en temps de calcul, nous présentons un plan d'expériences actif, récemment mis au point et basé sur les processus Gaussiens, qui réduit considérablement le coût de calcul. Après l'exécution de deux cents simulations, nous constatons que dans aucun des cas considérés l'île n'offre une protection à la zone côtière derrière elle. Au contraire, nous avons mesuré une amplification du run-up sur la plage derrière elle par rapport à une position latérale sur la plage non directement affectée par la présence de l'île. Cette amplification a atteint un facteur maximal de 1.7. Ainsi, les petites îles à proximité du territoire continental agissent comme des amplificateurs des vagues longues dans la région directement derrière elles et non comme des obstacles naturels comme il était communément admis jusqu'ici. The initial N -wave is located at (x * , y * ) = (27, 34)km and is tilted 5 • to be consistent with [START_REF] Synolakis | The slump origin of the 1998 Papua New Guinea tsunami[END_REF]. Sissano Lagoon, where most of the damage was observed, is located approximately r = 25km away from the initial wave location between the 30 • -45 • lines. Triangles over the 30 • and 45 • lines show the focusing points when p 0 = 5. (e) Maximum wave height envelopes for the source lengths L = 1 (solid line), 10 (dashed line), 20 (dash-dotted line) and 30 (dotted line). Triangles represent MOST numerical results.

List of Figures

(f ) Maximum wave height envelopes p 0 = 2 (dash-dotted line), 5 (dashed line) and 15 (solid line). When p 0 and L are parametrized, ε is modified to ensure the same maximum wave amplitude as the original initial wave profile. . . . . . . . . . . . . . . . . . . . . . . . 52 3.7 (a) Maximum wave envelope and time evolution at (b) t = 5 and (c) t = 30 for the PNG source, using linear non-dispersive (solid line), linear dispersive (dashed line), non-linear non-dispersive (triangles) and weakly non-linear weakly dispersive (squares) theories. Refer to the caption of Fig. 3.6 for the initial wave parameters. . . . . . . 53 3.8 (a) Maximum wave envelope and time evolution at (b) t = 12 and (c) t = 24 for the PNG source with the steepness parameter p 0 = 2 using linear non-dispersive (solid line), linear dispersive (dashed line), non-linear non-dispersive (triangles) and weakly non-linear weakly dispersive (squares) theories. Refer to the caption of Fig. The relevant region R is shown in light green (without edges). The first query x 0 is the maximizer of f + . We show in dashed line the upper and lower bounds with the update of σ after having selected x 0 . The second query x 1 is the one maximizing the uncertainty inside R. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Relationship between the simple regret r K t , unknown in a real situation, and the rank correlation ρ Xv (π t-1 , π t ) (in log-scale), for the synthetic function Himmelblau. The stopping threshold ρ 0 was set to 10 -4 , and the lag ℓ to 4, the algorithm stopped at iteration 12, after having found a good candidate at iteration

Introduction

General Characteristics of Tsunamis

Tsunami is a japanese word meaning "harbor wave". Despite its oriental etymological origins, the first historical reference to a tsunami is conjectured to be the wave created after the volcanic eruption of Thera in ancient Greece, around 1500 -1450 BC. Others support that tsunami history dates even further back to around 6100 BC, when a tsunami was supposedly triggered in the Norwegian Sea When tsunamis reach coastal waters, they increase their amplitude, while their propagation speed decreases.

Generation

The life of a tsunami is usually divided into three phases: the generation, the propagation and the inundation. For a tsunami to be generated, a large volume of water must be displaced. The most common generation mechanism is an earthquake-triggered uplift or subsidence of the oceanic floor. Convergent adjoint tectonic plates subduct the one below the other. At the interplate boundaries, strain is accumulated between earthquakes as the plates are locked together in their motion. During an earthquake, the plates break free and the overriding plate springs seaward in an upward motion, which displaces a large amount of water and thus creates a tsunami (Fig. 1.1).

In addition to earthquakes, landslides can also generate tsunamis and several recent

Propagation

Tsunami propagation is the phase which is best understood. Being a water wave problem, one should consider Laplace's equation

∇ 2 ϕ + ∂ 2 ϕ ∂z 2 = 0, (1.1)
together with the free surface (z = η(x, y, t)) and bottom boundary (z = -h(x, y))

conditions 1 KFSBC: ∂η ∂t + ∇ϕ • ∇η - ∂ϕ ∂z = 0 at z = η(x, y, t) , (1.2)
DFSBC:

∂ϕ ∂t + 1 2 |∇ϕ| 2 + 1 2
( ∂ϕ ∂z ) 2 + gη = 0 at z = η(x, y, t) , (1.3)
BBC: ∇ϕ • ∇h + ∂ϕ ∂z = 0 at z = -h(x, y) , (1.4)
where

∇ . = ( ∂ ∂x , ∂ ∂y )
is the horizontal gradient operator, ϕ is the velocity potential and h is the undisturbed water depth. Once the previous set of equations is combined with an initial condition, one can solve for ϕ and η. However, as stated earlier, tsunamis are long waves and thus their wavelength is much larger than the water depth, while their amplitude is much smaller than the water depth when traveling in the ocean. Hence, several approximations can be made to the full water wave problem. The most commonly used ones are the Boussinesq approximation, which retains nonlinearity and dispersion to a certain degree, and the shallow water equations, which are non-dispersive and can be either nonlinear or linear in the simplest case. In the rest of the thesis we will mostly consider the system of nonlinear shallow water equations (NSWE)

∂η ∂t + ∂ ∂x (u(h + η)) + ∂ ∂y (v(h + η)) = 0 , (1.5) ∂u ∂t + u ∂u ∂x + v ∂u ∂y + g ∂η ∂x = 0 , (1.6
)

∂v ∂t + u ∂v ∂x + v ∂v ∂y + g ∂η ∂y = 0 , (1.7)
where the Coriolis and viscous forces have been neglected. The above system of equations is used in various tsunami models because it is proved to capture well the main features of oceanic long waves. Among the several numerical NSWE models, we find the Method of Splitting Tsunami (MOST) [START_REF] Titov | Numerical modeling of tidal wave runup[END_REF] [2011b], which we will use in some of the chapters that follow. Only a few tsunami codes use the Boussinesq equations which take into account dispersion. One of the most known ones is FUNWAVE [START_REF] Kirby | FUNWAVE 1.0, Fully nonlinear Boussinesq wave model documentation and user's manual[END_REF] which was developed at the University of Delaware and is constantly upgraded.

Inundation

The last stage of a tsunami is the inundation or run-up phase. Wave run-up is defined as the maximum vertical extent of wave uprush on a beach above still water level [START_REF] Sorensen | Basic Coastal Engineering[END_REF]. During this phase, the tsunami decelerates and its amplitude grows significantly (Fig. 1.2), usually reaching 4 -15m, depending on the source and the local bathymetry. Tsunamis at this late phase are usually associated with strong currents which are responsible for infrastructural damages and life losses. Due to the high concentration of debris and the uncertainty that governs several parameters such as the drag coefficient of structures, the flow becomes very complex and thus it is difficult to simulate or study.

From a mathematical point of view, the difficulty lies in the moving shoreline.

If one considers a no-slip condition at the intersection of the free surface with the seafloor, then the shoreline will not move, which certainly is not the case in reality. Despite these contradictions, progress has been made, especially for the 1-D problem of a linearly varying bottom. The linear theory can predict the free surface elevation at the standing shoreline, which gives a good estimate of the run-up as it was shown by [START_REF] Synolakis | The runup of solitary waves[END_REF]. Nonlinear shallow water theory has also yielded results, this time on the moving shoreline by the use of the ingenious hodograph transformation applied by [START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF]. Later, several other studies followed, which used this transformation, but we will provide more information in Chapter 2.

Tsunami Engineering

The Boxing Day tsunami alerted the world about the fierceness and global reach of tsunamis, and initiated an effort to better understand their nature, develop warning systems and achieve higher levels of preparedness. Scientists and engineers worked together to pioneer ways to protect people and infrastructure from these destructive waves. Tsunami engineering mainly focuses on: (i) tsunami detection, warning and forecasting; (ii) regional tsunamigenic risk assessment and vulnerability; (iii) construction of efficient defenses such as sea walls and tsunami forests, as well as establishment of adequate design criteria for critical infrastructure and coastal structures under tsunami threat.

Detection, Warning and Forecasting

Due to the complexity and uncertainty as to whether an undersea earthquake has the potential to generate a tsunami, measurements of free surface elevation are necessary to determine whether a tsunami has been generated indeed, and thus reduce the risk of false alarms. To make these measurements, deep ocean tsunami detection buoys have been developed. These buoys comprise of two parts; a pressure sensor anchored to the sea floor and a surface buoy (Fig. 1.3). The sensor measures changes in pressure and associates them to changes in water height above it. This water column height is communicated to the surface buoy by acoustic telemetry and then relayed via satellite to the tsunami warning centre.

Tsunami warning centers, in addition to the buoy measurements, gather information on earthquake parameters such as epicenter location, depth, moment magnitude, as well as the focal mechanism (e.g. through the W-Phase centroid moment tensor by [START_REF] Kanamori | Source inversion of W phase: speeding up seismic tsunami warning[END_REF]). Tsunami warning centers have precomputed scenarios for their operational region, based on numerical simulations and historical data. A linear combination of these precomputed scenarios is used to best fit the new observational data, which gives estimates of the tsunami characteristics in deep water. These characteristics are used as initial conditions for site-specific inundation calculations. According to the calculations and to a decision matrix, the corresponding warning is issued.

Tsunami Hazard Assessment

Apart from real time warnings in case of a tsunami, governmental organizations and local communities need to be informed about the tsunami threat they face.

The evaluation of tsunamigenic risk at a specific area depends on the objective of its use and can be both deterministic and probabilistic. For instance, tsunami evacuation maps have been derived from inundation maps based on the maximum credible tsunami for a given region (i.e. deterministic scenarios) [START_REF] Geist | Probabilistic Analysis of Tsunami Hazards[END_REF]]. On the other hand, insurance applications focus on the annual probability of exceedance of a given flood depth. However, even the deterministic approach includes a notion of probability when one considers the possible events (e.g. tsunamis generated by volcanic eruptions or even submarine landslides are usually not considered) or the range of input source parameters. Therefore, the most critical point of the deterministic approach is the inclusion of a large amount of data, and thus a long time perspective of tsunami catalogues is needed to assess the risk of catastrophic events like the 1755 tsunami that destroyed Lisbon.

Paleo-tsunami science can help extend these catalogues and our understanding of recurrence intervals by studying paleo-tsunami deposits along the coasts in the sub-seafloor [e.g. [START_REF] Baptista | Revision of the Portuguese catalog of tsunamis[END_REF].

Probabilistic tsunami hazard assessment (PTHA) is inspired by the probabilistic seismic hazard assessment (PSHA) [START_REF] Cornell | Engineering seismic risk analysis[END_REF]. PTHA can be beneficial for region-wide analyses or in areas where there are few historical data. The PTHA is decomposed in three steps: (1) specification of the earthquake source parameters and associated uncertainties; (2) numerical propagation and run-up of the sources; and (3) the probabilistic calculations [START_REF] Geist | Probabilistic Analysis of Tsunami Hazards[END_REF]]. The first step (tsunami generation) is to determine the maximum tsunami amplitude (η max (r 0 , ψ s )) at a particular source location given by the location vector r 0 and a parameter space ψ s [START_REF] Geist | Probabilistic Analysis of Tsunami Hazards[END_REF]Parsons, 2006, Ward, 2001]. For earthquake sources, η max (r 0 , ψ s ) is taken to be equal to the coseismic vertical displacement field for offshore regions. In the second step, a numerical solver is used to propagate the initial wave to the coastal area of interest defined by the vector r. After running several simulations with different initial conditions, one can find the set of source parameters ψ crit s for which η max (r, ψ crit s ) > η crit . In the third step, we perform the statistical analysis. If the rate at which a given source occurs ṡ(r 0 , ψ s ) is known (following the Gutenberg-Richter law for example), then the total rate of this source that results in wave amplitudes exceeding η crit is [START_REF] Geist | Probabilistic Analysis of Tsunami Hazards[END_REF]]

Ṡ(r, r 0 , η crit ) = ∫ ψ max s ψ crit s ṡ(r 0 , ψ s ) dψ s . (1.8)
Then, by integrating over all possible sources, we can find the total number of tsunamis per year occurring at the location of study with amplitude greater than η crit , Ṡ(r, η crit ) = ∫ Ṡ(r, r 0 , η crit ) dr 0 . Another practical measure one can find is the probability a tsunami with amplitude grater than η crit occurs in a period T .

To do that one has to assume a distribution for the probability of occurrence of tsunamis in time.

Coastal Defenses

If a certain region identifies that it is exposed to a considerable tsunami threat, then the authorities may decide to take protective actions. Some of the better known protective works are the sea walls (Fig. 1.4) which were used for example in Japan (approximately 40% of the coastline is covered with sea walls), and the so-called tsunami forests [Danielsen et al., 2005, Harada and[START_REF] Harada | Effects of coastal forest on tsunami hazard mitigationa preliminary investigation[END_REF]].

Yet, the design criteria of these preventive works should be adequate and thus for tsunamis one needs to estimate the risk over very large time periods spanning several centuries or even millennia. Therefore, the study of paleotsunamis could enhance the evaluation of risk. Speaking about design criteria, tsunamigenic risk must also be incorporated in the design of critical infrastructures and coastal structures but still the tsunami-induced forces on these structures is poorly studied. The tremendous destruction in Japan clearly showed that some aspects of the tsunamis are not well understood and/or modeled. More specifically, one of the identified tsunami stages where there is high uncertainty is the description of the source. Many studies after each tsunami event are carried out to define the source which best matches the run-up observations. In some cases, to get a good fit, scientists make extra assumptions which are not clearly justified: For example the assumption that the earthquake initiated a submarine landslide. Then they try to find the size, the location and the motion of this landslide which best describes the observations. Even though it is true that an earthquake can ignite a landslide, the frequent lack of GPS data which would measure the characteristics of the landslide, if any, turns the assumption into a guess. Before we jump to the conclusion of a landslide, we first need to rule out the possibility that something more fundamental is not taken into account. For tsunami generation, apart from the earthquake rupture mechanism, two aspects that have already been identified

to play a role on the wave profile are the dynamic generation (i.e. the dependence of the seafloor motion on time) [START_REF] Hammack | Tsunamis -A Model of Their Generation and Propagation[END_REF][START_REF] Kervella | Comparison between three-dimensional linear and nonlinear tsunami generation models[END_REF], as well as the presence of sediments in the generation area [START_REF] Dutykh | Influence of sedimentary layering on tsunami generation[END_REF]. What is not studied yet, is the effect of seamounts or shelf breaks (i.e. extreme bathymetric features) inside the generation area and whether the passive generation is suitable in this case.

Tsunami propagation is the phase best understood. However, for long, tsunamis were described as solitary waves, whereas now it is clear that they are not [Tadepalli and Synolakis, 1994]. The initial free-surface deformation obtained by Okada's solution [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF] resembles a finite-crested N-shape profile, which inspired Tadepalli andSynolakis [1994, 1996] to introduce the notion of N-waves. While nowadays this model is widely accepted, the propagation of these waves is only scarcely studied, especially compared to the extensive investigation of the behavior of solitary waves.

One tantalizing question about tsunami propagation and run-up is why the folklore believes that the third wave is the largest and not the leading wave as scientists think so far. Is it due to reflections and bathymetric focusing from ridges and other bathymetric features which act a diffractive lenses [START_REF] Berry | Focused tsunami waves[END_REF] ? Is it due to dispersion or some type of resonance? Is it actually the third wave that is the largest? Could we predict which wave will be the largest at a given location?

Probably the situation is more complex and a simple answer cannot be given, but all previous questions are still open and can be addressed separately and in an idealized setting, which will give us some insight and will try to answer these questions in a more fundamental way.

Another widely spread belief among people living in coastal communities is that small islands in the vicinity of the shore offer protection to the regions behind them. However, simulations for the 2010 Mentawai tsunami [START_REF] Hill | The 2010 Mw 7.8 Mentawai earthquake: Very shallow source of a rare tsunami earthquake determined from tsunami field survey and near-field GPS data[END_REF] (Fig. 1.5) showed amplified run-up in these areas. Therefore it is under question whether this was an isolated observation or it is a more general phenomenon.

And if this is more general, which are the parameters that control it and which combination is the worst? The last question inherently introduces the notion of exploration of a parameter space and optimization at the same time. This is a much more general problem but in tsunami science it has been rarely applied, with the PTHA being the most obvious application. So far Monte Carlo methods dominate the scene in PTHA. However, when a large number of parameters are involved, typical in natural sciences, the computational cost rises prohibitively. Hence, it is interesting to introduce new methods which are more computationallyprudent from the exploding field of Machine Learning.

Summary of Contributions

Resonant run-up amplification of transient long waves

Motivated by the belief that it is not always the first tsunami wave the largest, we investigated numerically the run-up of long waves in the framework of the nonlinear shallow water theory. For long, the run-up of long waves has been treated as a standing wave problem, but the period until a progressive wave train becomes a standing wave due to reflection was disregarded so far. To study this, we considered the boundary value problem and initially we forced monochromatic waves on an infinite sloping beach.

Regardless of the beach slope and the length of the beach, we found that resonant run-up amplification could occur when the incoming wavelength is 5.2 times larger than the beach length [START_REF] Stefanakis | Local run-up amplification by resonant wave interactions[END_REF]. This corresponds to the first zero of the Bessel function J 0 (2ω √ L/gα), as it was first briefly noted by [START_REF] Billingham | Wave Motion[END_REF] in the context of linear shallow water theory, where ω is the cyclic frequency of the wave, L is the beach length, g is the gravitational acceleration and α is the beach slope. Resonant regimes can also be observed for subsequent zeros of the Bessel function J 0 , but with diminishing amplification factor. These results were confirmed experimentally by [START_REF] Ezersky | Physical simulation of resonant wave run-up on a beach[END_REF], who further distinguished this type of resonance from wavemaker resonance.

According to linear theory mentioned above, in resonance, the seaward boundary is a node of the fully developed standing wave. However, when the free surface is initially undisturbed, and thus the wave is transient, it takes a finite time for the node to travel to the seaward boundary. During that time, resonant run-up amplification can be significant.

Taking one step further, we tested the persistence of the resonant run-up amplification to modal perturbations and we found that both bichromatic and polychromatic waves can demonstrate such behavior when one of their frequencies is in the resonant regime. Moreover, resonant run-up amplification can also be observed for cnoidal waves, but the picture is more complicated.

Run-up resonance is not restricted to infinite sloping beaches. Piecewise linear, synthetic and real bathymetric profiles also exhibit resonant frequencies, which can either be evaluated analytically (piecewise linear) or numerically. In the particular case of a sloping beach connected to a constant depth region, resonance can only be observed when the incoming wavelength is larger than the distance from the undisturbed shoreline to the seaward boundary. Wavefront steepness is also found to play a role in wave run-up, with steeper waves reaching higher run-up values.

Given that tsunamis arrive at the shore as a series of waves, resonant run-up amplification may help explain why it is not always the first wave the most destructive.

Focusing of finite-crested N-waves over constant depth

The propagation of tsunamis from a finite strip source over constant depth is studied analytically, by means of linear shallow water theory. We derive a new solution based on separation of variables and a double Fourier transform in space [START_REF] Kânoğlu | Focusing of long waves with finite crest over constant depth[END_REF]. This solution is exact and allows the study of realistic waveforms, such as N-waves. We find that N-wave type initial displacements feature a focusing point where the wave amplitude is amplified beyond expected.

The focusing takes place in the leading depression side as two smaller waves emerge from the sides of the main wave and focus towards the center. The location of the focusing point and the wave amplification, which can be significant, depend on the geometric characteristics of the initial disturbance. N-wave focusing is studied analytically using linear non-dispersive and linear dispersive theory, and numerically using nonlinear non-dispersive and nonlinear weakly dispersive theory.

By comparing the previous theories, only small differences are found for the wave amplitude envelope and we can conclude that focusing is an intrinsic characteristic of N-waves.

Our results imply that when a shore lies on the radius of the focusing point, run-up on this area can potentially be significantly higher. Of course, the real bathymetry is never flat but this means that in certain circumstances, bathymetric focusing [START_REF] Berry | Focused tsunami waves[END_REF] can amplify run-up even more.

Long wave generation above a cylindrical sill

The physical problem of tsunami generation triggered by tectonic displacements involves substantial motions of the seafloor, which is not of constant depth. The typical practice is to ignore seafloor features, and simply assign an initial condition based on the elastic solutions which are calculable through seismic inversions, strictly applicable only for constant depth generation. When bathymetric gradients are small, then the flat seabed approximation appears reasonable. However, this is seldom the case; in real subduction zones, where tsunamis are often forced in regions with steep gradients and large features (Fig. 1.6), the question remains how large localized bathymetric features such as seamounts may affect the generated tsunami.

To study the influence of seamounts in the generation of tsunamis, we approach the problem analytically. We use the crude approximation of a seamount as a submerged cylindrical sill and we solve the forced linear shallow water equations (LSWE) in cylindrical coordinates. We divide the physical space in the near field, which is the area above the sill and the far field, which is the area beyond the sill. To solve the equations we can choose between applying the Laplace transform in time, which is the integral transform of choice for transient phenomena, or the Fourier transform. In the transformed space we solve the forced LSWE in each subregion and we match the solutions at the common boundary by ensuring continuity of the free surface elevation and continuity of the radial fluxes. Going back to physical space requires the inverse transform. Both integral transforms yield the same result, but the Fourier transform is more straightforward because it does not involve the sometimes cumbersome contour integral evaluation.

Our results correlate well with those obtained by [START_REF] Hammack | Tsunamis -A Model of Their Generation and Propagation[END_REF] in the absence of the sill, even though our solution is non dispersive. Dispersive effects become important as the bottom displacement tends to instantaneous. The presence of the sill is found to reduce the amplitude of the wave in the far field, which even though intuitive since less volume of water is displaced has not been proved. In addition, due to partial wave trapping above the sill, smaller trailing waves are observed in the far field, the amplitudes of which grow with increasing sill height. During wave generation, above the sill, the maximum wave height increases with increasing sill height but is always bounded by the maximum bottom displacement.

On the other hand, above the sill, the wave is partially trapped and it becomes more apparent with increasing sill height. As the generated wave propagates towards the edge of the sill, it gets partially reflected back and directional focusing is observed above the center of the sill. As stated in the previous section (Section 1.3), there is a widespread belief that small islands in the vicinity of the mainland will protect the area behind them in the unfortunate event of a tsunami. This misconception contradicts the recent post-tsunami survey data and the simulations from the Mentawai 2010 tsunami.

In order to shed more light, we need to check whether the Mentawai tsunami was a particular case, or a more general amplification mechanism.

In this part, we study the run-up amplification of a long wave on a plane beach behind a conical island, with respect to the run-up on an adjacent location on the beach, not directly affected by the presence of the island. The geometrical setup is defined by five physical parameters, namely, the island slope, the beach slope, the distance between the island and the beach, the water depth and the cyclic frequency of the wave. The long wave is modeled by a single wave of fixed amplitude.

We treat this model as a nonconvex optimization procedure, where we seek to find the maximum run-up amplification with the least computational cost. This is a more general problem and thus can be applied to a wide spectrum of disciplines. A set of points in the parameter space is pre-selected by a classic (static) experimental design technique, the Latin Hypercube Sampling [START_REF] Mckay | Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF], which aims to better fill the parameter space. Initially we run a first batch of randomly selected simulations. Based on the observations obtained, we construct an emulator, or statistical model, hypothesizing that our problem can be approximated by a Gaussian Process. Then, we build an algorithm for the optimization problem which is parallelized and therefore batches of simulations can be ran simultaneously. This algorithm guides the selection of the query points in the parameter space.

After each batch, the emulator is updated based on the new observations, which accounts for its active nature. Our algorithm also takes into account the exploration/exploitation tradeoff which is essential to avoid getting "stuck" in a local maximum. Finally, we present a stopping criterion, independent of the algorithm, that serves as a signal for the end of the optimization, before the query of all the points in the input space. We test our algorithm versus other commonly used strategies not only on the database produced by the tsunami experiment, but also on synthetic databases. We find that it outperforms the other strategies in most cases.

In terms of physical results, we find that in all cases considered, the island focuses the wave behind it, which results in amplified run-up in the area behind it. The maximum amplification factor is ≈ 1.7 and the median amplification factor is 1.3.

Our algorithm reduces the computational cost by more than 60% compared to a classic experimental design and it is orders of magnitude more cost efficient than a regular grids approach -3 orders of magnitude for our 5-dimensional problem.

Perspectives

The present thesis sheds light to some aspects of tsunamis not well understood so far, but other questions still remain. To begin with, tsunami induced currents and the forces they exert on coastal structures and buildings, not addressed here, is a topic of particular interest due to its obvious consequences.

Another issue related to resonant run-up (Chapter 2), is the development of a boundary condition which will be more transparent on the treatment of the reflected wave component. The recent research by [START_REF] Antuono | Solving the nonlinear shallow-water equations in physical space[END_REF] may lead towards this end. Moreover, run-up resonance for obliquely incident waves was not studied, neither was studied whether this phenomenon would manifest for dispersive N-waves, the prevailing model for tsunamis.

Regarding the effect of bathymetry on tsunami generation (Chapter 4), further research is needed to investigate whether wave trapping can be observed above more general forms of seamounts, other than sills. Given the complexity of the bathymetry, classical analytical solutions might be limiting, but other formulations such as conformal mappings may give new insight in a more realistic setting.

Finally, the integration of statistical emulators in tsunami science can yield important improvements in fast warnings at coastal areas close to the generation zone.

In addition, they can be used for computationally efficient optimization problems (Chapter 5), but a more robust stopping criterion should be developed, which would either be derived from the optimization algorithm, or it will dependent on the dimensionality of the problem.

Chapter 2

Run-up amplification of transient long waves

Long wave run-up, the maximum elevation of wave uprush on a beach above still water level, is difficult to observe in nature in real time due to the large physical dimensions of the phenomenon and to the catastrophic consequences it usually leads to, since the most famous representation of a long wave is that of a tsunami.

Tidal waves, meteotsunamis [START_REF] Rabinovich | Meteorological tsunamis: Atmospherically induced destructive ocean waves in the tsunami frequency band[END_REF] and storm surges are also long waves.

Most observational data concerning run-up are collected during post-tsunami surveys. Nevertheless, this data does not offer any information, by itself, on the time history of the event, which leads field scientists to rely on interviews with eye-witnesses, who, in some cases, have reported that it is not always the first wave which results in the worst damage. Moreover, unexpected extreme localized run-up values have been measured during several tsunami events, such as in Java Stefanakis et al. [2011] showed that for a given plane beach slope there exist wave frequencies that lead to resonant long wave run-up amplification by non-leading waves. These results were confirmed experimentally in a wave tank by Ezersky et al.

[2013] who distinguished the frequency that leads to resonant run-up from the resonant frequency of the wavemaker. They also observed a secondary resonant regime which was not identified before. The authors also recognized that the resonant state occurs when the Bessel function J 0 ( √ 4ω 2 L/g tan θ) = 0, as predicted by the linear theory [START_REF] Lamb | Hydrodynamics[END_REF], where ω is the angular frequency of the wave, tan θ is the beach slope, L is the horizontal distance from the undisturbed shoreline to the point where the wave amplitude is imposed and g is the acceleration due to gravity. Several other possible explanations for the observed extreme run-up values are also available.

Miles [1971] described the conditions for harbor resonance and the importance of the Helmholtz mode to tsunami response and later [START_REF] Kajiura | Local behaviour of tsunamis[END_REF] A thorough review of the long wave run-up problem with additional results on its relation to the surf-similarity is given by [START_REF] Madsen | Run-up of tsunamis and long waves in terms of surf-similarity[END_REF]. The aforementioned theoretical results do not exhibit any resonant regimes and were reproduced numerically by [START_REF] Madsen | Run-up of tsunamis and long waves in terms of surf-similarity[END_REF] by placing a relaxation zone close to the wave generation region, which absorbs the reflected wavefield.

These sponge layers are widely used because the combination of incoming and outgoing waves at the boundary still remains poorly understood. These are artifacts and are not part of the governing wave equations.

In the present chapter we first provide an overview of the theory behind long wave run-up on a plane beach1 and we confirm the resonance results of Stefanakis et al.

[2011] with more geophysically relevant bottom slopes. We also prove their robustness to modal perturbations. The case of a piecewise linear beach follows where we show both analytically and computationally the existence of resonant states. Then we explore whether resonance can be observed when a sloping beach is connected to a constant depth region and we test the effect of wave nonlinearity and how the results relate to the theory. Finally, we discuss the effect the boundary condition has on the resonant run-up amplification.

Statement of problem and method of analytical solution

In the following we present a review of the analytical solution. Consider a wave propagation problem described by the one-dimensional NSWE

∂η * ∂t * + ∂ ∂x * [(h + η * ) u * ] = 0, ∂u * ∂t * + u * ∂u * ∂x * + g ∂η * ∂x * = 0 (2.1)
where z * = η * (x * , t * ) is the free surface elevation, h(x * ) is the water depth, u * (x * , t * ) is the depth-averaged horizontal velocity and g is the acceleration due to gravity. Consider a topography consisting of a sloping beach with unperturbed water depth varying linearly with the horizontal coordinate, h(x * ) = -αx * , where α = tan θ is the bottom slope (see Fig. 2.1).

In order to solve equations (2.1), appropriate initial and boundary conditions must be supplied. In most wave problems, one must provide the initial conditions η * (x * , 0) and u * (x * , 0) (for tsunamis, it is usually assumed that u * (x * , 0) = 0).

The boundary condition far from the tsunami source area ("left boundary") is If the tsunami source is far from the shore, it is convenient not to include the source area in the fluid domain and apply the following "left" (incoming wave) boundary condition at some point x * = x * 0 :

η * (x * , t * ) → 0, u * (x * , t * ) → 0 (x * → -∞). (2.2) z * x * θ u * (x * , t * ) η * (x * , t * )
u * (x * 0 , t) = √ g/h(x * 0 )η * (x * 0 , t), (2.3)
which corresponds to the tsunami wave approaching the shore. Another boundary condition is the boundedness of all functions on the unknown moving boundary,

h(x * ) + η * (x * , t * ) = 0 , (2.4)
which determines the location of the moving shoreline. The condition (2.4) is the main difference from the classical formulations of the Cauchy problem for hyperbolic systems.

There is an analytical method for solving this system, based on the use of the Riemann invariants. These invariants for a plane beach are

I ± = u * ± 2 √ g(-αx * + η * ) + gαt * , (2.5)
and the system (2.1) can be rewritten as

∂I ± ∂t * + ( 3
4 I ± + 1 4 I ∓ -gαt * ) ∂I ± ∂x * = 0. (2.6)
It is important to mention that this approach is applied for water waves on a beach of constant slope, and there are no rigorous results for arbitrary depth profiles h(x * ). The existence of Riemann invariants in the general case is an open mathematical problem.

Then the hodograph transformation can be applied to the system (2.6), assuming that the determinant of the Jacobian J = ∂ (x * , t * ) / ∂ (I + , I -) does not vanish (this determinant vanishes when the wave breaks; we note that [START_REF] Synolakis | The runup of solitary waves[END_REF] has argued that this point is simply where the hodograph transformation becomes singular and the interpretation is that the wave breaks, and in fact corresponds to breaking during the rundown, at least for solitary waves). As a result, the following set of equations is derived:

∂x * ∂I ∓ - ( 3
4 I ± + 1 4 I ∓ -gαt * ) ∂t * ∂I ∓ = 0. (2.7)
These equations are still nonlinear but they can be reduced to a linear equation by eliminating x * (I + , I -) :

∂ 2 t * ∂I + ∂I - + 3 2 (I + -I -) ( ∂t * ∂I - - ∂t * ∂I + ) = 0. (2.8)
It is convenient to introduce the new variables

λ = 1 2 (I + + I -) = u * + gαt * , σ = 1 2 (I + -I -) = 2 √ g(-αx * + η * ) . (2.9)
Then the system (2.8) takes the form

σ ( ∂ 2 t * ∂λ 2 - ∂ 2 t * ∂σ 2 ) -3 ∂t * ∂σ = 0. (2.10)
Expressing the time t * from Eq. (2.9),

t * = λ -u * gα , (2.11)
and substituting

u * = 1 σ ∂Φ ∂σ , (2.12)
where Φ is a dependent variable named 'potential' in [START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF],

we finally rewrite Eq. (2.10) in the form of the classical cylindrical wave equation

∂ 2 Φ ∂λ 2 - ∂ 2 Φ ∂σ 2 - 1 σ ∂Φ ∂σ = 0. (2.13)
All physical variables can be expressed through the function Φ(σ, λ). In addition to the time t * (2.11) and the velocity u * (2.12), the horizontal coordinate x * and the water displacement η * are given by

x * = 1 2gα ( ∂Φ ∂λ -u * 2 - σ 2 2 ) , (2.14)
η * = 1 2g ( ∂Φ ∂λ -u * 2
) .

(2.15) So, the initial set of nonlinear shallow water equations has been reduced to the linear wave equation (2.13) and all physical variables can be found via Φ using simple operations. The main advantage of this form of the nonlinear shallowwater system is that the moving shoreline corresponds to σ = 0 (since the total depth h(x * ) + η * (x * , t * ) = 0 ) and therefore Eq. (2.13) is solved in the half-space Meanwhile, the typical situation in tsunamis is that the wave approaches the shore from deep water where the wave can be considered as linear. In this case it is possible to find the function Φ without using the implicit formulas of the hodograph transformation. Let us consider the linear version of the shallow water system:

0 ≤ σ < ∞
∂u * ∂t * + g ∂η * ∂x * = 0, ∂η * ∂t * + ∂ ∂x * (-αx * u * ) = 0, (2.16)
and apply the linearized version of the hodograph transformation

η * = 1 2g ∂Φ l ∂λ l , u * = 1 σ l ∂Φ l ∂σ l , x * = - σ 2 l 4gα , t * = λ l gα , (2.17)
where the subscript l denotes quantities derived from linear theory. In this case the system (2.16) reduces naturally to the same linear cylindrical wave equation

∂ 2 Φ l ∂λ 2 l - ∂ 2 Φ l ∂σ 2 l - 1 σ l ∂Φ l ∂σ l = 0, (2.18)
which has the same form as in nonlinear theory (2.13). If the initial conditions for the wave field are determined far from the shoreline, where the wave is linear, then the initial conditions for both equations (2.13) and (2.18) are the same, and therefore, their solutions will be the same,

Φ(σ, λ) = Φ l (σ l , λ l ), (2.19)
after replacing the arguments. So the function Φ can be found from linear theory.

From the operational point of view, it is important to know the extreme runup characteristics like run-up height, rundown amplitude, onshore and offshore velocity, and these characteristics can be calculated within the framework of linear theory. This surprising result, also noted by [START_REF] Synolakis | The runup of solitary waves[END_REF], can be explained as follows. Indeed, it follows from Eq. (2.19) that extreme values of Φ and its derivatives are the same. But for a moving shoreline (σ = 0) in extreme points of run-up or rundown, the velocity is zero, and the expressions of the hodograph transformations (2.15) and (2.17) coincide. So, it is believed that the extreme characteristics of tsunami run-up which determine the flooding zone can be found from linear theory despite the real nonlinear character of the wave process in the nearshore area, and this is an important result for tsunami engineering.

Moreover, the nonlinear dynamics of the moving shoreline (σ = 0) can be easily derived using linear theory. It follows from (2.11) that the moving shoreline velocity is

u * (λ, σ = 0) = λ -gαt * , (2.20)
or in equivalent form

u * (t * ) = u * l (t * + u * /gα), (2.21)
where the function u * l (λ) is found using the known function Φ. Therefore it can be found from linear theory (it is the velocity at the point x * = 0). Similarly, the water displacement should be found first from linear theory (at the point x * = 0)

z * l (t * ) = η * l (x * = 0, t * ) = α ∫ u * l (t * )dt * , (2.22)
and then one can find the "real" nonlinear vertical displacement of the moving shoreline,

z * (t * ) = η * (σ = 0) = α ∫ u * (t * )dt * = z * l (t * + u * /αg) -u * 2 (t * )/2g. (2.23)
As can be seen from these formulas, the extreme values of functions in linear and nonlinear theories coincide as we pointed it out already. The manifestation of nonlinearity is in the shape of the water oscillations on shore due to the nonlinear transformation (2.21).

z * x * θ u * (x * , t * ) η * (x * , t * ) x * = -L
η * (x * , t * ) = η R J 0 ( √ 4ω 2 |x * | gα ) cos ωt * , (2.24)
with η R an arbitrary constant. Using asymptotic expressions for the Bessel function J 0 and matching with the solution of the mild slope equation [see [START_REF] Madsen | Run-up of tsunamis and long waves in terms of surf-similarity[END_REF] one finds that the wave field far from the shoreline consists of the linear superposition of two waves propagating in opposite directions and having the same amplitudes (a standing wave):

η * (x * , t * ) = 2η 0 ( L |x * | ) 1/4 cos ( 2ω √ |x * | gα + φ ) cos ωt * , (2.25)
where the incident wave amplitude has been fixed to η 0 at the coordinate x * = -L.

The coefficient of wave amplification in the run-up stage is found to be

η R η 0 = 2 ( π 2 ω 2 L gα ) 1/4 = 2π √ 2L λ 0 , (2.26)
where λ 0 = 2π √ gαL/ω is the wavelength of the incident wave.

A more realistic example

The rigorous theory described above is valid for the waves in a wedge of constant slope. For all other depth profiles rigorous analytical results are absent. Real bathymetries, which are complex in the ocean, can be approximated by a beach of constant slope in the vicinity of the shore only. If the "matching" point is relatively far from the shoreline, the linear theory of shallow water can be applied for waves in a basin of complex bathymetry except in the nearshore area. Within this approximation, and arguing as in [START_REF] Synolakis | The runup of solitary waves[END_REF], the 1D linear wave equation

∂ 2 η * ∂t * 2 - ∂ ∂x * ( c 2 (x * ) ∂η * ∂x * ) = 0, c 2 (x * ) = gh(x * ) (2.27)
should be solved analytically or numerically, and then its solution should be matched with the rigorous solution of the run-up problem described above. A popular example of such matching is given for the geometry presented in Fig. 2.2, which is often realized in laboratory experiments. The elementary solution of the wave equation (2.27) for a basin of constant depth h 0 is the superposition of incident and reflected waves

η * (x * , t * ) = η 0 exp[iω(t * -x * /c)] + A r exp[iω(t * + x * /c)] + c.c., c = √ gh 0 (2.28)
with η 0 real and A r complex, and this solution should be matched with (2.24) at the point x * = -L using the continuity of η * (x * ) and dη * /dx * . As a result, the unknown constants A r and η R can be calculated from the boundary conditions at

x * = -L, and the run-up amplitude is

η R η 0 = 2 √ J 2 0 (χ) + J 2 1 (χ) , χ = 2ωL c = 4π L λ 0 . ( 2

Numerical results

The solutions described in the previous section are standing waves. If the motion starts from scratch, one does not have a standing wave at the beginning. A (2.24) with it at x * = -L yields

η * (x * , t * ) = 2η 0 J 0 ( √ 4ω 2 |x * | gα ) /J 0 ( √ 4ω 2 L gα ) cos ωt * . (2.30)
There is indeed the possibility of a resonance, which occurs when 2ω √ L/gα is a zero of J 0 and the solution (2.30) is then unbounded. Another way to look at this resonance is to consider Figure 9(a) in [START_REF] Madsen | Run-up of tsunamis and long waves in terms of surf-similarity[END_REF]. The resonance occurs when ones tries to force the wave amplitude to a finite value at one of the nodes of the solution (2.24). We will see below that this resonance can occur in numerical as well as laboratory experiments.

Waves on a plane beach

We present some results on the resonant long wave run-up phenomenon on a plane resonant run-up by non-leading waves (Fig. 2.4a) when the nondimensional wavelength is λ 0 /L ≈ 5.2, where λ 0 = 2π √ gαL/ω is the incident wavelength, which comes as a direct consequence of Eq. (2.30) when the wave at the seaward boundary is specified, since resonant states can be identified with the roots of the Bessel function J 0 . Resonance is also observed for the maximum horizontal velocities assumed by the waves (Fig. 2.4b), which are found at the shoreline. Since and

η * (x * , t * ) = η R J 0 (σ) cos ωt * , u * (x * , t * ) = - 2ωη R σ tan θ J 1 (σ)
lim σ→0 J 1 (σ) σ = 1 2 , (2.33)
then the shoreline velocity is

u * s (t * ) = - ωη R tan θ sin ωt * , (2.34)
which implies that when the run-up is resonant, so is the shoreline velocity. The resonance mechanism was found to rely on a synchronization between incident and receding waves, but should be distinguished from wavemaker resonance since the computational domain is not closed (Fig. 2.5), as it would be in a laboratory setting, and we can observe strong inflow-outflow during run-up and run-down.

Furthermore, the experiments of Ezersky et al.

[2013] confirmed this claim, by observing that the resonant frequency of the system is different from the frequency that leads to the resonant run-up.

The spatio-temporal behavior of the non-dimensional horizontal velocity is shown in Fig. 2.6. For visualization purposes, we plot only the last 500 m of the beach to the left of the initial shoreline position. We can observe that in the resonant regime, after the rundown induced by the leading wave, during run-up of subsequent waves, a fixed spatial point undergoes an abrupt change of velocity, from highly negative to highly positive values. Furthermore, the maximum absolute velocity increases over time in the resonant regime, while this is not true in the non-resonant case.

To give a feeling of dimensions, imagine a plane beach with tan θ = 0.01 and an incoming wave of amplitude η 0 = 1 m at L = 10, 000 m offshore where the water depth is h 0 = 100 m. For that wave to be in the resonant regime, its wavelength has to be approximately 52, 000 m. If there is a run-up amplification close to 40, this means that max u s ≈ 15 m/s according to Eq. (2.34).

In order to increase our confidence in the numerical solver that we use, we ran In order to further investigate the effects of modal interactions in the resonant regime, we tested incoming waves of bichromatic modal structure. To remain consistent with the monochromatic case, each mode had half the amplitude of the equivalent monochromatic wave (η 0 ). Our computations (Fig. 2.8) show that no important new interactions occur. When one of the two frequencies is resonant, the run-up is dominated by it, while the other does not alter the dynamics. If both frequencies are resonant, their constructive interference is small overall and does not differ significantly from the equivalent monochromatic resonant state.

Therefore, this result indicates that the resonant run-up mechanism is robust and is not restricted to monochromatic waves only.

In order to further investigate the robustness of the resonant run-up mechanism, we introduced 10 semi-random perturbation components to the monochromatic wave signal. The amplitude of the perturbations followed a normal distribution with zero mean and standard deviation much less than the wave amplitude η 0 . The perturbation frequency followed the lognormal distribution. The monochromatic wave and a corresponding semi-randomly perturbed signal in physical and Fourier spaces are shown in Fig. 2.9. By running the simulation in the resonant regime when the slope α = 0.02 and L = 5000 m we obtained the same run-up timeseries that we would obtain with the unperturbed monochromatic wave (Fig. 2.10).

ω 1 (L/g tan θ) 1/2 ω 2 (L/g tan θ) 1/2
Therefore we can conclude that the resonant run-up mechanism is robust and the resonant frequency dominates the run-up.

Piecewise linear bathymetry

Kânoğlu and Synolakis [1998] developed a general methodology to study the problem of long wave run-up over a piecewise linear bathymetry and applied it to the study of solitary wave run-up, but in their formulation, the last offshore segment of the bathymetry consisted of a flat bottom. Here we will only use two uniformly sloping regions as in Fig. 2.11. Following the steps of [START_REF] Lamb | Hydrodynamics[END_REF], we take the linearized form of the NSWE (2.1) and search for solutions of the type (2.37)

η * (x * , t * ) = Z(x * ) cos ωt * , u * (x * , t * ) = V (x * ) sin ωt * . ( 2
Since the bathymetry is piecewise linear, at each segment we have h i (x * ) = -α i x * + c i , where α i = tan θ i ̸ = 0, c i is a constant and the subscript i is indicative of the segment number from the shoreline to the seaward boundary. In that case, Eq.

(2.36) becomes a standard Bessel equation of order zero and the general solution for Z i and V i is

Z i (x * ) = A i J 0 (σ) + B i Y 0 (σ), V i (x * ) = 2ω σα i dZ i dσ , (2.38)
where A i and B i are linear coefficients, J n and Y n are the nth order Bessel functions of the first and second kind respectively and

σ = 2ω √ g √ -x * + c i α i α i . (2.39)
In order to solve this problem we require continuity of the free surface elevation and of the horizontal fluxes at two adjacent segments and we prescribe a wave amplitude Z i = η 0 cos ωt * at the seaward boundary. Here for simplicity we will focus on the case of two segments but it can be generalized to an arbitrary number of segments as shown in Kânoğlu and Synolakis [1998]. In the first segment boundedness of the free surface elevation at the shoreline (x * = 0) requires B 1 = 0.

Since J 0 (0) = 1, A 1 represents the run-up and therefore we will name it η R . Hence, we have the following linear system of equations:

J 0 (σ 1 )η R -J 0 (σ 2 )A 2 -Y 0 (σ 2 )B 2 = 0 (2.40) J 1 (σ 1 )η R -J 1 (σ 2 )A 2 -Y 1 (σ 2 )B 2 = 0 (2.41) J 0 (σ 3 )A 2 + Y 0 (σ 3 )B 2 = η 0 (2.42)
where

σ 1 = 2ω √ L 1 gα 1 (2.43) σ 2 = 2ω α 2 √ α 1 L 1 g (2.
44)

σ 3 = 2ω α 2 √ α 2 (L 2 -L 1 ) + α 1 L 1 g (2.45) Then η R = 0 -J 0 (σ 2 ) -Y 0 (σ 2 ) 0 -J 1 (σ 2 ) -Y 1 (σ 2 ) η 0 J 0 (σ 3 ) Y 0 (σ 3 ) / J 0 (σ 1 ) -J 0 (σ 2 ) -Y 0 (σ 2 ) J 1 (σ 1 ) -J 1 (σ 2 ) -Y 1 (σ 2 ) 0 J 0 (σ 3 ) Y 0 (σ 3 ) . (2.46)
Therefore, when the determinant in the denominator vanishes, the run-up becomes resonant (Fig. 2.12). Furthermore, the shoreline velocity is now given by Eq.

(2.38) as

V s = lim σ→0 - 2ωη R σα 1 J 1 (σ) ⇒ V s = - ω α 1 η R (2.47)
which indicates that when the run-up is resonant, so is the shoreline velocity. 2 The same argument of course applies to the case of the infinite sloping beach (see Fig.

2.4b

). Numerical simulations performed in this setting with α 1 = 0.02, α 2 = 0.01, L 1 = 5000 m and L 2 = 6000 m (η 0 = 0.1 m) agree with the above analytical solution and again resonant wavelengths can be identified (Fig. 2.13).

Plane beach connected to a flat bottom

A more characteristic bathymetric profile consists of a constant depth region connected to a sloping beach, hereafter referred as the canonical case (Fig. 2.2). Using this profile [START_REF] Madsen | Run-up of tsunamis and long waves in terms of surf-similarity[END_REF] showed very good agreement between theory and their computations for a range of wavelengths 1 < λ 0 /L < 7 and wave nonlinearity 0.001 < η 0 /h 0 < 0.01, in which even a tsunami at h 0 = 100 m is fairly linear. 1998] but it is somewhat artificial. The length of the relaxation zone should be comparable to the wavelength, but the resonant wavelength is found to be greater than the beach length, which is the reason why we could not employ it in the infinite slope case.

For the current bathymetric profile, our objectives were to investigate both if resonance would occur and whether the existence of the relaxation zone would play any role on the run-up. Hence we examined the run-up of monochromatic waves of amplitude η 0 = 1.25 m on a plane beach with slope tan θ = 0.02, which reached a maximum depth h 0 = 100 m. We performed simulations with and without the relaxation zone. The setup without relaxation zone is more natural since no artificial filtering is used but at the same time is more computationally demanding, due to the increased length of the constant depth region. In Fig. 2.14 we observe that both with and without the relaxation zone, the computations predict slightly higher maximum run-up values than the ones predicted by the theory in the non-breaking regime, but the qualitative behavior is the same.

The discrepancies between theory and computations are higher when the use of a relaxation zone is avoided. Like Pelinovsky and Mazova [1992], one can introduce the breaking number Br = ω 2 η R /gα 2 . When Br = 1, or η R /η 0 = g tan 2 θ/(η 0 ω 2 ), the analytical solution breaks down. When Br > 1, the wave breaks. According to [START_REF] Mei | Theory and applications of ocean surface waves[END_REF], this criterion can only be used as a qualitative criterion.

When waves are close to breaking, the run-up amplification reaches its maximum.

However, we cannot observe any significant resonance as we did in the infinite It is well known that in the context of NSWE, as waves propagate over a flat bottom, the wavefront tends to become steeper and the higher the wave nonlinearity, the faster the wave steepening. [START_REF] Synolakis | Evolution of maximum amplitude of solitary waves on plane beaches[END_REF] have shown that offshore and far from breaking the wave evolves with Green's law, while closer to breaking the evolution is more rapid, and they named this region the Boussinesq regime. In the previous case the discrepancies observed with and without the relaxation zone could be attributed to the different lengths of the constant bottom region (hereafter L 0 ). Before, the wave nonlinearity was η 0 /h = 0.0125 and in order to increase the effect of wave steepening we decided to double the incoming wave amplitude. Hence, we tested three different cases, namely the same two as before, one without a relaxation zone and L 0 = 4λ 0 , one with a relaxation zone and L 0 = 2λ 0 and finally one with a relaxation zone but now the constant depth region has a length equal to 4 wavelengths.

In Fig. 2.15 we see that the influence of L 0 is important and hence the wavefront steepness is critical to the run-up amplification. The existence of the relaxation zone does not affect the results when the constant depth region has a fixed length.

The longer L 0 and therefore the wavefront steepness, the higher the run-up amplification, which in this case differs significantly from the theoretical curve (Fig.

2.16

), which is calculated for symmetric monochromatic waves.

In the previous cases, we only considered waves that were shorter than the distance from the undisturbed shoreline to the seaward boundary. However, in the piecewise linear bathymetry (Fig. 2.11), we found that resonance is possible for wavelengths larger than the distance mentioned above. The canonical case which we study in this section can be seen as the limiting example of the piecewise linear bathymetry as θ i → 0 , i > 1. Therefore, we decided to perform simulations using a plane beach with slope tan θ = 0.02 connected to a region of constant depth (h 0 = 100 m), . The horizontal extent is two wavelengths offshore from the toe of the beach (λ 0 /L = 3.14, η 0 = 2.5 m). Steeper wavefronts are observed when L 0 = 4λ 0 . Run-up timeseries (d). Waves with steeper wavefront run-up higher. which has a length L 0 = 3000 m. This means that the distance from the initial shoreline to the seaward boundary is L t = 8000 m. We used very small amplitude waves (η 0 /h 0 = 0.001) and we did not put a relaxation zone close to the generation region. For each simulation we sent four non-breaking waves. We can observe in Fig. 2.17 that resonance is possible for wavelengths larger than L t and this result is closer to our observations from the piecewise linear bathymetry.

Discussion

In this chapter, based on the findings of Stefanakis et al.

[2011], we reproduced their results of run-up amplification using milder, more geophysically relevant bottom slopes and we showed that resonant run-up amplification on an infinite sloping beach is found for several waveforms and is robust to modal perturbations.

In Appendix A we also present resonant run-up amplification of cnoidal waves as well as resonant run-up of monochromatic waves on a transect of the Mentawai bathymetry.

Resonant run-up was confirmed by the laboratory experiments of Ezersky et al.

[2013] for monochromatic waves and they also distinguished the resonant run-up frequencies from the natural frequencies of the system. The first resonant regime (λ 0 /L = 5.2, where λ 0 is the incoming wavelength and L is the horizontal beach length) was achieved for non-breaking waves as in [START_REF] Ezersky | Physical simulation of resonant wave run-up on a beach[END_REF] . Moreover, it is also interesting to note that our findings present similarities to those of Bruun and Johannesson [1974] or [START_REF] Bruun | Stability of sloping structures in relation to ξ = tanα/ √ H/L 0 risk criteria in design[END_REF] who described a resonance phenomenon of short wave run-up on sloping structures. They described it as wave breaking taking place at the point of maximum run-down simultaneously with the arrival of the subsequent wave. Here, we do not see wave breaking, but there is a synchronization between the maximum run-down of a wave and the arrival of the next wave (Fig. the impression that tsunamis break just before they reach the shoreline.

The same resonant mechanism is found when the bathymetry is piecewise linear.

However, when the beach is connected to a constant depth region, the picture is different. No resonant regimes are observed when the incoming wavelength is smaller than the distance between the initial shoreline and the seaward boundary.

The maximum run-up amplification is found close to the breaking limit for nearly symmetric low amplitude waves. In that case the linear theory is in close agreement with the results for non-breaking waves. Nevertheless, the steepness of the wavefront plays an important role on run-up, with increasing steepness leading to higher run-up. It is not clear though if it is the wavefront steepness which is responsible for the increase of run-up values or the wave asymmetry (skewness).

Increasing the incoming wavelength more than the wave propagation distance to the undisturbed shoreline results in observing resonant regimes similar to those found in the piecewise linear bathymetry example, which can be thought as the limiting case when the angles θ i → 0 , i > 1. 4 It is of interest to report that [START_REF] Keller | Water wave run-up on a beach[END_REF] tried to reproduce numerically their analytical solution and found a peak which corresponds to the resonant frequency in our simulations.

However, they dismissed these results by saying that their computational scheme was not good enough. On the theoretical side, we can say that in linear theory, the existence or not of resonance depends on the geometry the bathymetry has at the seaward boundary.

The discrepancies of the results using the two bathymetric profiles raise questions Once generated, tsunamis evolve substantially through spreading as seen in Fig. 3 Nonlinear effects become important only as the tsunami evolves over nearshore topography, through shoaling and refraction.

In our analysis, we will not develop exact results for onland inundation; when needed to model few interesting past events, we will use numerical methods to calculate the run-up. We note that, in terms of 1+1 dimensional run-up, Synolakis

[1987] developed an exact solution to both linear and non-linear shallow-water theory for the canonical problem of a non-periodic long wave climbing up a beach.

Synolakis and Skjelbreia [1993] then discussed the evolution of the maximum during shoaling. [START_REF] Carrier | Tsunami propagation from a finite source[END_REF] developed an analytical solution based on the methodology defined by [START_REF] Carrier | Tsunami propagation from a finite source[END_REF] to evaluate the propagation of finite crest length sources of Gaussians over flat bathymetry and discussed the directivity.

Their solution involves computation of complete elliptic integrals of the first kind, with singularities, as also in [START_REF] Carrier | Tsunami run-up and draw-down on a plane beach[END_REF], a difficulty subsequently resolved by [START_REF] Kânoğlu | Nonlinear evolution and runup-rundown of long waves over a sloping beach[END_REF] and [START_REF] Kânoğlu | Initial value problem solution of nonlinear shallow water-wave equations[END_REF]. Moreover, in its current form, the Carrier and Yeh [2005] model can only be used for Gaussians, and cannot even be applied for more standard long wave models such as solitary waves.

Their solution can be extended for N -waves by superposing Gaussians, but the calculation of the singular integrals is challenging and involves approximations.

We note that Tadepalli and Synolakis [1994] proposed a paradigm change for analytical studies of the impact of long waves and introduced N -waves as more realistic initial waveforms for tsunamis. Indeed, the initial waveform of real tsunamis is dipole shaped, with leading-elevation or -depression waves depending on the polarity of the seafloor deformation and the observation location. Marchuk and Titov [1989] described the process of tsunami wave generation by rectangular positive and negative initial ocean surface displacements, and their results suggested unusual amplification.

In summary, we will first present a new general analytical solution for linear shallow water-wave equation for propagation of a finite crest length source without any restriction on the cross section of the initial profile, i.e., Gaussian, solitary, or Nwaves, and without singular elliptic integrals [START_REF] Aydın | Analytical Solutions of Shallow-water Wave Equations[END_REF]. We will then show that focusing points exist for N -wave shaped sources and persist in the corresponding dispersive solutions. Last, we will apply our analytical solution to the 17 July 1998 Papua New Guinea, the 17 July 2006 Java, and the 11 March 2011 Japan events to explain some extreme run-up observations.

Here, we examine focusing and large local amplification, not by considering the effects of underwater diffractive lenses, which anyway are calculable since Berry

[2007], but by considering the dipole nature of the initial profile. We do not of course purport to explain specific patterns of real transoceanic tsunamis, but only to suggest that, in addition to the Berry focusing from bathymetric lenses, large amplification might result from focusing exclusively dependent on the shape and orientation of the initial wave.

We note that our analytical solution is not intended to replace numerical models, that are necessary for identifying the impact of scenario events or historic tsunamis.

Just as analytical results for idealized problems help establish the scaling of natural phenomena, far easier than repeated numerical computations over large parameter ranges, our basic wave theory analysis helps interpret puzzling field observations.

Analytical solution

We use the linear shallow-water wave equation to describe a propagation problem over a constant water depth d as a governing equation. In terms of the free surface elevation η * (x * , y * , t * ), the dimensional governing equation is

η * t * t * -g d (η * x * x * + η * y * y * ) = 0, (3.1)
where g is the gravitational acceleration. Dimensionless variables are introduced as

(x, y) = (x * , y * ) l 0 , η = η * d 0 , and t = t * t 0 . (3.2)
Here, l 0 (d 0 ) = d and t 0 = l 0 / √ gd 0 = √ d/g are the characteristic length (the depth), and the time scales respectively. The dimensionless form of the governing equation then takes the form

η tt -η xx -η yy = 0, (3.3)
with an initial surface profile η(x, y, t = 0) = η 0 (x, y) and zero initial velocity

η t (x, y, t = 0) = 0.
The Fourier transform pair over the space variables (x, y) is

η = ∫ ∞ -∞ ∫ ∞ -∞ η e -i(kx+ly) dx dy (a), η = 1 (2π) 2 ∫ ∞ -∞ ∫ ∞ -∞
η e i(kx+ly) dk dl (b).

(3.4) η = η(k, l, t), η = η(x, y, t), and k and l are the wave numbers in the x and y directions respectively. We transform the governing equation (3.3) into

ηtt + (k 2 + l 2 )η = 0, (3.5)
and the initial conditions to η(k, l, t = 0) = η0 (k, l) and ηt (k, l, t = 0) = 0, using (3.4a). The solution of (3.5) under these conditions is now straightforward;

η(k, l, t) = η0 (k, l) cos ωt where ω = √ k 2 + l 2 . Back-transformation through (3.4b) gives η(x, y, t) = 1 (2π) 2 ∫ ∞ -∞ ∫ ∞ -∞
η0 (k, l) e i(kx+ly) cos ωt dk dl.

(3.6)

We express the finite crested initial waveform as the product of two independent functions,

η 0 (x, y) = f (x) g(y), (3.7)
as in Carrier and Yeh [2005]. f (x) describes the transverse extent of the initial wave profile, g(y) represents the streamwise lateral cross-section of the source, such as a solitary wave or an N -wave (Fig. 3.2). Representation of the initial wave in the form of (3.7) is advantageous, because it allows evaluation of the Fourier transforms of f (x) and g(y) independently;

η0 (k, l) = ∫ ∞ -∞ ∫ ∞ -∞ f (x) g(y) e -i(kx+ly) dx dy = [∫ ∞ -∞ f (x) e -ikx dx ] [∫ ∞ -∞ g(y) e -ily dy ] = f (k) ĝ(l). (3.8)
Unlike [START_REF] Carrier | Tsunami propagation from a finite source[END_REF] who used Gaussians, we prefer hyperbolic functions to define lateral cross-section profiles g(y) such as solitary or N -waves. We also use a hyperbolic function in the transverse direction to define finite crest length

f (x), i.e., f (x) = 1 2 [tanh γ(x -x 0 ) -tanh γ(x -(x 0 + L))]. (3.9)
In (3.9), x 0 is the starting point of the source and L is its crest length, as shown in Fig. 3.2. The parameter γ in (3.9) is determined by the lateral cross-section of the initial wave; it is either γ = γ s for solitary waves or γ = γ n for N -waves. The factor 1 2 is included so that the amplitude of f (x) is equal to unity, in the limit x 0 -→ -∞ and x 0 + L -→ +∞, for a given γ. In that case, the problem reduces to a single propagation direction, and (3.7) represents an infinitely long source.

Given that ∫ ∞ -∞ tanh γx e -ikx dx = -i π γ cosech π 2γ k, (3.10)
(see Appendix B for details), the transform of (3.9) takes the following form

f (k) = i π 2γ (e -ikL -1) e -ikx 0 cosech π 2γ k. (3.11)
Spectra for cross-section profiles g(y) have been given by [START_REF] Synolakis | The runup of solitary waves[END_REF] for solitary waves, and in Tadepalli and Synolakis [1994] for generalized N -waves.

A solitary wave with amplitude H can be described by g s (y) = Hsech 2 γ s (yy 0 ) with γ s = √ 3H/4. Its transform is given by Synolakis (1987) as ĝs (l) = (4π/3) l e -ily 0 cosech α s l, with α s = π/(2γ s ). Consequently, an initial finite crest wave with solitary wave cross-section can be described with

η s (x, y) = H 2 [tanh γ s (x -x 0 ) -tanh γ s (x -(x 0 + L))] sech 2 γ s (y -y 0 ), (3.12)
and its transform is given by ηs (k, l) = i 4 π 3 α s l (e -ikL -1) e -i(kx 0 +ly 0 ) cosech α s k cosech α s l.

(3.13)

The generalized N -wave profile is defined by Tadepalli and Synolakis [1994] as

g n (y) = εH(y -y 2 ) sech 2 γ n (y -y 1 )
, where ε is a scaling parameter which ensures that the initial wave amplitude is H. The steepness of the wave is controlled by the parameter p 0 in γ n = √ 3Hp 0 /4. The locations of depression and elevation parts of an N -wave are controlled by y 1 and y 2 . The transform of the generalized N -wave is given by Tadepalli and Synolakis [1994] as

ĝn (l) = (4εH/π) α 2 n e -ily 1 [(y 1 -y 2 ) l + i(1 -α n l coth α n l)] cosech α n l 1 with α n = π/(2γ n ).
As a result, an initial finite crest wave with an N -wave cross-section is

η n (x, y) = εH 2 [tanh γ n (x -x 0 ) -tanh γ n (x -(x 0 + L))] ×(y -y 2 ) sech 2 γ n (y -y 1 ), (3.14)
with the corresponding transform,

ηn (k, l) = i 4 ε H π α 3 n (e -ikL -1) e -i(kx 0 +ly 1 ) [(y 1 -y 2 ) l + i (1 -α n l coth α n l)] × cosech α n k cosech α n l. (3.15)
In addition, we follow the linear dispersive analytical solution of Kervella et al.

[2007] for the potential flow equation. 2 The solution is similar to (3.6) except the dispersion relation and it is given by

η(x, y, t) = 1 (2π) 2 ∫ ∞ -∞ ∫ ∞ -∞ η0 (k, l) e i(kx+ly) cos ϖt dk dl, (3.16)
where ϖ = √ ω tanh ω, and, again η0

(k, l) = f (k) ĝ(l) with ω = √ k 2 + l 2 .
We have to note that another solution to the same Cauchy problem is given in cylindrical coordinates by the Poisson equation [Dias and Pelinovsky]:

η * ( ⃗ r * , t * ) = 1 2πc ∂ ∂t * ∫ ∫ D η * 0 (⃗ ρ) d⃗ ρ √ c 2 t * 2 -| ⃗ r * -⃗ ρ | 2 , ( 3.17

Results and discussion

In what follows, we will use realistic initial wave profiles to show the existence of focusing points. Then, we will discuss possible consequences of focusing points for the 17 July 1998 Papua New Guinea tsunami. We will then compare our results obtained with linear non-dispersive theory with results from linear dispersive, nonlinear non-dispersive and weakly non-linear weakly dispersive theories. When a typical N -wave initial source (Fig. 3.3a) propagates over a constant depth, the initial wave profile splits into two outgoing waves as presented in Fig. 3.3bc, i.e., leading-elevation N -wave (LEN) and leading-depression N -wave (LDN). This is consistent with the inferences of Tadepalli andSynolakis [1994, 1996] and also with field observations, for example those after the 26 December 2004 Great Sumatran tsunami. In Male, Maldives the tsunami manifested itself as an LEN, as elsewhere to the west of the Sumatran subduction zone. In Phuket, Thailand, it manifested itself as an LDN, as elsewhere to the east of the subduction zone In Fig. 3.6c, we present the maximum wave height envelopes along different directions for the wave propagating into the leading-depression side, where we observe the focusing points in each direction. Sissano Lagoon was the area with the maximum impact where the maximum number of casualties was reported. This region is approximately r = 25km away from the source, between 30 • -45 • from the source orientation close to the 30 • radial line (Fig. 3.6d). As seen from Fig. 3.6c, along 30 • , the shoreline might had faced a wave 1.6 times larger than if the shoreline was r = 50km away. Using the run-up formalism of [START_REF] Synolakis | The runup of solitary waves[END_REF] with R ∼ H 5/4 , the equivalent run-up is 1.8 times greater, than had there been no focusing.

We then investigate the effect of the source crest length (L) has on the location and amplitude of the focusing point (Fig. 3.6e). As L increases, the focusing point moves further away, as expected, since focusing is an effect of the finite crest length. Overall, the maximum wave height along the leading-depression wave side is higher than along the direction of the leading-elevation wave side, at any point past the focusing point. Also, increasing L up to a certain value increases the maximum wave height at the focusing points, then it is constant. In Fig. 3.6f, we present the effect of the steepness parameter over the location of the focusing point and the maximum wave height value at the focusing point. Decreasing the wave steepness translates the focusing point further away. However, it does not change the maximum wave heights at the focusing points.

In addition, we investigate the effects nonlinearity and dispersion have on focusing points, using the PNG initial condition as a test case. When the initial wave is steep (p 0 = 15), the wave is more dispersive, and we observe a slight increase of the maximum at the focusing point (Fig. 3.7). However, all four approximations of shallow-water wave theory -linear non-dispersive, linear dispersive, non-linear nondispersive, weakly non-linear weakly dispersive-produce almost identical results when p 0 = 2, as seen in Fig. 3.8.

Conclusions

We considered three-dimensional long wave propagation over a constant depth basin. We solved the linear shallow-water wave equation as an initial value problem subject to realistic initial conditions, such as N -waves. We showed the existence of focusing point for evolution of dipole sources.

We then discussed unusually pronounced run-up observations from recent events which might have been exasperated by focusing, using the 1998 Papua New Guinea tsunami as particular example. Our results strongly imply that focusing increases the shoreline amplification of the tsunami. Also, preliminary examination of the field survey data and analysis of the 25 October 2010 Mentawai Islands, Sumatra, Indonesia tsunami [START_REF] Hill | The 2010 Mw 7.8 Mentawai earthquake: Very shallow source of a rare tsunami earthquake determined from tsunami field survey and near-field GPS data[END_REF] suggests a similar focusing effect.

We note that the three-dimensional focusing identified here supplements the focusing by bathymetric lenses proposed by [START_REF] Berry | Focused tsunami waves[END_REF] and points to a very rich and complex evolutionary process for long waves in the ocean. While for most practical applications, numerical solutions will be used, event-specific computations can seldom help identify basic wave phenomena. Both focusing mechanisms suggest extreme care when attempting to further regularize ill-posed seismic inversions from run-up measurements, using simplistic multiplicative factors to scale sources. Further, our analysis suggests that the occurrence of extreme run-up at a given local should not always automatically lead to the inference of offshore landslides, at least not before focusing is eliminated.

Chapter 4

Tsunami generation above a sill

Modeling earthquake-triggered tsunami generation and propagation is now standard for hazard analysis of vulnerable coastlines. Over the past fifteen years, it has been possible to numerically calculate the evolution of tsunamis across oceans, develop inundation maps and provide operational warnings. The seafloor displacement is calculated from seismic characteristics based on classic elastic dislocation theory applied on a flat half-space. The motion is considered instantaneous, and it is usually translated to the free surface to define an initial condition for the hydrodynamic problem. This kind of generation is referred to as passive. When the seafloor evolves and interacts continuously with the water surface, the generation is referred to as active or dynamic, as for example is believed to be the case with submarine landslides. However, not all earthquake motions obey standard scaling laws in terms of their rupture times; the so-called slow earthquakes -e.g. the Aleutian 1946, the Nicaragua 1992 [START_REF] Kanamori | The 1992 Nicaragua earthquake: a slow tsunami earthquake associated with subducted sediments[END_REF], & Java 2006 [START_REF] Fritz | Extreme runup from the 17 July 2006 Java tsunami[END_REF]. All have produced larger tsunamis than if they were of normal duration.

Hammack [1972,1973] wave increased with increasing seafloor deformation rate, the maximum free surface displacement was found to be bounded by the maximum bottom displacement that triggered it.

While the rupture time appears to be well constrained from seismic analysis, the rise time is not [START_REF] Synolakis | Damage, conditions of East Java tsunami of 1994 analyzed[END_REF]. Rupture velocities are of the order of thousands of m/sec, rise velocities of m/sec. Dynamic tsunami generation in terms of the seafloor rise but with instantaneous rupture was further studied by [START_REF] Dutykh | Linear theory of wave generation by a moving bottom[END_REF], solving a linear problem using a double Fourier transform in space combined with a Laplace transform in time. They also considered more general seafloor displacements, as used by seismologists, usually based on Okada

[1992] solution for an elastic half space, and always in a uniform depth setting.

They showed that differences exist between passive and dynamic seafloor motions.

Kervella et al.

[2007] compared nonlinear and linear theories for tsunami generation and concluded that in most cases linear theory is sufficient to capture the initial dynamics. They further showed that even for impulsive bottom displacements, passive and dynamic generation produce differences in the initial wave field, a result which was not a priori obvious given the small rise times. gradients are small, then the flat seabed approximation appears reasonable. However, this is seldom the case; in real subduction zones, where tsunamis are often forced in regions with steep gradients and large features. The question arises how large localized bathymetric features such as seamounts may affect the generated tsunami. As an example, the Nias, Indonesia 2005 M8.7 earthquake triggered a small tsunami in comparison to its size, because the highest deformation took place under Nias; in essence, the earthquake uplifted the island along with the surrounding seafloor. While now obvious, in the days following the earthquake it was not [START_REF] Kerr | Model shows islands muted tsunami after latest Indonesian quake[END_REF].

Here we investigate how seamounts or other large seafloor features affect tsunami generation, by solving analytically the linear shallow water equations, for the simplified bathymetry of a cylindrical sill initially sitting on a flat horizontal seabed.

While our model is idealized, our purpose is to provide better qualitative understanding of the dynamics of seafloor displacements, and whether the standard practice of ignoring them in numerical studies of tsunami evolution is always adequate.

In Section 4.1 we briefly introduce the LSWE and then solve it in section 4.2, in the transformed space. We return to the physical space in section 4.3, and discuss results for various geometries and dynamical properties of the seafloor displacement in section 4.4.

The field equation

We consider waves forced over a bathymetry (Figure 5.1) consisting of a radius r * c cylinder of height h * c (starred quantities are dimensional), resting on a horizontal seafloor of constant depth h * 2 . Over the sill surface, when r * < r * c , the depth is

h * 1 .
The seafloor is forced underneath the cylinder (r * < r * c ). The forced LSWE in polar coordinates (r * , θ) is

ζ * t * t * -gh * ( 1 r * ζ * r * + ζ * r * r * + 1 r * 2 ζ * θθ ) -gh * r * ζ * r * = f * t * t * , (4.1)
where ζ * is the free surface elevation, h * (r * ) is the water depth from the deforming seafloor to the initially undisturbed water surface, not the entire flow depth, g is the gravitational acceleration, t * is the time, r * is the radial distance from the center of the cylinder and f * is the forcing term, which corresponds to the time-dependent bottom displacement. Subscripts denote differentiation of the relevant variable with respect to the subscript. Now, consider the following non-dimensionalization:

r = r * /r * c , t = t * √ g h * 2 r * c , ζ = ζ * /f * 0 , f = f * /f * 0 , h = h * /h * 2 , (4.2)
where f * 0 is the maximum vertical displacement of the seafloor. The forced nondimensionalized LSWE becomes

ζ tt -h ( 1 r ζ r + ζ rr + 1 r 2 ζ θθ ) -h r ζ r = f tt . (4.3)
We use f (r, t) = (1-e -γt )H(r c -r)H(t), t > 0 as forcing; r c = 1, H is the Heaviside function, and the parameter γ defines the rate of the displacement. Other forms of reasonable forcing can also be used, with minimal change in our development.

We require the solution to be bounded inside the whole fluid domain, at all times: |ζ(r, θ, t)| < ∞. Initially, the free surface is at rest: ζ(r, θ, 0) = 0 and ζ t (r, θ, 0) = 0.

To solve Eq. ( 4.3), we apply the Laplace transform pair We now split the fluid domain into two subregions, the near field where r * < r * c or r < 1, and the far field where r * > r * c or r > 1. We will solve separately in each subregion the forced long-wave equation (Eq. 4.5), and then we will match the solutions at the common boundary r = 1.

ζ(r, θ; ω) = ∫ ∞ 0 ζ(r, θ, t)e -ωt dt, ζ(r, θ, t) = 1 2πi ∫ c+i∞ c-i∞ ζ(
( 1 + r h r h ) ζr - ω 2 r 2 h ζ + ζθθ = - r 2 h ftt , ( 4

Solution in the transformed space

4.2.1 The near field (r < 1)

In the near field, h * = h * 1 , therefore h = h * 1 /h * 2 .
= c d . Furthermore, since both the bathymetry and the forcing are axisymmetric, ζ θθ = 0. in Eq. (4.5). Given that the water depth is constant in the near field, h r = 0. Therefore, Eq. (4.5) becomes,

r 2 ζrr + r ζr - ω 2 r 2 c d ζ = - r 2 c d ftt , (4.7)
which is an inhomogeneous second order partial differential equation. We solve it with the method of variation of parameters, by considering first the relevant homogeneous equation,

r 2 ζrr + r ζr - ω 2 r 2 c d ζ = 0 . (4.8)
With the change of variables χ = ωr/ √ c d , the modified Bessel equation of zeroth order is revealed, whose two independent solutions are the modified Bessel functions I 0 (χ) and K 0 (χ). Therefore, the general solution to Eq. (4.8) is From [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], the Wronskian of the two independent solutions of the homogeneous equation (Eq. 4.8) is

ζh (r; ω) = α 1 I 0 ( ωr √ c d ) + β 1 K 0 ( ωr √ c d ) , ( 4
W {I 0 (χ); K 0 (χ)} = - 1 χ . (4.11)
Thus, the solution of the inhomogeneous (Eq. (4.8)) is

ζ(r; ω) = α 1 I 0 ( ωr √ c d ) -P (r; ω) , (4.12)
where

P (r; ω) = ∫ r 0 ftt (ω) c d W (ρ) [ I 0 ( ωρ √ c d ) K 0 ( ωr √ c d ) -I 0 ( ωr √ c d ) K 0 ( ωρ √ c d )] dρ (4.13)
is the particular solution, with the Wronskian of the two homogeneous solutions W (ρ) = -1/ρ and ftt (ω) = γω/(γ + ω) is the forcing term. Integration yields the simpler expression,

P (r; ω) = ftt (ω) ω 2 ( I 0 ( ωr √ c d ) - 1
) . (4.14)

The far field (r > 1)

In the far field, the flow is axisymmetric and h * = h * 2 or h = 1 and there is no direct forcing. Consequently, Eq. (4.5) can be simplified:

r 2 ζrr + r ζr -ω 2 r 2 ζ = 0 . (4.15)
The above equation is similar to Eq. (4.17)

The coefficients α 1 and β 2 will be obtained by the matching conditions at r = 1.

Matching at r = 1

The near and far field solutions are matched at the common boundary r = 1. We require continuity of the free surface elevation ζ and the radial fluxes h ζr . These assumptions have been shown by [START_REF] Bartholomeusz | The reflexion of long waves at a step[END_REF] to give correct results for the reflexion coefficient of long waves incident on a step, even though the vertical acceleration might not be small locally. Therefore,

α 1 I 0 ( ω √ c d ) -P (1; ω) = β 2 K 0 (ω) , c d α 1 ω √ c d I 1 ( ω √ c d ) -c d P r (1; ω) = -β 2 ω K 1 (ω) . (4.18)
By solving the linear system of equations (4.18), we obtain the expressions for α 1 and β 2 (the dependence on c d is omitted for brevity):

α 1 (ω) = ω K 1 (ω) P (1; ω) + c d K 0 (ω) P r (1; ω) ω K 1 (ω) I 0 ( ω √ c d ) + ω √ c d K 0 (ω) I 1 ( ω √ c d ) , (4.19)
β 2 (ω) = -ω √ c d I 1 ( ω √ c d ) P (1; ω) + c d I 0 ( ω √ c d ) P r (1; ω) ω K 1 (ω) I 0 ( ω √ c d ) + ω √ c d K 0 (ω) I 1 ( ω √ c d ) . (4.20)
By replacing the expression for P (r; ω) (Eq. 4.14) and by noting that

P r (r; ω) = ftt (ω) ω 2 ω √ c d I 1 ( ωr √ c d ) , (4.21)
we find

α 1 (ω) = ftt (ω) ω 2   1 - K 1 (ω) K 1 (ω) I 0 ( ω √ c d ) + √ c d K 0 (ω) I 1 ( ω √ c d )   (4.22)
and

β 2 (ω) = ftt (ω) ω 2 √ c d I 1 ( ω √ c d ) K 1 (ω) I 0 ( ω √ c d ) + √ c d K 0 (ω) I 1 ( ω √ c d ) = ftt (ω) ω 2 β(ω) . (4.23)
The expression of the transformed free surface elevation in the near field can be further simplified:

ζ(r; ω) = ftt (ω) ω 2   1 - K 1 (ω) I 0 ( ωr √ c d ) K 1 (ω) I 0 ( ω √ c d ) + √ c d K 0 (ω) I 1 ( ω √ c d )   = ftt (ω) ω 2 [ 1 -α(ω) I 0 ( ωr √ c d )] . (4.24)
Finally, the transformed free surface elevation in the complex half-plane

| arg(ω)| < π/2 is ζ(r; ω) =        ftt(ω) ω 2 [ 1 -α(ω) I 0 ( ωr √ c d )] for r < 1 ftt(ω) ω 2 β(ω) K 0 (ωr) for r > 1 (4.25)

Wave description

To obtain results in the physical space, we apply the inverse Laplace transform.

We first consider the free-surface elevation in the far field (r > 1).

4.3.1

The far field (r > 1)

Inverse transforming (4.25), in the far field, yields

ζ(r, t) = 1 2πi ∫ c+i∞ c-i∞ ftt (ω) ω 2 β(ω) K 0 (ωr) e ωt dω . (4.26)
For t > 0, contour integration in the complex plane (Appendix C) results into

∫ c+i∞ c-i∞ ftt (ω) ω 2 β(ω) K 0 (ωr) e ωt dω = i ∫ ∞ -∞ - ftt (-is) s 2 β(-is) K 0 (-isr) e -ist ds . (4.27)
Further replacing s with ω yields

ζ(r, t) = 1 2π ∫ ∞ -∞ - ftt (-iω) ω 2 β(-iω) K 0 (-iωr) e -iωt dω , (4.28)
in the far field. By noting that

K n (x) = π 2 i n+1 H n (ix) I n (x) = 1 i n J n (ix)      when -π < arg(x) ≤ π 2 (4.29)
where H n is the Hankel function of the first kind and order n and J n is the Bessel function of the first kind and order n, the free surface elevation in the far field is

ζ(r, t) = 1 2π ∫ ∞ -∞ B(ω) H 0 (ωr) e -iωt dω , (4.30)
where

B(ω) = - iωγ γ -iω √ c d ω 2 J 1 ( ω √ c d ) H 1 (ω) J 0 ( ω √ c d ) - √ c d H 0 (ω) J 1 ( ω √ c d ) . (4.31)
For large distances, we can find the form of the leading wave by the asymptotic expansion of the Hankel function for large arguments

H n (x) = √ 2 πx e i(x-πn/2-π/4) . (4.32)
We should also note that the wave celerity in the far field is C = 1, and thus, r = t. By replacing the previous in Eq. ( 4.30) we get

ζ l (t) = 1 √ 2π 3 t ∫ ∞ -∞ B(ω) 1 √ ω e -i π 4 dω , (4.33)
which suggests that the leading wave amplitude decays as O(t -1/2 ), similarly to a tsunami generated by a Gaussian-shaped landslide on a conical island [START_REF] Renzi | Landslide tsunamis propagating around a conical island[END_REF]. For the special case of an initially flat seabed (c d = 1), the previous expression can be further simplified since

H 1 (ω) J 0 (ω) -H 0 (ω) J 1 (ω) = -πω/(2i).
We also replace e -iπ/4 = (1i)/ √ 2, and we obtain

ζ l (t) = 1 4 √ πt ∫ ∞ -∞ γ(1 -i) γ -iω J 1 (ω) √ ω dω .

The near field (r < 1)

Following similar steps, we take the inverse Laplace transform of the free-surface elevation

ζ(r, t) = 1 2πi ∫ c+i∞ c-i∞ ftt (ω) ω 2 [ 1 -α(ω) I 0 ( ωr √ c d )] e ωt dω , (4.35)
which, as shown in Appendix C, simplifies to

ζ(r, t) = 1 2π ∫ ∞ -∞ - ftt (-iω) ω 2 [ 1 -α(-iω) I 0 ( -iωr √ c d )] e -iωt dω . (4.36)
By using the relations (4.29) we obtain

ζ(r, t) = 1 2π ∫ ∞ -∞ - ftt (-iω) ω 2 [ 1 -A(ω) J 0 ( ωr √ c d )] e -iωt dω , (4.37)
where

A(ω) = H 1 (ω) H 1 (ω) J 0 ( ω √ c d ) - √ c d H 0 (ω) J 1 ( ω √ c d ) . (4.38)
Above the center of the sill, the free surface elevation is

ζ(0, t) = 1 2π ∫ ∞ -∞ - ftt (-iω) ω 2 [1 -A(ω)] e -iωt dω , (4.39)
and in the special case of an initially flat seafloor (c d = 1), it can be further simplified to

ζ(0, t) = 1 2π ∫ ∞ -∞ iγ γ -iω 1 ω [ πω 2i H 1 (ω) + 1 ] e -iωt dω . (4.40)
In Appendix D, we present a solution using Fourier transforms. While counterintuitive, [START_REF] Mei | The applied dynamics of water waves[END_REF] has applied both integration kernels for transient progressive waves. As we show in Appendix D, both integral transforms produce mathematically identical results, which is reassuring.

Results

To validate our solution, we compare it with the solution obtained by Hammack [1972], for the limiting case c d = 1, which is the case of no sill in place. This corresponds to uplift of a circular region within a flat seafloor which otherwise remains stationary.

In Fig. 4.2, we compare the time series of free surface elevations obtained by our solution and Hammack's, at r = 0 and r = 1, for three different values of the rate of bottom deformation γ. We also include results from nonlinear, nondispersive theory, using the NSWE solver VOLNA [Dutykh et al., 2011b], to check whether nonlinearity is important. As the figures suggest, it is not, for this choice of parameters, with one exception we will discuss later. We note that Hammack's solution is based on potential flow theory and it is dispersive, ours is not. The figure suggests that dispersion is more apparent, as the bottom displacement becomes more rapid. The amplitude of the leading wave increases with γ, until it reaches unity (i.e. equal to the maximum bottom displacement). High values of γ correspond to impulsive generation. Moreover, to leading order, the predictions of the far field heights of the forced wave are identical among the three approximations: linear nondispersive, nonlinear nondispersive and linear dispersive. In the near field, the predictions for the maximum height are identical for the impulsive motion, which is of greater physical interest.

The effect of non uniform water bathymetry on the wave formation is depicted in Fig. 4.3. The higher the sill (lower c d ), the more energy is trapped above the sill, which results in higher maximum values of free surface elevation at r = 0, in the transitional regime (Fig. 4.3 (a)) and wider wave crests in the impulsive regime (Fig. 4.3 (d)). In addition, when the bathymetry is not flat, we observe secondary crests of the free surface elevation above the center of the sill, which become larger with increasing sill height and faster bottom motion (larger γ). These secondary crests are formed as the wave travels from r = 0 to r = 1, and is partially reflected back from the edge of the sill. Thus the time lag between the leading crest and the second is t ≈ 2/ √ c d .

Increasing the sill height, results in deeper backdrop of the free surface after the formation of the leading crest. In the far field (Fig. 4.3 (c,f)), the leading wave has lower amplitude and longer wavelength as c d decreases, again due to wave trapping above the sill. The differences are neither spectacular nor negligible.

Wave trapping is visible in Fig. 4.4 (b), which is an extreme case with the sill occupying 90% of the total water depth. In this case, directional focusing of the waves reflected back from the edge of the sill towards the centre results in amplification of the free surface elevation at r = 0.

Discussion

We have explored analytically long wave generation above a cylindrical sill by solving the forced linear shallow water equations. We used both Fourier and Laplace transforms and have gotten identical expressions for the free surface elevation.

Our results are in good agreement with those of [START_REF] Hammack | Tsunamis -A Model of Their Generation and Propagation[END_REF] who derived a solution without a sill, even though our solution is non-dispersive. Dispersion becomes important when the bottom displacement is instantaneous. According to our model, increasing the height of the sill, results in reduced maximum wave height in the far field. This is intuitive, because as the sill height increases, a smaller body of water is set into motion. This had not yet been shown analytically and helps explain the observation of Arcas and Synolakis reported in [START_REF] Kerr | Model shows islands muted tsunami after latest Indonesian quake[END_REF] for the Nias earthquake.

Moreover, due to partial wave trapping above the sill, smaller trailing waves are observed in the far field, the amplitudes of which grow with increasing sill height.

During wave generation, above the surface of the sill, the maximum wave height increases with increasing sill height and is always less or equal to the maximum bottom displacement.

However, due to wave trapping, directional focusing from the edge of the sill to the center is observed. When the sill height is large and the generation tends to be impulsive, the wave is amplified above the center of the sill. Nevertheless, further investigation is needed for this case, to assess the effect of nonlinearity. The behavior is reminiscent of a weak Anderson localization, which has been observed in other classic wave theory phenomena [e.g. [START_REF] Schwartz | Transport and Anderson localization in disordered two-dimensional photonic lattices[END_REF], if indeed the analogy is more than superficially correct.

Our results imply that when calculating the initial condition for tsunami generation, care must be taken to account for the effects of large submarine features, such as sills. Initializing computations only with the seafloor deformation obtained from standard elastic models for dislocations, as is universally the practice in tsunami modeling [START_REF] Synolakis | Validation and verification of tsunami numerical models[END_REF], may, on occasion, underestimate the initial condition, which is the sum of the elastic solution and that of the sill. Given that seamounts are invariably of much smaller area than the entire deformed region, this effect is likely small, at geophysical scales. However, if extremely nearshore, seafloor features may affect the run-up, and this needs to be carefully quantified.

Chapter 5

Can small islands protect nearby coasts from tsunamis? An active experimental design approach

In recent years we have witnessed the dreadful damage tsunamis caused in coastal Island [START_REF] Yeh | The Flores island tsunamis[END_REF][START_REF] Yeh | Tsunami catastrophe in Babi island[END_REF] focused scientists attention on tsunami run-up on islands and the studies that followed, which included both laboratory experiments , 1998], showed that long waves can cause extensive run-up on the lee side of a conical island.

Earlier studies [START_REF] Homma | On the behaviour of seismic sea waves around circular island[END_REF][START_REF] Lautenbacher | Gravity wave refraction by islands[END_REF][START_REF] Longuet-Higgins | On the trapping of wave energy round islands[END_REF][START_REF] Smith | Scattering of surface waves by a conical island[END_REF][START_REF] Vastano | Tsunami response for islands: Verification of a numerical procedure[END_REF]] have given some insight on the behavior of long waves around conical islands, but did not deal with run-up. The big conclusion of all the aforementioned studies is the fact that long waves do not behave as wind generated waves and that small islands which would act as natural barriers in normal sea conditions, transform into amplifiers of wave energy in areas believed to be protected and where coastal communities thrive. Furthermore, recent findings [START_REF] Hill | The 2010 Mw 7.8 Mentawai earthquake: Very shallow source of a rare tsunami earthquake determined from tsunami field survey and near-field GPS data[END_REF] have shown enhanced tsunami run-up in areas which lied behind small islands in the vicinity of the mainland and therefore were supposedly protected.

In recent years, the developments in computer science and the increase of computational power in combination with the smaller associated cost compared to laboratory experiments, have led scientists to more and more rely on numerical simulations. However, each simulation has a computational cost, which increases with model complexity and spatiotemporal resolution. Therefore, a series of experiments which have a specific objective, such as maximization/minimization of an output, should be carefully designed in order to reach the desired conclusion with the least number of experiments. Thus, finding the argmax x f (x) where f (x)

is the output of the experiment depending on the parameters x is not trivial since we do not know the analytical expression of f (x) and therefore it should be approximated. The difficulty of the problem increases with the number of parameters on which the output depends and the "naive" approach to create regular grids and test all the points becomes prohibitively expensive.

For this reason, sampling techniques have been developed which aim to reduce the number of points by finding a representative sample of the input space. These techniques are commonly referred to as "Experimental Design" -e.g. Sacks et al.

[1989] -and are static, meaning that the design (sampling) is made initially, before the execution of the experiments and the selection of the future query points is not guided by the experimental results. At the end all the points are queried. This is already a great advancement compared to the regular grids approach.

More recently, adaptive design [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] From a physical point of view, the current study aims to elucidate the tsunami run-up on a plane beach behind a small conical island compared to an adjacent lateral location on the beach not directly influenced by the presence of the island.

To achieve that we use numerical simulations of the nonlinear shallow water equations (NSWE). Moreover, we will present a newly developed method for Active

Experimental Design [START_REF] Contal | Parallel Gaussian Process optimization with pure exploration[END_REF], which we will apply to our problem and we will also discuss its advantages and limitations in a more general setting.

Finally, we will present some metrics that can be used for the comparison of the performance of different learning strategies and an empirical stopping criterioni.e. a criterion which will signal the achievement of the optimization objective.

Experimental Configuration

Simulations

The simplified bathymetric profile consists of a conical island sitting on a flat bottom and a plane beach behind the island (Fig. 5.1). The height of the crest of the island above still water level is always fixed at 100m. The distance between the seaward boundary and the toe of the island is also fixed at 7600m. A single wave profile is prescribed as forcing at the seaward boundary, having the form η 0 (t) = 1.5 sech 2 (ωt -2.6). We use this formulation because we want to avoid the solitary wave link between the water depth and the wave amplitude as is discussed in Madsen et al. [2008] and in [START_REF] Madsen | Analytical solutions for tsunami runup on a plane beach: single waves, N-waves and transient waves[END_REF]. The problem is governed by 5 physical parameters, namely the island slope, the plane beach slope, the water depth, the distance between the island and the beach and the prescribed incident wavelength which is controlled by ω (Table 5.1).

The numerical simulations were performed using VOLNA Dutykh et al. [2011b] which solves the NSWE. VOLNA can simulate the whole life cycle of a tsunami from generation to run-up. It uses a Finite Volume Characteristic Flux scheme [START_REF] Ghidaglia | Une méthode volumes-finis à flux caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de conservation[END_REF][START_REF] Ghidaglia | On the numerical solution to two fluid models via cell centered finite volume method[END_REF]] with a MUSCL type of reconstruction for higher order terms [START_REF] Kolgan | Application of the minimum-derivative principle in the construction of finite-difference schemes for numerical analysis of discontinuous solutions in gas dynamics[END_REF][START_REF] Kolgan | Finite-difference schemes for computation of three dimensional solutions of gas dynamics and calculation of a flow over a body under an angle of attack[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme V: a second order sequel to Godunov' method[END_REF]] and a third order Runge-Kutta time discretization. The code uses an unstructured triangular mesh, which can handle arbitrary bathymetric profiles and can also be refined in areas of interest.

The mesh resolution that we used varied from 500m at the seaward boundary to 2m at the areas where we measured run-up (Fig. 5.2).

The run-up was measured on the plane beach exactly behind the island and on a lateral location on the beach, which was far enough from the island and thus was not directly affected by its presence (Fig. 5.2). To compute run-up, 11 equally spaced virtual wave gauges were positioned at each location. The gauge location has an inherent uncertainty due to spatial discretization, which is minimized with the use of higher resolution around these locations (Fig. 5.2). The actual horizontal spacing of the wave gauges was dependent on the beach slope. The minimum height of the gauges was the still water level and the maximum height was selected to be 5.5m above the undisturbed water surface. The run-up never exceeded this height in any of the simulations. The maximum run-up is defined as the maximum recorded wave height at the highest wave gauge. When the wave did not reach the height of a gauge, then that gauge did not record any signal.

Experimental Design

In order to fill the input parameter space we had to choose the input points in such a way that maximal information is obtained with a moderate number of points.

This procedure is known as "Experimental Design" and it is a passive approach as we described in the Introduction. This is the first step to reduce the computational cost. For this purpose we used Latin Hypercube Sampling [START_REF] Mckay | Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] with maximization of the minimum distance between points. When using the Latin Hypercube Sampling (LHS) of a function of M variables, the range of each variable is divided into N equally probable, non-overlapping intervals. Then one value from each interval is randomly selected for every variable. Finally, a random combination of N values for M variables is formed. The maximization of the minimum distance between points is added as an extra constraint. The LHS is found to lead to better predictions than regular grids when used with multivariate emulators [START_REF] Urban | A comparison of Latin hypercube and grid ensemble designs for the multivariate emulation of an Earth system model[END_REF]. In order to accurately cover the input space, we ran 200 simulations selected by LHS. For comparison a regular grid approach with 10 grid points in each dimension would require 10 5 points/simulations. However, we can further improve the performance of this approach with the so called "Active Experimental Design", which can suggest the order in which the points will be queried and is result driven, as we describe in the following section. We should clarify that the Active Experimental Design does not require to be initialized with LHS and it can work on any set of points. The LHS was used for evaluation purposes, to better fill the parameter space with a limited number of points, because we do not have the luxury to employ a large number of random points. Finally, we need to stress that in terms of statistical learning, our strategy is not restrained to the specific tsunami research problem and thus can be applied to a wide spectrum of disciplines where the objective is scalar optimization with cost constraints.

Active Experimental Design

Active Batch Optimization

Let f : K → R be an unknown function, with K a compact subset of R d , which

we can evaluate for a fixed cost at any location x ∈ X , where X is a finite subset of K. We address the problem of sequentially finding the maximizer of f ,

x ⋆ = argmax x∈X f (x) ,
in the lowest possible number of evaluations. The arbitrary choice of formulating the optimization problem as a maximization is without loss of generality, as we can obviously take the opposite of f if the problem is a minimization one. We consider the special case where K evaluations of f can be acquired with no increase in cost. For example when f is the result of a numerical experiment which can be computed in parallel on a machine with K cores, and the cost to minimize is computation time. At each iteration t, we choose a batch of K points in X called the queries {x k t } 0≤k<K , and then observe simultaneously the observations of f at these points, potentially noisy,

y k t = f (x k t ) + ϵ k t ,
where the ϵ k t is independent Gaussian noise N (0, σ 2 ). The stochasticity in the measurements is due to the discretization, both spatial and temporal and not due to the accuracy of the evaluation model, since we consider it to be deterministic.

Objective

Assuming that the number of iterations allowed, hereafter called horizon T , is unknown, a strategy has to be good at any iteration. Care must be taken to tackle the exploration/exploitation tradeoff, that is balance learning the function f globally with focusing around the predicted maximum. We aim to minimize the cumulative regret [START_REF] Bubeck | Regret analysis of stochastic and nonstochastic multi-armed bandit problems[END_REF],

R K T = ∑ t<T ( f (x ⋆ ) -max k<K f (x k t )
) .

The loss r K t incurred at iteration t is the simple regret for the batch {x k t } k<K [START_REF] Bubeck | Pure exploration in multi-armed bandits problems[END_REF], defined as

r K t = f (x ⋆ ) -max k<K f (x k t ) .
A strategy is said to be "no-regret", when

R K T T ---→ T →∞ 0 .

Gaussian Processes

In order to analyze the efficiency of a strategy, we have to make some assumptions on f . We want extreme variations of the function to have low probability. Modeling f as a sample of a Gaussian Process (GP) is a natural way to formalize the intuition that nearby locations are highly correlated. It can be seen as a continuous extension of multidimensional Gaussian distributions. We said that a random process f is Gaussian with mean function m and non-negative definite covariance C given by m and k. That is, ∀n < ∞, ∀x 1 , . . . , x n ∈ X ,

(f (x 1 ), . . . , f (x n )) ∼ N (µ, C) , with µ[x i ] = m(x i )
and C[x i , x j ] = k(x i , x j ) .

If we have the prior knowledge that f is drawn from a GP with known kernel function k, then, based on the observations of f after T iterations, the posterior distribution remains a Gaussian process, with mean µ T and variance σ 2 T , which can be computed via Bayesian inference by [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF],

µ T (x) = k T (x) ⊤ C -1 T Y T (5.1)
and

σ 2 T (x) = k(x, x) -k T (x) ⊤ C -1 T k T (x) .
(5.2)

X T = {x k t } t<T,k<K , Y T = [y k t ] x k t ∈X T
are the set of queried locations and the vector of noisy observations respectively.

k T (x) = [k(x k t , x)] x k t ∈X T
is the vector of covariances between x and the queried points and

C T = K T + σ 2 I with K T = [k(x, x ′ )]
x,x ′ ∈X T the kernel matrix and I stands for the identity matrix.

The three most common kernel functions are:

• the polynomial kernels of degree α ∈ N,

k(x 1 , x 2 ) = (x ⊤ 1 x 2 + c) α , c ∈ R.
• the (Gaussian) Radial Basis Function kernel (RBF or Squared Exponential)

with length-scale l > 0,

k(x 1 , x 2 ) = exp ( - ∥x 1 , x 2 ∥ 2 2l 2 ) , (5.3)
• the Matérn kernel, of length-scale l and parameter ν,

k(x 1 , x 2 ) = 2 1-ν Γ(ν) ( √ 2ν ∥x 1 , x 2 ∥ l ) ν K ν ( √ 2ν ∥x 1 , x 2 ∥ l ) , (5.4)
where K ν is the modified Bessel function of the second kind and order ν.

The Bayesian inference is represented on

Parallel Optimization Procedure

Now that we have set the statistical background, we can describe the learning strategy that we used, namely the Gaussian Process Upper Confidence Bound with Pure Exploration algorithm [START_REF] Contal | Parallel Gaussian Process optimization with pure exploration[END_REF].

Confidence Region

A key property from the GP framework is that the posterior distribution at a location x has a normal distribution N ( µ t (x), σ 2 t (x)). We can then define a upper confidence bound f + and a lower confidence bound f -, such that f is included in the interval with high probability:

f + t (x) = µ t (x) + √ β t σ t-1 (x) and
(5.5) f + and f -are illustrated in Fig. 5.3 by the upper and lower envelope of the grey area respectively. The region delimited in that way, the high confidence region, contains the unknown f with high probability.

f - t (x) = µ t (x) - √ β t σ t-1 (x) (5.6)

Relevant Region

We define the relevant region R t being the region which contains x ⋆ with high probability. Let y • t be our lower confidence bound on the maximum,

x • t = argmax x∈X f - t (x)
and

y • t = f - t (x • t ) .
y • t is represented by the horizontal dotted green line on Fig. 5.4. R t is defined as, t . The relevant region R is shown in light green (without edges). The first query x 0 is the maximizer of f + . We show in dashed line the upper and lower bounds with the update of σ after having selected x 0 . The second query x 1 is the one maximizing the uncertainty inside R.

R t = { x ∈ X | f + t (x) ≥ y • t } . ( 5

GP-UCB-PE

We present here the Gaussian Process Upper Confidence Bound with Pure Exploration algorithm (GP-UCB-PE), a very recent algorithm from [START_REF] Contal | Parallel Gaussian Process optimization with pure exploration[END_REF] combining two strategies to determine the queries {x k t } k<K for batches of size K. The first location is chosen according to the UCB rule,

x 0 t = argmax x∈X f + t (x) . (5.8)
This single rule is enough to deal with the exploration/exploitation tradeoff. The value of β t balances between exploring uncertain regions (high posterior variance σ 2 t (x)) and focusing on the supposed location of the maximum (high posterior mean µ t (x)). This policy is illustrated with the point x 0 in Fig. 5.4.

The K -1 remaining locations are selected via Pure Exploration restricted to the region R t . We aim to maximize I t , the information gain about f granted by the K -1 points [START_REF] Cover | Elements of Information Theory[END_REF]]. This can be efficiently approximated by the greedy procedure which selects the K -1 points one by one and never backtracks 1 .

1 Formally, It is the gain in Shannon entropy H when knowing the values of the observations at those points, conditioned on the observations we have seen so far, It(X) = H(Y) -H(Y | Xt) . Finding the K -1 points that maximize It is known to be intractable [START_REF] Ko | An exact algorithm for maximum entropy sampling[END_REF] and thus an approximation is required.

The location of the single point that maximizes the information gain is easily computed by maximizing the posterior variance. Our greedy strategy selects for each 1 ≤ k < K the following points one by one,

x k t = argmax x∈Rt σ (k) t (x) , (5.9)
where σ

(k) t

is the updated variance after choosing {x k ′ t } k ′ <k . We use here the fact that the posterior variance does not depend on the values y k t of the observations, but only on their position x k t . One such point is illustrated with x 1 in Fig. 5.4.

These K -1 locations reduce the uncertainty about f , improving the guesses of the UCB procedure by x 0 t . The overall procedure is shown in Algorithm 1.

Theoretical Guarantees

With the Gaussian process assumption on f , we can adjust the parameter β t such that f will be contained by the high confidence region with high probability.

Under this condition, Contal et al.

[2013] prove a general bound on the regret achieved by GP-UCB-PE, making this strategy a "no-regret" algorithm. The order of magnitude of the cumulative regret R K T we obtained with a linear kernel is √ T K log T K, and for RBF Kernel √ T K (log T K) d , up to polylog factors. When K ≪ T , these probabilistic bounds with parallel queries are better than the ones incurred by the sequential procedure by an order of √ K.

Stopping criterion

One challenging problem we face in practice is to decide when to stop the iterative strategy. The theoretical analysis only gives general expected bounds when T tends to infinity, it does not provide estimations of the constant and short term factors. We have to design an empirical, yet robust, criterion. One trivial solution is to fix the number of iterations (or computation time) allowed by a predefined limit. This is not a suitable solution for the general case, as one does not know precisely the amount of exploration needed to be confident about the maximum of f .

Other rules like the criteria based on the local changes of the queries in the target space (improvement-based criteria) as well as in the input space (movement-based criteria), come from the convex optimization literature. These are not suitable in our global, nonconvex setting, as they stops after a local maximum is found. t , unknown in a real situation, and the rank correlation ρ Xv (π t-1 , π t ) (in log-scale), for the synthetic function Himmelblau. The stopping threshold ρ 0 was set to 10 -4 , and the lag ℓ to 4, the algorithm stopped at iteration 12, after having found a good candidate at iteration 7 and the true maximum at iteration 11.

Our approach is to stop when the procedure ceases to learn relevant information on f . We attempt to measure the global changes in the estimator µ t between two successive iterations, with more focus on the highest values. The algorithm then stops when these changes become insignificant for a short period. The change between µ t and µ t+1 is measured by the correlation between their respective values on a finite validation data set X v ⊂ X .

We

x∈Xv µ t (x)) = n v .
We then define the discounted rank dissimilarity d Xv and the normalized rank correlation ρ Xv as

d Xv (π t , π t+1 ) = ∑ x∈Xv ( π t+1 (x) -π t (x) ) 2 ( π t+1 (x) ) 2 , ρ Xv (π t , π t+1 ) = 1 - d Xv (π t , π t+1 ) max π + ,π -∈Sn v d Xv (π + , π -) .
The denominator in the definition of ρ Xv represents the discounted rank dissimilarity between two reversed ranks π + and π -. The normalized rank correlation for such two ranks will therefore be equal to 0, whereas for any rank π, this correlation with itself will be ρ Xv (π, π) = 1. The correlation ρ Xv can be seen as a modified Spearman's rank correlation coefficient, where the squared distances are weighted by their position in the new rank π t+1 . If we note π{n ↔ m} the inversion of the n th and m th ranks in π, for all rank π, we remark that

ρ Xv (π, π{2 ↔ 3}) > ρ Xv (π, π{1 ↔ 2}).
We show in Fig. 5.5 the observed relationship between the regret incurred at iteration t and the rank correlation ρ Xv (π t-1 , π t ). The algorithm stops when ρ Xv (π t , π t+1 ) stays below a given threshold ρ 0 for ℓ iterations in a row. The value of this threshold has to be fixed empirically. In Fig. 5.5, ρ 0 = 10 -4 and ℓ = 4, the algorithm stopped at iteration 12.

In Fig. 5.6, we can see the distribution of the final number of iterations T (lower is better) together with the final gap G T = min t<T R K t (lower is better) for 4 different thresholds ρ 0 . The value ρ 0 = 10 -4 appears to be a good threshold, since the final regret is always 0, and the number of iterations remains low. Further reduction of this threshold will again guarantee zero regret, but with higher computational cost. On the other hand, greater values of ρ 0 will reduce the computational cost, but with an increased risk to miss the maximum.

Experiments

Synthetic data sets

Apart from the tsunami experiment, which is 5-dimensional and we do not know a priori the form of the response surface, in order to test the performance of the active learning algorithm, it is wise to use some synthetic data sets. These data sets can be easily visualized (2-dimensional) and we can attribute to them some desired properties, such as several local maxima or background noise, which aim to test the algorithm and will give us a direct feedback of its behavior.

Himmelblau function

The Himmelblau data set is a nonconvex function in dimension 2. We compute a slightly tilted version of the Himmelblau function, and take the opposite to match the challenge of finding its maximum.

Gaussian mixture

This synthetic function comes from the addition of three 2-D Gaussian functions. at (0.2, 0.5), (0.9, 0.9), and the maximum at (0.6, 0.1). We then perturb these Gaussian functions with smooth variations generated from a Gaussian Process with Matérn Kernel (Eq. 5.4) and very few noise. It is shown on Fig. 5.7(b). The highest peak being thin, the sequential search for the maximum of this function is quite challenging.

Assessment

We verify empirically the performance of GP-UCB-PE by measuring the decay of the regret obtained on several synthetic functions. For the sake of convenience we do not report the cumulative regret R K t on the figures, but the gap between the maximum point discovered and the true maximum, defined as the minimum regret so far G K t = min t ′ ≤t r k t .

First, we show in Fig. 5.9 the impact of the size of the batch K on the minimum regret G K t . It is shown that the sequential (the red curve K = 1) performs better than the rest, without however being extremely outperforming (orders of magnitude). Therefore in a situation where the final number of queries is the restrictive factor, one would better choose small batch sizes. On the other hand, if total time of the optimization is the restrictive factor, then one could choose larger batch sizes without sacrificing too much computational cost.

We then compare our approach to three other strategies,

• Random, which chooses the next queries {x k t } k<K at random,

• Exploration, which attempts to maximize the information gain on f at each iteration,

x k t = argmax x∈X σ (k) t (x) ,
• Exploitation, which only focuses on the predicted maximum,

x k t = argmax x∈X \{x k ′ t } k ′ <k µ t (x) .
For all data sets and algorithms, the batch size K was set to 10 and the learners were initialized with a random subset of t init = 20 observations (x i , y i ). The curves in Fig. 5.8 show the evolution of the gap G K t in term of iteration t. We report the average value with the confidence interval over 64 experiments (random initializations). The kernel function used was always an RBF kernel (Eq. 5.3).

The parameters of the algorithm, like the length-scale of k (represented by l), were chosen as the best parameters found by validation on a random subsample of the data.

Our learning algorithm is shown to outperform the rest strategies on the synthetic data sets. Exploitation in these data sets looses time because it gets stuck in a local maximum, while Exploration and Random strategies will asymptotically find the global maximum on average. On the Tsunami data set, Exploitation perform slightly better than GP-UCB-PE probably due to the simplicity of the run-up function, which even though is 5-dimensional does not seem to pose any serious challenges, probably because it has only one maximum.

The Effect Of The Conical Island

After running 200 simulations, we have found that in none of the situations considered the island did offer protection to the coastal area behind it. On the contrary, we have measured amplified run-up on the beach behind it compared to a lateral location on the beach, not directly affected by the presence of the island (Fig. 5.10).

This finding shows that small islands in the vicinity of the mainland will act as amplifiers of long wave severity at the region directly behind them and not as natural barriers as it was commonly believed so far. The maximum amplification achieved was ∼ 70% more than were the island absent and the median amplification factor is 1.3. The island focuses the wave on its lee side, while far from it the

((a)) ((b)) ((c)) ( (d))
((e)) ((f))

Figure 5.11: Snapshots of the free surface elevation measured in meters as the wave passes the island and runs up the beach behind it. The island focuses the wave on its lee side and the amplified wave propagates towards the beach. The colorbar is in logarithmic scale for visualization purposes. In the present case the run-up amplification is 1.59 .

wave propagates unaffected (Fig. 5.11). The amplified wave propagates towards the beach and causes higher run-up in the region directly behind the island.

One of the key questions is which parameters control the run-up amplification (RA)

and in what way. To answer these questions, we can use the statistical model. Figure 5.12: Local sensitivity of the maximum run-up amplification on (a) the island and beach slopes, (b) the distance between the island and the beach and the water depth, (c) the cyclic frequency of the wave. The range of the above parameters can be found in Table 5.1.

We perform a local sensitivity analysis around the maximum RA by fixing all parameters except one each time at the value which corresponds to the maximum RA and we vary the excluded parameter across the whole range of its input space (Fig. 5.12). We can observe that some parameters vary more than others and thus are more important. These are the water depth, the beach slope and the cyclic frequency of the wave. Having said that, one would wonder why the parameters of the island do not seem to be that important. The answer might not be simple, because dependencies could be hidden in the correlations of the input parameters.

To better understand these dependencies, it is of interest to recombine the input parameters in order to obtain nondimensional but physically interpretable measures. In Fig. 5.13 we express the RA as a function of the ratio of the wavelength over the island radius at its base λ 0 /r 0 and the Iribarren number J computed using the beach slope and normalized with the relative wave amplitude H 0 /h 0 . We see that the RA strongly depends on the ratio λ 0 /r 0 and that the highest values are attained when the wavelength is almost equal to the island radius. The normalized Iribarren number gives a satisfactory classification, with smaller values qualitatively leading to higher RA. Of course the complexity of the problem is superior and cannot be completely explained by the previous two measures. Nevertheless, Fig. 5.13 shows that a comprehensive knowledge of the system can give better insight than pure statistics and implies that interdisciplinary problems like this one should be treated with close collaboration between the various fields.

Conclusions

We examined the effect the presence of a small conical island has on the long-wave run-up on a plane beach behind it. Using a simplified geometry dependent on five physical parameters, we wanted to find the combination of parameters which will give us the maximum run-up amplification with a minimal computational cost.

To achieve that, we employed an active experimental design strategy based on a Gaussian Process surrogate model. The strategy, which is parallelizable, can handle efficiently the tradeoff between exploration of the input space and focusing on the region where the argmax x f (x) is believed to reside in.

Even though our algorithm is asymptotically convergent, we are interested in its behavior for a finite time horizon T . Comparing our strategy to other commonly used ones, we showed that it performs better in most cases. Overall, the active experimental design approach can reduce the computational cost more than 60% compared to a classic experimental design (LHS) and potentially much higher (e.g. Fig. 5.8(c)). In addition, the computational gain is orders of magnitude smaller than a regular grid approach -3 orders of magnitude for a 5-dimensional problem.

Moreover, a stopping criterion was presented, which can signal the achievement of the optimization objective and thus the end of the experiments. The development of such a criterion is essential in real applications where only a small number of experiments is allowed due to cost constraints and thus the theoretical asymptotic convergence is useless. Our stopping criterion is based on the difference in the ranking of the predictions of the surrogate model between two consecutive iterations. Even though, it is shown to correlate well with the regret r K t (Fig. 5.5), it depends on an empirically set threshold. Therefore, more research is needed to develop a more robust stopping criterion, which will either be derived directly from the learning algorithm or will relate the threshold to the dimensionality of the problem.

The active learning strategy is not restricted to tsunami research and can be applied to a wide range of problems and disciplines where the optimization should be balanced with a reasonable computational or actual cost. Another interesting perspective is to incorporate in the optimization not only physical parameters, but also numerical ones, such as the spatial discretization, the placement of the (virtual) sensors and others. The inclusion of these numerical parameters can be handled by the GP-UCB-PE algorithm. Finally, further research is needed for the development of active learning algorithms for multi-objective optimization and pareto front tracking.

From a physical point of view, our results show that for the given setup and range of input parameters, the island instead of protecting the beach behind it, as it was widely believed so far, it acts as a focusing lens of wave energy on its lee side.

Until now, the prevailing practice in studying maximum run-up for civil defense applications has been that a plane beach provides the worst possible condition for wave amplification and thus, small offshore islands were believed to offer protection to coastal areas in their vicinity.This finding is of fundamental importance for the correct education of coastal communities and thus their preparedness in case of a tsunami. In order to evaluate the integral on the semicircle C δ , we need to evaluate the behavior of ζ(r; -is) for small arguments, using the small argument approximations for I 0,1 and K 0,1 (Eq. C.12):

ζ(r; -is) ≈ iγ γ -is 1 s [ 1 - 2 2 + s 2 ln(-is) ] (C.25)
Consequently, by making use of the parametric transform s = δe iϕ , ϕ ∈ (0, π) and letting δ → 0, we get and because the water depth is constant, its radial derivative is zero. Therefore Eq. (D.2) can be simplified:

∫ C δ ζ(r; -is) e -ist ds = = lim δ→0 ∫ 0 π - γ γ -iδe iϕ
r 2 ζrr + r ζr + ω 2 r 2 c d ζ = - r 2 c d ftt (D.4)
which is an inhomogeneous second order partial differential equation. In order to solve it, we will apply the method of variation of parameters. Thus, we first consider the relevant homogeneous equation Similarly as before, in the far field we have axial symmetry in the bathymetry and the water depth does not vary with r and is h * = h * 2 or h = 1. Moreover, there is no direct influence from the forcing term on the fluid motion. Consequently, Eq. (D.2) can be simplified:

r 2
ζrr + r ζrω 2 r 2 ζ = 0 (D.12)

The above equation is a standard Bessel equation of zeroth order. The two independent solutions could be again the Bessel functions J 0 (ωr) and Y 0 (ωr) or the Hankel functions of the first and second kind and order zero, H

0 (ωr) and H

D.2 Wave Description

In order to go back to physical space, we need to apply the inverse Fourier transform, which is straightforward: The last two expressions are identical to Eqs. (4.37) and (4.30) respectively.

ζ(r, t) = ∫ +∞ -∞ ftt (ω)

1. 1 2 . 9

 129 The accumulation of strain at interplate boundaries and earthquake triggered tsunami generation. Source: Atwater [2005] . . . . . . . . 1.2 Shoaling of a tsunami and run-up. . . . . . . . . . . . . . . . . . . . 1.3 Schematic of the DART buoy system. Source: NOAA . . . . . . . . 1.4 Sea wall in Japan. Source: www.marineinsight.com . . . . . . . . . 1.5 Enhanced run-up behind small islands in the vicinity of the shore of Sibigau, Mentawai Islands during the 2010 tsunami. From the numerical computations of Hill et al. [2012] . . . . . . . . . . . . . . 1.6 The bathymetry off the coast of Mentawai Islands. A large seamount is apparent. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Geometry of the problem of the run-up of transient long waves along a sloping beach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The geometry of a plane beach connected to a region of constant depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Run-up height of a sine wave arriving from a basin of constant depth. The solid line is formula (2.29) and the dashed line is the result for an infinite beach of constant slope (2.26). . . . . . . . . . 2.4 Maximum run-up amplification R max /η 0 (a) and maximum horizontal velocity amplification (b) of monochromatic waves on a plane beach with respect to nondimensional wavelength for three different slopes, namely tan θ = 0.02 ; 0.05 ; 0.1 (L = 5000 m). Resonance is observed when the incoming wavelength is approximately 2.4 and 5.2 times the beach length. . . . . . . . . . . . . . . . . . . . . . . 2.5 Evolution of the volume of fluid V inside the computational domain during resonance (α = 0.02, L = 5000 m, λ 0 /L = 5.2). V i is the initial volume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Spatio-temporal behaviour of non-dimensional horizontal velocity u/(g tan θL) 1/2 in the resonant regime (a) and non-resonant regime (b). The black line describes the evolution of the shoreline position in time. In both cases tan θ = 0.05 and L = 5000 m. . . . . . . . . . 2.7 Maximum run-up of monochromatic waves on a plane beach as a function of nondimensional wavelength for two different slopes, namely tan θ = 0.13; 0.26 (L = 12.5 m). The results were obtained with VOLNA, a 2D finite volume solver of . . . . . . . . . . . . . . 2.8 Maximum run-up of bichromatic waves with respect to nondimensional frequency (tan θ = 0.13 , L = 12.5 m). . . . . . . . . . . . . . xi List of Figures xii Comparison of a typical monochromatic wave and a corresponding semi-randomly perturbed signal both in physical space (top) and Fourier space (bottom), where N is the Fourier mode and h(k) is the spectral amplitude. Time t is in seconds and the free-surface elevation η 0 (t) is in meters. . . . . . . . . . . . . . . . . . . . . . . . 2.10 Run-up timeseries of a perturbed resonant monochromatic wave when the slope is 0.02 and L = 5000 m. . . . . . . . . . . . . . . . . 2.11 Schematic of the piecewise linear bathymetry . . . . . . . . . . . . . 2.12 Plot of the zeros of the determinant in the denominator of Eq. (2.46) as a function of ω and α 2 when α 1 = 0.02, L 1 = 5000 m and L 2 = 6000 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Run-up amplification of monochromatic waves on a piecewise linear bathymetry consisting of two segments as a function of the nondimensional wavelength when α 1 = tan θ 1 = 0.02, α 2 = tan θ 2 = 0.01, L 1 = 5000 m and L 2 = 6000 m. . . . . . . . . . . . . . . . . . 2.14 Maximum run-up amplification as a function of nondimensional wavelength for the canonical case (η 0 = 1.25 m , h 0 = 100 m , tan θ = 0.02). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 Snapshots of free surface elevation over the constant depth region (a,b,c) . The horizontal extent is two wavelengths offshore from the toe of the beach (λ 0 /L = 3.14, η 0 = 2.5 m). Steeper wavefronts are observed when L 0 = 4λ 0 . Run-up timeseries (d). Waves with steeper wavefront run-up higher. . . . . . . . . . . . . . . . . . . . . 2.16 Maximum run-up amplification as a function of nondimensional wavelength for the canonical case (η 0 = 2.5 m, h 0 = 100 m, tan θ = 0.02). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Maximum run-up amplification as a function of nondimensional wavelength for the canonical case when L t = 8000 m is the distance from the undisturbed shoreline to the seaward boundary (η 0 = 0.1 m, h 0 = 100 m, tan θ = 0.02). . . . . . . . . . . . . . . . . . . . . . 3.1 A map with maximum wave amplitudes for the 11 March 2011 Japan tsunami based on a real time forecast (Tang et al., 2012) using the non-linear shallow-water wave equations solver MOST [Titov and Synolakis, 1998]. Colour-filled contours show predicted maximum tsunami amplitudes in deep water. The green star shows the earthquake epicentre location. . . . . . . . . . . . . . . . . . . . 3.2 Definition sketch: (left inset) three-dimensional and (right inset) top views. Not to scale. . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Definition sketch for focusing. Evolution of an N -wave source over a constant depth calculated using the MOST model; (a) initial wave, (b-c) evolution, and (d) maximum amplitude at each grid point. . . 3.4 (a) Two-and (b) three-dimensional evolution of an N -wave (top insets) with H = 0.001, L = 30, p 0 = 15, y 1 = 0, y 2 = 2.3, and ε = 0.04 (for the two-dimensional case ε = 0.07) over a constant depth. Time evolution results for two-and three-dimensional propagations are given at t = 20 and t = 60, including maximum wave envelopes -maximum wave height for the entire time at each spatial location-(thick lines). Note that the three-dimensional results are given along the x-bisector line. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.5 (a) Maximum wave envelopes and time evolution at (b) t = 20 and (c) t = 60 using linear non-dispersive (solid line), linear dispersive (dashed line), non-linear non-dispersive (triangles) and weakly nonlinear weakly dispersive (squares) theories. Refer to the caption of Fig. 3.4 for the initial wave parameters. . . . . . . . . . . . . . . . 50 3.6 (a) The PNG source is defined as in (3.14) with H = 0.01, p 0 = 15, L = 1, γ n = 0.34, y 1 = 50, y 2 = 50.2, x 0 = 49.5, and ε = 4.93. Dimensionless quantities are calculated using the reference depth 1600m. (b) Maximum wave amplitude at each grid point calculated using MOST. (c) Maximum wave amplitude envelopes along the β = 0 • -60 • lines, with 15 • increments. r originates from the point where the maximum initial wave height is located, i.e., (x, y) = (50, 52.5) in inset (a). Dots indicate locations of focusing points. (d)

3 . 6

 36 for the initial wave parameters. . . . . . . . . . . . . . . . . . . 54 4.1 Definition sketch: A radial transect of the fluid domain and the bathymetry in physical coordinates. . . . . . . . . . . . . . . . . . . 58 4.2 Comparison of our solution with predictions from the NSWE solver VOLNA and Hammack's solution when there is no sill (c d = 1), at r = 0 (a -c) and r = 1 (d -f) for different values of the rate of seafloor deformation γ. In the top panel γ = 0.1, which corresponds to the creeping regime, in the middle panel γ = 1 corresponds to the transitional regime and in the bottom panel γ = 10, which corresponds to the impulsive regime. . . . . . . . . . . . . . . . . . 4.3 The effect of bathymetry on the wave formation. Time histories of the free surface elevation for different heights of the sill at r = 0 (top panel), r = 1 (middle panel) and r = 5 (bottom) panel. The left column (a -c) corresponds to a bottom motion with γ = 1 , and the right column (d -f) corresponds to a bottom motion with γ = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Spatio-temporal plots of the free surface evolution when c d = 1 (a) and c d = 0.1 (b) for an impulsive bottom motion with γ = 10. When the sill is very high, wave trapping occurs and thus little amount of energy leaks to the far field. . . . . . . . . . . . . . . . . 5.1 Schematic of the geometry of the experimental setup. . . . . . . . . 5.2 The unstructured triangular grid. Colors represent bathymetric contours. The areas of high grid density on the beach, are the locations of run-up measurements. . . . . . . . . . . . . . . . . . . . 5.3 Gaussian Process inference of the posterior mean µ (blue line) and variance σ based on 4 realizations (blue crosses). The high confidence region (area in grey) is delimited by f + and f -. . . . . . . . 5.4 Two queries of GP-UCB-PE on the previous example. The lower confidence bound on the maximum is represented by the horizontal dotted green line at y • t .

7

 7 and the true maximum at iteration 11. . . . . . . . . . . . . . . . . . . . . . . . 5.6 Distribution of the final number of iterations T and the final gap (minimum regret) G T for 4 different thresholds. . . . . . . . . . . . 5.7 Visualization of the synthetic data sets used for assessment. . . . . 5.8 Experiments on several real and synthetic tasks. The curves show the decay of the mean of the simple regret r K t with respect to the iteration t, over 64 experiments. We show with the translucent area the 95% confidence intervals. . . . . . . . . . . . . . . . . . . . . . 5.9 Impact of the batch size K on the gap G K t . The curves show the mean of 64 experiments on the synthetic data Gaussian Mixture, with t init = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Histogram of the run-up amplification on the beach directly behind the island compared to the run-up on a lateral location on the beach, not directly affected by the presence of the island. . . . . . . . . . . 5.11 Snapshots of the free surface elevation measured in meters as the wave passes the island and runs up the beach behind it. The island focuses the wave on its lee side and the amplified wave propagates towards the beach. The colorbar is in logarithmic scale for visualization purposes. In the present case the run-up amplification is 1.59 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 Local sensitivity of the maximum run-up amplification on (a) the island and beach slopes, (b) the distance between the island and the beach and the water depth, (c) the cyclic frequency of the wave. The range of the above parameters can be found in Table 5.1. . . . . . . 5.13 Run-up amplification (RA) as a function of the wavelength to the island radius (at its base) ratio. The color code indicates the surf similarity (Iribarren number) computed with the beach slope and multiplied with the wave nonlinearity (wave height to water depth ratio). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 Maximum run-up amplification of cnoidal waves as a function of parameter m for an infinite sloping beach (L = 4000m, tan θ = 0.13). A.2 (a) Transect of the Mentawai bathymetry (Left Boundary: 99.8333 o E -3.2333 o N, Right Boundary: 100.2333 o E -2.8167 o N). (b) Amplification ratio as a function of non-dimensional wavelength. . . C.1 Contour integration path. . . . . . . . . . . . . . . . . . . . . . . . C.2 Contour plot of the zeros of the real (solid line) and imaginary (dashed line) parts of the denominator in the expression for β 2 as a function of complex ω. Since the two curves do not intersect, the denominator of β 2 does not have any complex zeros. . . . . . . . . . C.3 Transformed contour integration path. . . . . . . . . . . . . . . . . Chapter 1

  by the Storegga Slides, 100km north-west of the coast of Møre [Bondevik et al., 2003, Bryn et al., 2005, Dawson et al., 1988, O'Brien et al., 2013]. Tsunamis are long waves which have wavelength order magnitude of 10km and amplitude order magnitude 1m in the open ocean, where they travel with speeds of approximately 700 -800km/h. They can transfer energy to very far distances and it is not uncommon to observe tsunamis crossing the Indian or the Pacific Oceans, as happened during the most recent event in Japan on March 11, 2011.

Figure 1 . 1 :

 11 Figure 1.1: The accumulation of strain at interplate boundaries and earthquake triggered tsunami generation. Source: Atwater [2005]

  events have been attributed to (coseismic) landslides, such as the 1998 Papua New Guinea tsunami [Synolakis et al., 2002, Tappin et al., 2008]. A more unusual form of tsunamis are the so-called meteotsunamis [Monserrat et al., 2006, O'Brien et al., 2013], which are formed by atmospheric disturbances in the open ocean and amplify in coastal regions through resonance mechanisms. They are much less energetic than seismic tsunamis and thus are always local.

  , the operational code used by NOAA; the Cornell Multi-grid Coupled Tsunami model (COMCOT) [Liu et al., 1998] and VOLNA, the code developed by Dutykh et al.

Figure 1 . 2 :

 12 Figure 1.2: Shoaling of a tsunami and run-up.

Figure 1 . 3 :

 13 Figure 1.3: Schematic of the DART buoy system. Source: NOAA

Figure 1 . 4 :

 14 Figure 1.4: Sea wall in Japan. Source: www.marineinsight.com

Figure 1 . 5 :

 15 Figure 1.5: Enhanced run-up behind small islands in the vicinity of the shore of Sibigau, Mentawai Islands during the 2010 tsunami. From the numerical computations of Hill et al. [2012]

Figure 1 . 6 :

 16 Figure 1.6: The bathymetry off the coast of Mentawai Islands. A large seamount is apparent.

1994[

  Tsuji et al., 1995], Java 2006 [Fritz et al., 2007], Chile 2010 [Fritz et al., 2011] and Japan 2011 [Grilli et al., 2012]. Hence, a question rises whether these extreme run-up values are related to non-leading waves.

  introduced the notion of bay resonance. Agnon and Mei [1988], Grataloup and Mei [2003] studied the long wave resonance due to wave trapping and wave-wave interactions. Munk et al. [1964] and Rabinovich and Leviant [1992] studied wave resonance in the context of shelf resonance, which occurs when tidal waves have a wavelength four times larger than the continental shelf width. Fritz et al. [2007] suggested that the extreme run-up values measured after the Java 2006 tsunami could be explained by a submarine landslide triggered by the earthquake. All of the above underline the critical role of bathymetry and coastal geometry to long wave propagation and run-up. In a recent study, Kânoğlu et al. [2013] argued that finite-crest length effects may produce focusing. Nonetheless, resonant run-up has already been documented for the case of short waves [Bruun and Günbak, 1977] with an interesting description: "[Resonant run-up] occurs when run-down is in a low position and wave breaking takes place simultaneously and repeatedly close to that location." Similar observations have been made by Stefanakis et al. [2011].On a theoretical basis, the main mathematical difficulty of the run-up problem is the moving shoreline. Progress was made through the introduction of the Carrier and Greenspan (CG) transformation[START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF]] which leads to the reduction of the two Nonlinear Shallow Water Equations (NSWE) into one linear, but the ingenuity of this transformation is that in the transformed space the moving shoreline is static. With the aid of the CG transformation several significant contributions were made to the long wave run-up problem[START_REF] Antuono | The boundary value problem for the nonlinear shallow water equations[END_REF][START_REF] Brocchini | Integral flow properties of the swash zone and averaging[END_REF][START_REF] Carrier | Gravity waves on water of variable depth[END_REF][START_REF] Kânoğlu | Initial value problem solution of nonlinear shallow water-wave equations[END_REF][START_REF] Keller | Water wave run-up on a beach[END_REF][START_REF] Synolakis | The runup of solitary waves[END_REF][START_REF] Tadepalli | The run-up of N-waves on sloping beaches[END_REF]].

Figure 2 . 1 :

 21 Figure 2.1: Geometry of the problem of the run-up of transient long waves along a sloping beach

  with a fixed boundary, unlike the initial equations. The linear cylindrical wave equation (2.13) is well-known in mathematical physics, and its solution can be presented in various forms (Green's function, Hankel and Fourier transforms). Using its solution, the wave field in "physical" variables can be found from algebraic manipulations. Detailed analyses of the wave transformation and run-up have been performed for various initial conditions, see for instance Carrier et al. [2003], Kânoğlu and Synolakis [2006], Tadepalli and Synolakis [1996].

Figure 2 . 2 :

 22 Figure 2.2: The geometry of a plane beach connected to a region of constant depth

Figure 2 . 3 :

 23 Figure 2.3: Run-up height of a sine wave arriving from a basin of constant depth. The solid line is formula (2.29) and the dashed line is the result for an infinite beach of constant slope (2.26).

Figure 2 . 4 :

 24 Figure 2.4: Maximum run-up amplification R max /η 0 (a) and maximum horizontal velocity amplification (b) of monochromatic waves on a plane beach with respect to nondimensional wavelength for three different slopes, namely tan θ = 0.02 ; 0.05 ; 0.1 (L = 5000 m). Resonance is observed when the incoming wavelength is approximately 2.4 and 5.2 times the beach length.

Figure 2 . 5 :

 25 Figure 2.5: Evolution of the volume of fluid V inside the computational domain during resonance (α = 0.02, L = 5000 m, λ 0 /L = 5.2). V i is the initial volume.

Figure 2 . 6 :

 26 Figure 2.6: Spatio-temporal behaviour of non-dimensional horizontal velocity u/(g tan θL) 1/2 in the resonant regime (a) and non-resonant regime (b). The black line describes the evolution of the shoreline position in time. In both cases tan θ = 0.05 and L = 5000 m.

Figure 2 . 7 :

 27 Figure 2.7: Maximum run-up of monochromatic waves on a plane beach as a function of nondimensional wavelength for two different slopes, namely tan θ = 0.13; 0.26 (L = 12.5 m). The results were obtained with VOLNA, a 2D finite volume solver of

R max / η 0 Figure 2 . 8 :

 028 Figure 2.8: Maximum run-up of bichromatic waves with respect to nondimensional frequency (tan θ = 0.13 , L = 12.5 m).

Figure 2 . 9 :

 29 Figure2.9: Comparison of a typical monochromatic wave and a corresponding semirandomly perturbed signal both in physical space (top) and Fourier space (bottom), where N is the Fourier mode and h(k) is the spectral amplitude. Time t is in seconds and the free-surface elevation η 0 (t) is in meters.

Figure 2 .

 2 Figure 2.10: Run-up timeseries of a perturbed resonant monochromatic wave when the slope is 0.02 and L = 5000 m.

Figure 2 .

 2 Figure 2.11: Schematic of the piecewise linear bathymetry

Figure 2 .Figure 2 .

 22 Figure 2.12: Plot of the zeros of the determinant in the denominator of Eq. (2.46) as a function of ω and α 2 when α 1 = 0.02, L 1 = 5000 m and L 2 = 6000 m.

Fig. 2 .

 2 Fig.2.14 that theory and computations do not agree. However, the computations qualitatively follow the trend of the laboratory experiments presented byAhrens [1981] even though they refer to irregular wave run-up.

Figure 2 . 14 :

 214 Figure 2.14: Maximum run-up amplification as a function of nondimensional wavelength for the canonical case (η 0 = 1.25 m , h 0 = 100 m , tan θ = 0.02).

Figure 2 . 15 :

 215 Figure 2.15: Snapshots of free surface elevation over the constant depth region (a,b,c). The horizontal extent is two wavelengths offshore from the toe of the beach (λ 0 /L = 3.14, η 0 = 2.5 m). Steeper wavefronts are observed when L 0 = 4λ 0 . Run-up timeseries (d). Waves with steeper wavefront run-up higher.

Figure 2 . 16 :

 216 Figure 2.16: Maximum run-up amplification as a function of nondimensional wavelength for the canonical case (η 0 = 2.5 m, h 0 = 100 m, tan θ = 0.02).

Figure 2 . 17 :

 217 Figure 2.17: Maximum run-up amplification as a function of nondimensional wavelength for the canonical case when L t = 8000 m is the distance from the undisturbed shoreline to the seaward boundary (η 0 = 0.1 m, h 0 = 100 m, tan θ = 0.02).

3

 3 about the role boundary conditions play both physically and numerically, and more importantly about the character of the flow (stationary vs transient). The problem of long wave run-up has been attacked primarily from a stationary point of view in the past. The well-known solution (2.24) is a standing wave solution and as such does not exhibit any net propagation of energy over time. The solutions we investigated numerically are transient and as such can exhibit an amplification of energy over time.ChapterFocusing of finite-crested N-waves over constant depthThe devastating effects of tsunamis, near-and far-field, became widely recognized following the 2004 Boxing Day tsunami. Run-up measurements over the periphery of the Indian Ocean showed considerable variation in near-and far-field impact.This variation, while it can be inferred from direct numerical simulations of the evolution of the initial wave, remains largely counterintuitive[START_REF] Synolakis | Runup measurements of the December 2004 Indian Ocean tsunamis[END_REF]].

Figure 3 . 1 :

 31 Figure 3.1: A map with maximum wave amplitudes for the 11 March 2011 Japan tsunami based on a real time forecast (Tang et al., 2012) using the non-linear shallowwater wave equations solver MOST [Titov and Synolakis, 1998]. Colour-filled contours show predicted maximum tsunami amplitudes in deep water. The green star shows the earthquake epicentre location.

. 1 ,

 1 for the 2011 Japan event similar to the 2004 tsunami. In their map of global propagation patterns of the 26 December 2004 tsunami,[START_REF] Titov | The global reach of the 26 December 2004 Sumatra tsunami[END_REF] observed two main factors affecting directionality: the configuration of the source region[START_REF] Marchuk | Source configuration and the process of tsunami waves forming[END_REF]] and the waveguide structure of mid-ocean ridges[START_REF] Koshimura | Propagation of obliquely incident tsunamis on a slope, Part 1: Amplification of tsunamis on a continental slope[END_REF]. Continental shelves also act as waveguides[START_REF] González | Edge wave and non-trapped modes of the 25 April 1992 Cape-Mendocino tsunamis[END_REF] and likely caused the persistent ringing for the Pacific coasts of South and North America in the 2004 tsunami. Similar effects were observed during the 2010 Chilean and the 2011 Japan tsunamis; in both cases, strong currents persisted in California ports for days[START_REF] Barberopoulou | Reassessing the tsunami risk in major ports and harbors of California I: San Diego[END_REF][START_REF] Wilson | Observations and impacts from the 2010 Chilean and 2011 Japanese tsunamis in California (USA)[END_REF].Most studies of past tsunamis have concentrated on modelling specific historic or scenario events, and have invariably focused on numerically estimating coastal effects and devastation. The standard practice is for initial conditions to be established directly from estimates of the seismic parameters, then their subsequent evolution is deterministic. Coastal inundation involves often supercritical overland flow depths[START_REF] Fritz | Indian Ocean tsunami flow velocity measurements from survivor videos[END_REF][START_REF] Fritz | The 2011 Japan tsunami current velocity measurements from survivor videos at Kesennuma Bay using LiDAR[END_REF], as observed in the numerous videos of the 11 March 2011 Japan tsunami, and remains the most temperamental aspect of computations, as small coastal features affect flooding patterns to first order[START_REF] Kânoğlu | Long wave runup on piecewise linear topographies[END_REF]].It is thus not surprising why analyses of directivity and focusing in the open ocean remain few. During the 2011 Japan tsunami, Guam, 2750km from the source, experienced maximum run-up of less than 0.60cm, while in Irian Jaya, 4500km away, the reported run-up reached 2.6m. It is our objective here to supplement the few existing substantive studies of the physics of deep-sea evolution of tsunamis and suggest that some counterintuitive observations may be explainable through the classic field theory. Tsunamis triggered by submarine earthquakes have a finite crest (strip) length, which is believed to be calculable adequately from estimates of the seismic parameters and scaling relationships. The initial profile is dipole-shaped, directly reflecting regions of uplift and subsidence of the seafloor, yet the generated wavefield exhibits finger-like radiation patterns, a process often referred to as directivity.Ben-Mehanem [1961] defined a directivity function using the source length and the rupture velocity. Later, Ben-Menahem and Rosenman[1972] used linear theory to calculate the radiation pattern from an underwater moving source and showed that tsunami energy radiates primarily at a direction normal to a rupturing fault.[START_REF] Okal | A field survey of the 1946 Aleutian tsunami in the far fields[END_REF] reported field observations of the 1946 Aleutian tsunami in the far-field, and concluded that a large slow earthquake and a landslide must have occurred concurrently to have caused the observed far-field distribution and nearfield run-up.[START_REF] Okal | Normal mode energetics for far-field tsunamis generated by dislocations and landslides[END_REF] then identified differences in directivity patterns between tsunamis from landslides and dislocations.Directivity arguments alone, however, cannot explain the complexity of the radiated patterns in oceans with trenches and seamounts.[START_REF] Berry | Focused tsunami waves[END_REF] discovered how such underwater features may concentrate tsunamis into cusped caustics, causing large local amplifications at specific focal points. He used linear dispersive theory which describes fairly accurately the evolution of tsunamis in the open ocean.

Figure 3 . 2 :

 32 Figure 3.2: Definition sketch: (left inset) three-dimensional and (right inset) top views. Not to scale.

  )where ⃗ r * and ⃗ ρ are two-dimensional vectors in the (x * , y * ) plane, c = √ gd 0 is the long wave celerity, and integration is performed over the disk D centered at ⃗ r * andof radius ct * | ⃗ r * -⃗ ρ | 2 < c 2 t * 2 . (3.18)This solution can be solved analytically for axisymmetric initial displacements with use of Hankel transforms, but for more general wave profiles it has to be solved numerically.

Figure 3 . 3 :

 33 Figure 3.3: Definition sketch for focusing. Evolution of an N -wave source over a constant depth calculated using the MOST model; (a) initial wave, (b-c) evolution, and (d) maximum amplitude at each grid point.

Figure 3 . 4 :

 34 Figure 3.4: (a) Two-and (b) three-dimensional evolution of an N -wave (top insets) with H = 0.001, L = 30, p 0 = 15, y 1 = 0, y 2 = 2.3, and ε = 0.04 (for the twodimensional case ε = 0.07) over a constant depth. Time evolution results for two-and three-dimensional propagations are given at t = 20 and t = 60, including maximum wave envelopes -maximum wave height for the entire time at each spatial location-(thick lines). Note that the three-dimensional results are given along the x-bisector line.

[

  [START_REF] Synolakis | Tsunami science before and beyond Boxing Day 2004[END_REF]]. More interestingly, however, as the LEN and LDN travel in opposite directions, in the path of the LDN, a positive wave from the centre of elevation part and two positive waves from the sides of depression arrive simultaneously at a point along the bisector line, as shown in Fig.3.3b-c. This is the focusing point, and in its vicinity abnormal tsunami wave height is observed (Fig.3.3d)3 .

Figure 3 . 3 . 2 . 1

 3321 Figure 3.5: (a) Maximum wave envelopes and time evolution at (b) t = 20 and (c) t = 60 using linear non-dispersive (solid line), linear dispersive (dashed line), non-linear non-dispersive (triangles) and weakly non-linear weakly dispersive (squares) theories.Refer to the caption of Fig.3.4 for the initial wave parameters.

Figure 3 . 6 :

 36 Figure 3.6: (a) The PNG source is defined as in (3.14) with H = 0.01, p 0 = 15, L = 1, γ n = 0.34, y 1 = 50, y 2 = 50.2, x 0 = 49.5, and ε = 4.93. Dimensionless quantities are calculated using the reference depth 1600m. (b) Maximum wave amplitude at each grid point calculated using MOST. (c) Maximum wave amplitude envelopes along the β = 0 • -60 • lines, with 15 • increments. r originates from the point where the maximum initial wave height is located, i.e., (x, y) = (50, 52.5) in inset (a). Dots indicate locations of focusing points. (d) The initial N -wave is located at (x * , y * ) = (27, 34)km and is tilted 5 • to be consistent with Synolakis et al. [2002]. Sissano Lagoon, where most of the damage was observed, is located approximately r = 25km away from the initial wave location between the 30 • -45 • lines. Triangles over the 30 • and 45 • lines show the focusing points when p 0 = 5. (e) Maximum wave height envelopes for the source lengths L = 1 (solid line), 10 (dashed line), 20 (dash-dotted line) and 30 (dotted line). Triangles represent MOST numerical results. (f ) Maximum wave height envelopes p 0 = 2 (dashdotted line), 5 (dashed line) and 15 (solid line). When p 0 and L are parametrized, ε is modified to ensure the same maximum wave amplitude as the original initial wave profile.

Figure 3 . 7 :

 37 Figure 3.7: (a) Maximum wave envelope and time evolution at (b) t = 5 and (c) t = 30 for the PNG source, using linear non-dispersive (solid line), linear dispersive (dashed line), non-linear non-dispersive (triangles) and weakly non-linear weakly dispersive (squares) theories. Refer to the caption of Fig. 3.6 for the initial wave parameters.

Figure 3 . 8 :

 38 Figure 3.8: (a) Maximum wave envelope and time evolution at (b) t = 12 and (c) t = 24 for the PNG source with the steepness parameter p 0 = 2 using linear non-dispersive (solid line), linear dispersive (dashed line), non-linear non-dispersive (triangles) and weakly non-linear weakly dispersive (squares) theories. Refer to the caption of Fig. 3.6for the initial wave parameters.

  studied the dynamic generation of long waves from an uplifting or subsiding constant depth region, by analytically solving Laplace's equation in a fluid domain of uniform depth. He considered one-dimensional evolution forced from a uniform bottom displacement, then analyzed axisymmetric uniform bottom displacements of finite size. Hammack employed a Laplace transform in time and a Fourier transform in space for the 1D case and Hankel transform for the 2D case. Depending on the speed of the uplift, Hammack identified three regimes, which he named impulsive, transitional and creeping. While the forced

Figure 4 . 1 :

 41 Figure 4.1: Definition sketch: A radial transect of the fluid domain and the bathymetry in physical coordinates.

  r, θ; ω)e ωt dω. (4.4) ω = c+ia is the non-dimensional transform parameter, and c > 0 is a real constant, large enough for the integrals in Eq. (4.4) to exist. The transformed Eqn. (4

f

  tt (r, θ, t)e -ωt dt .(4.6)

. 9 )

 9 which has a branch-cut on the negative real axis of the complex plane ω introduced by K 0 . Due to the positive real part c of ω = c + ia, arg(ω) ∈ (-π/2, π/2). When |arg(χ)| < π/2, the modified Bessel function of the second kind K 0 (χ) ≈ln(χ) as χ → 0. Boundedness of the free surface elevation at r = 0 requires β 1 = 0. Hence, ζh (r; ω) = α 1 I 0

  (4.8) and is a standard modified Bessel equation of zeroth order. The two independent solutions are the modified Bessel functions I 0 (ωr) and K 0 (ωr). The general solution of the governing equation (4.15) isζ(r; ω) = α 2 I 0 (ωr) + β 2 K 0 (ωr) ,(4.16)with a branch-cut on the negative real axis of the complex plane ω introduced by K 0 . Due to the positive real part c of ω = c + ia, arg(ω) ∈ (-π/2, π/2). When | arg(χ)| < π/2, the modified Bessel function of the first kind I 0 (χ) ≈ e χ / √ 2πχ as χ → |∞|. Therefore, boundedness of the free surface elevation as r → |∞| requires α 2 = 0. Hence ζ(r; ω) = β 2 K 0 (ωr) .

Figure 4 . 2 :

 42 Figure 4.2: Comparison of our solution with predictions from the NSWE solver VOLNA and Hammack's solution when there is no sill (c d = 1), at r = 0 (a -c) and r = 1 (d -f) for different values of the rate of seafloor deformation γ. In the top panel γ = 0.1, which corresponds to the creeping regime, in the middle panel γ = 1 corresponds to the transitional regime and in the bottom panel γ = 10, which corresponds to the impulsive regime.

Figure 4 . 3 :Figure 4 . 4 :

 4344 Figure 4.3: The effect of bathymetry on the wave formation. Time histories of the free surface elevation for different heights of the sill at r = 0 (top panel), r = 1 (middle panel) and r = 5 (bottom) panel. The left column (a -c) corresponds to a bottom motion with γ = 1 , and the right column (d -f) corresponds to a bottom motion with γ = 10.

  areas around the globe. Especially during the last decade two of the most catastrophic tsunamis ever recorded, the December 2004 tsunami in Indonesia [Liu et al., 2005, Titov et al., 2005] and the most recent March 2011 event in Japan [Fritz et al., 2012, Fujii et al., 2011, Ide et al., 2011, Mori et al., 2011] spread panic and pain combined with a huge economic loss at the damaged sites. On the positive side, increased public attention to tsunamis has raised awareness and preparedness, which is the only effective countermeasure and has saved lives, like during the Chilean tsunami in March 2010 [Peachey, 2010]. By better understanding the generation, evolution and run-up of tsunami waves, scientists should ultimately provide early warnings and education to coastal communities. Run-up is defined as the maximum wave uprush on a beach or structure above still water level. Since the 1950's tsunami run-up on a plane beach has been extensively studied by Antuono and Brocchini [2010], Brocchini and Peregrine [1996], Carrier and Greenspan [1958], Didenkulova and Pelinovsky [2008], Keller and Keller [1964], Synolakis [1987], Tadepalli and Synolakis [1994] and Stefanakis et al. [2011] among others. All these studies deal with the mathematical description of long wave run-up on uniform sloping beaches. The catastrophe in Babi

Figure 5 . 1 :

 51 Figure 5.1: Schematic of the geometry of the experimental setup.

  and machine (or statistical) learning learning algorithms have been developed for the "Active Experimental Design", which uses the existing experimental results as a guide for the selection of the future query points, such as Gramacy and Lee[2009]. To achieve that, a statistical model f (x) of the experiment is built, which is constantly updated as new results arrive. Using the predictions of this statistical model (emulator) the future query points are selected according to the objective of the experiment until we can confidentially say that the objective has been achieved. This dynamic approach can further reduce the computational cost. Moreover, building an emulator has further advantages, the most important one being the ability to use it instead of the actual simulator since it is much less computationally demanding to evaluate and thus can be applied very rapidly, especially in cases where someone needs a quick forecast. Depending on the emulator, it also possible to perform a sensitivity analysis of the model output to the several input parameters. A recent and probably the first example of an emulator built in the context of tsunami research is that of[START_REF] Sarri | Statistical emulation of a tsunami model for sensitivity analysis and uncertainty quantification[END_REF], who emulated landslide-generated tsunamis on a plane beach based on the theoretical model of[START_REF] Sammarco | Landslide tsunamis propagating along a plane beach[END_REF].

Figure 5 . 2 :

 52 Figure 5.2: The unstructured triangular grid. Colors represent bathymetric contours. The areas of high grid density on the beach, are the locations of run-up measurements.

(Figure 5 . 3 :

 53 Figure 5.3: Gaussian Process inference of the posterior mean µ (blue line) and variance σ based on 4 realizations (blue crosses). The high confidence region (area in grey) is delimited by f + and f -.

  Fig. 5.3 in a sample unidimensional problem. The posteriors are based on four observations of a Gaussian Process. The vertical height of the grey area is proportional to the posterior deviation at each point.

. 7 ) 1 Figure 5 . 4 :

 7154 Figure5.4: Two queries of GP-UCB-PE on the previous example. The lower confidence bound on the maximum is represented by the horizontal dotted green line at y • t . The relevant region R is shown in light green (without edges). The first query x 0 is the maximizer of f + . We show in dashed line the upper and lower bounds with the update of σ after having selected x 0 . The second query x 1 is the one maximizing the uncertainty inside R.

Figure 5 . 5 :

 55 Figure5.5: Relationship between the simple regret r K t , unknown in a real situation, and the rank correlation ρ Xv (π t-1 , π t ) (in log-scale), for the synthetic function Himmelblau. The stopping threshold ρ 0 was set to 10 -4 , and the lag ℓ to 4, the algorithm stopped at iteration 12, after having found a good candidate at iteration 7 and the true maximum at iteration 11.

(Figure 5 . 6 :

 56 Figure 5.6: Distribution of the final number of iterations T and the final gap (minimum regret) G T for 4 different thresholds.

Figure 5 . 7 :

 57 Figure 5.7: Visualization of the synthetic data sets used for assessment.

Figure 5 . 8 :

 58 Figure5.8: Experiments on several real and synthetic tasks. The curves show the decay of the mean of the simple regret r K t with respect to the iteration t, over 64 experiments. We show with the translucent area the 95% confidence intervals.

Figure 5

 5 Figure5.10: Histogram of the run-up amplification on the beach directly behind the island compared to the run-up on a lateral location on the beach, not directly affected by the presence of the island.

Figure 5

 5 Figure5.13: Run-up amplification (RA) as a function of the wavelength to the island radius (at its base) ratio. The color code indicates the surf similarity (Iribarren number) computed with the beach slope and multiplied with the wave nonlinearity (wave height to water depth ratio).

Figure A. 1 :Figure A. 2 :provided 2 ) 2 ( 6 )Figure C. 2 :Figure C. 3 : 1 ( 1 ) 2 [ 1 -

 12226231121 Figure A.1: Maximum run-up amplification of cnoidal waves as a function of parameter m for an infinite sloping beach (L = 4000m, tan θ = 0.13).

For t > 0 2 (

 02 we take the contour Γ as before: ω) e ωt dω (C.20) We have seen that application of Cauchy's theorem to the first integral gives ∮ Γ ζ(r; ω) e ωt dω = 0. By following the same procedure as we have done for the far field, we can prove that ∫ Γ +,- ζ(r; ω) e ωt dω = lim δ→∞ ∫ ϵ,c c,ϵ ζ(r; s ± iδ) e (s±iδ)t ds = 0 (C.21) Therefore ∫ c+i∞ c-i∞ ζ(r; ω) e ωt dω = ∫ ϵ+i∞ ϵ-i∞ ζ(r; ω) e ωt dω (C.22) with ϵ small. By doing the same rotation as before, ω = e -iπ/2 s, the previous expression becomes ∫ c+i∞ c-i∞ ζ(r; ω) e ωt dω = i ∫ +∞+iϵ -∞+iϵ ζ(r; -is) e -ist ds . = iL 2 (s) (C.23) The new integration path in the limit ϵ → 0 is deformed due to the pole at s = 0 as shown in Fig. C.3, which yields L ; -is) e -ist ds (C.24)

δ 2 e 1 ) 2 .

 212 2iϕ ln(-iδe iϕ ) ] e -iδe iφ dϕ = ; -is) e -ist ds (C.27) D.1 Solution in Transformed Space D.1.1 The near field (r < In the near field h * = h * 1 and therefore h = h * 1 /h * = c d . Furthermore, since both the bathymetry and the forcing are axisymmetric, the θ-term is zero in Eq. (D.2)

) (D. 7 )D. 1 . 2

 712 of variables χ = ωr/ √ c d , the above equation becomes a standard Bessel equation of zeroth order, whose two independent solutions are the two Bessel functions J 0 (χ) and Y 0 (χ). Therefore the general solution to Eq. (D.5) isζh (r; ω) = α 1 J 0 ( ωr √ c d ) + β 1 Y 0 ( ωr √ c d ) (D.6)However, Y 0 (χ) ≈ 2 π ln(χ) as χ → 0. Therefore, boundedness of the free-surface elevation at r = 0 requires β 1 = 0. Henceζh (r; ω) = α 1 J 0 ( ωr √ c dIn order to find the solution of the forced ordinary differential equation, the method of variation of parameters requires the Wronskian of the two independent solutions of the homogeneous equation (Eq. (D.5)). From[START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] we haveW {J 0 (χ); Y 0 (χ)} = 2 πχ (D.8)Finally, the solution of the inhomogeneous ordinary differential equation (Eq. (D.5)) is ζ(r; ω) = α 1 J 0 solution, W (ρ) = 2/(πρ) is the Wronskian of the two homogeneous solutions and ftt (ω) = -iγω/(γiω) is the forcing term. The integration yields P (r; ω) = -ftt (ω) The far field (r > 1)

-c d α 1 ω √ c d J 1 16 ))H 1 ω 2 B 1 ]

 116121 ). However, from all these options, only the Hankel function of the first kind satisfies the Sommerfeld radiation condition and thus the solution to Eq. (D.12) isζ(r; ω) = α 2 H 0 (ωr) (D.13)where H 0 is the Hankel function of the first kind and order zero. The coefficients α 1 and α 2 will be obtained by the matching conditions at r = 1.D.1.3 Matching at r = 1The near and far field solutions are matched at the common boundary r = 1. We require continuity of the free-surface elevation ζ and the radial fluxes h ζr , which yield respectivelyα 1 J 0 ( ω √ c d ) -P (1; ω) = α 2 H 0 (ω) ( ω √ c d ) c d P r (1; ω) = -α 2 ω H 1 (ω) (D.14)By solving the linear system of equations (D.14), we obtain the expressions for α 1 and α 2 (the dependence on c d is omitted for brevity):α 1 (ω) = ω H 1 (ω) P (1; ω) + c d H 0 (ω) P r (1; ω) ω H 1 (ω) J 0 (By replacing the expression for P (r; ω) (D.11) and by noting thatP r (r; ω) = ftt (ω) (ω) J 0 ( ω √ c d ) -√ c d H 0 (ω) J 1 ( F (ω) . (D.19)The expression of the transformed free-surface elevation in the near field can be further simplified: for r < 1 ftt(ω) ω 2 B F (ω) H 0 (ωr) for r > 1 (D.21)

ω 2 [ 1 ]ω 2 B

 212 A F (ω) J 0 ( ωr √ c d) e -iωt dω for r < 1 (D.22) F (ω) H 0 (ωr) e -iωt dω for r > 1 (D.23)

3

  For their computations, they placed a relaxation zone close to the wave generation area. It is applied so that no reflected waves from the beach interact with the forcing boundary because there is no clear understanding of how to impose both incoming and outgoing waves at a boundary. It is convenient because it allows for a reduction of the computational cost and has been used successfully in several other studies [e.g.[START_REF] Lu | Numerical simulation of two-dimensional overtopping against seawalls armored with artificial units in regular waves[END_REF], Madsen et al., 2002, Mayer et al.,

  2.6).Run-up resonance in the laboratory experiments of[START_REF] Ezersky | Physical simulation of resonant wave run-up on a beach[END_REF] was achieved for breaking waves as well, but for these cases they did not comment on the location where the breaking takes place. Hence, probably wave breaking is not the key factor to the resonant mechanism. Long wave breaking in the context of NSWE and its physical demonstration is a subtle issue. As noted by Synolakis[1987], the NSWE tend to predict wave breaking sooner than it actually happens in nature. Still there is an open question about whether tsunamis break when they shoal up a beach.Madsen et al. [2008] suggest that the main flood wave does not break but instead short waves riding on top of the main tsunami do break, giving

Table 5

 5

	.1: Physical parameter ranges
	tan θ i tan θ b d h ω	0.05 -0.2 0.05 -0.2 0 -5000m 100 -1000m 0.01 -0.1rad/s

[START_REF] Briggs | Laboratory experiments of tsunami runup on a circular island[END_REF]

and analytical models [Kânoğlu and Synolakis

  T = O(log T ) defined in Contal et al. [2013]. The factor β t regulates the width of the confidence region.

	Algorithm 1 GP-UCB-PE
	t ← 0 stop ← false while stop = false do
	Compute µ t and σ t with eq. (5.1) and eq. (5.2)
	stop ← stopping_criterion( µ t-1 , µ t ) x 0 t (x) t ← argmax x∈X f + Compute R t with eq. (5.7)
	for k = 1, . . . , K -1 do Compute σ (k) with eq. (5.2) t x k (k) t (x) t ← argmax x∈Rt σ end for
	Query {x k t } k<K t ← t + 1 end while
	with β

  Impact of the batch size K on the gap G K t . The curves show the mean of 64 experiments on the synthetic data Gaussian Mixture, with t init = 10.

	Gap G K t	0.4 0.6 0.8 1						. . . K = 20 . . . K = 10 . . . K = 5 . . . K = 2 . . . K = 1
		0.2
		. . . . . . . . . . . . 0 . . . . . . . .	20	40	60	80	100	120	140
				Number of queries t init + tK
	Figure 5.9:

KFSBC: Kinematic Free Surface Boundary Condition; DFSBC: Dynamic Free Surface Boundary Condition; BBC: Bottom Boundary Condition

Since the theory has been developed over the last five decades, it is useful to provide a short review of the major advances in a condensed form. That way, the alternative approach proposed in the present chapter will appear more clearly.

During run-up, the maximum shoreline velocity is not reached when the wave reaches its maximum run-up. So the joint resonance is not as obvious as it may look. In the literature, there are much less results on velocities than on run-up. One exception is the paper by[START_REF] Madsen | Run-up of tsunamis and long waves in terms of surf-similarity[END_REF].

With h0 = 100m, the corresponding interval for η0 is 0.1m < η0 < 1m.

Even though the length of the computational domain is not a physical parameter, it is of importance from an operational point of view, when one wants to predict run-up elevation based on recorded wave signals at an offshore location. In that case, the incoming wavelength can be either larger or smaller than the beach length.

Note that Tadepalli and Synolakis[1994] gives ĝn(l) = (4εH/π) α

n e -ily 1 [(y1 -y2) l -i(1αnl coth αnl)] cosech αnl, because of the choice of Fourier transform pair.2 Instead, a variant of shallow-water wave theory which captures the effect of dispersion to the lowest order of the Boussinesq approximation[START_REF] Carrier | Tsunami propagation from a finite source[END_REF] can also be used, i.e., ηtt -ηxx -ηyy -1

3 (ηxx + ηyy)tt = 0. It produces almost identical results as the potential flow solution.

Similar observations were made by numerically solving the Poisson equation (3.17) with a dipolar initial disturbance[Dias and Pelinovsky]

By weakly non-linear, we mean the model which retains the lowest-order non-linear terms. Similarly, weakly dispersive refers to the model that considers the lowest-order dispersive terms.

Acknowledgements

Appendix A

More cases of resonant run-up

A.1 Cnoidal waves

Cnoidal wave profiles were also investigated as boundary forcing when L = 4000m and tan θ = 0.13:

)

where sn is a Jacobian elliptic function with parameter m (0 < m < 1), and Appendix C

Contour integral evaluation

Here we will evaluate the complex integrals that appear in the inverse Laplace transform.

C.1 The far field (r > 1)

Consider the integral The resulting governing equation in transformed space is

As before, we will split the fluid domain into two subregions, namely the near field where r * < r * c or r < 1 and the far field where r * > r * c or r > 1. We will solve separately in each subregion the forced long-wave equation (Eq. (D.2)) and then we will match the solutions at the common boundary r = 1.