## Abstract

Tsunamis are rare events with severe consequences. This generates a high demand on accurate simulation results for planning and risk assessment purposes because of the low availability of actual data from historic events. On the other hand, validation of simulation tools becomes very difficult with such a low amount of real-world data. Tsunami phenomena involve a large span of spatial and temporal scales—from ocean basin scales of to local coastal wave interactions of or even , or from resonating wave phenomena with durations of to rupture with time periods of . The scale gap of five orders of magnitude in each dimension makes accurate modelling very demanding, with a number of approaches being taken to work around the impossibility of direct numerical simulations. Along with the mentioned multi-scale characteristic, the tsunami wave has a multitude of different phases, corresponding to different wave regimes and associated equation sets. While in the deep ocean, wave propagation can be approximated relatively accurately by linear shallow-water theory, the transition to a bore or solitary wave train in shelf areas and then into a breaking wave in coastal regions demands appropriate mathematical and numerical treatments. The short duration and unpredictability of tsunami events pose another challenging requirement to tsunami simulation approaches. An accurate forecast is sought within seconds with very limited data available. Thus, efficiency in numerical solution processes and at the same time the consideration of uncertainty play a big role in tsunami modelling applied for forecasting purposes.

## 1. Introduction

After the Sumatra tsunami of 26 December 2004, Bernard *et al.* [1] emphasized that one of the greatest contributions of science to society is to serve it purposefully, as when providing forecasts to allow communities to respond before a disaster strikes.

In order to provide forecasts, a deep knowledge of tsunami run-up and inundation is required. To acquire this knowledge, it is necessary to investigate the geological evidence and impacts of past tsunamis, as well as generation sources. It is often difficult or even impossible to collect data for a tsunami event owing to the fact that tsunamis are rare events. The small number of available observations is the reason for the high uncertainty in the description of the source at the generation phase, even after many years of scientific investigations in this area.

Therefore, computer models, also referred to as simulators, are used to study past and future events. Tsunami simulations can describe past events and model potential hazard scenarios [2]. The computer models can be run to generate evacuation maps that can help in mitigating the hazard. To produce evacuation maps, the maximum possible coastal inundation for a specific location must be computed. Additionally, the computer models can be invaluable to insurance and reinsurance companies.

Simulators are extensively used for analyses of complex physical phenomena, as they can simulate real-world physical processes. Unfortunately, owing to the uncertainties about many tsunami features, it is challenging to simulate and analyse the water motion. The reliability standards for real-time predictions of natural hazards are very high and the development of an accurate model that can represent reality faithfully is very challenging. For the development of the simulators, several choices are required. These choices depend on the available level of expertise and on budget constraints. The latter affects the number of experimental data that can be obtained from laboratory tests. Each of these choices affects the computer model performance on representing reality and hence the amount of uncertainty in the predictions. However, for operational purposes, validation standards are widely accepted [3].

Computer models are used by the Tsunami Warning Centres (TWCs) at locations that are in high danger of tsunami events, for example, in the Pacific Ocean, the Indian Ocean or the Mediterranean (see [1] for a short history of TWCs). Recently opened TWCs in Europe (CENALT in France, KOERI in Turkey and NOA in Greece [4]) are examples of such operations. They help to reduce the risk of future tsunami events by assessing the hazard through early detection and evaluation of tsunamigenic earthquakes, model-based forecasting of wave propagation and coastal impact, and connection to disaster management agencies for evacuation and mitigation measures. Unrealistic estimates and unnecessary false alarms are not desirable, as they can cause needless panic and evacuations that can be costly. A combination of computer model evaluations and historical data of past events is used by the TWCs to create precomputed event scenarios. However, the historical data are usually incomplete and unreliable, rendering the task of early warning very challenging. Precomputed scenario data are combined with observations from the event itself in order to obtain forecasts and alert levels and also to generate inundation maps [5].

## 2. Source modelling

Simulating a tsunami wave can be formulated mathematically as an initial value problem with lateral boundary conditions. While the lateral boundaries are defined by coastlines or open ocean boundaries, where either reflecting, inundating or wave radiation boundary conditions are to be fulfilled, the initial values are less obviously specified.

Tsunamis can have diverse sources and they all require individual treatment. Common sources for tsunami waves are mega-thrust-type earthquakes, submarine landslides, rock slides, volcanic eruptions and meteoroid impacts.

Most commonly studied source mechanisms are earthquake-induced sea surface anomalies, which eventually develop into tsunami waves. Because, in the far field (i.e. more than two average wavelengths away from the source), the exact description of the source mechanism is less sensitive with respect to the inundation and wave height characteristics at the coast, many simplified source parametrizations exist. Two extreme examples of these tsunami initial conditions are shown in figure 1. Simplified analytical descriptions of the shape of the uplift area (e.g. [6,8]) are not common any longer in tsunami simulation. The most commonly used—yet still simplified—source description originates from [9] and is usually applied in the form of the well-known Okada parameter set [10]. This model assumes that the Earth's crust is a homogeneous elastic material. More sophisticated approaches combine a number of sub-plates to form a complex and possibly time-dependent uplift pattern [11]. Finally, high-fidelity models of rupture, with physics-based rate and state friction laws, are at the forefront of complex rupture modelling approaches [7,12,13]. These high-end simulations of the rupture process consume more computing time than the tsunami simulation itself [14].

Mass-movement-induced tsunamis have gained more and more scientific interest. In particular, rock slides and submarine slides are of interest, as their tsunamigenic characteristics can be quite distinct from those of earthquake-induced tsunamis [15]. Locally, they can exhibit extreme run-ups; for example, the highest known run-up of any tsunami occurred after a rock slide in Lituya Bay in Alaska [16]. An overview of simulations with slide-induced tsunamis is given in [17,18].

It is important to note that an accurate approximation of the tsunami source is particularly important in near-field tsunami modelling [19]. An assessment of the uncertainty with respect to the source in the near field is given in [20]. Therefore, studying accurate simulation techniques for the tsunami sources has become an important research direction in recent years.

## 3. Modelling approaches

Modelling tsunami wave propagation and inundation requires a useful description of the physical phenomenon in terms of mathematical equations as well as appropriate methods to solve these equations. In this section, we will introduce a number of approaches to formulate the tsunami wave behaviour as a set of hydrodynamical equations or conservation laws. In order to simplify the description, we will present equations only in Cartesian coordinates, even though a number of formulations use polar or spherical coordinates.

### (a) Linear and nonlinear shallow-water equations

Linear (and nonlinear) shallow-water theory has traditionally been used in tsunami modelling and forecasting. It has been a very successful description of the phenomenon, since up to a water depth corresponding to near-shore regions (say 50 m) the simplifying assumptions of this approximation to the general equations of wave behaviour are well matched. In particular, it is assumed that the wave height *h* is small compared with the mean water depth *H* (*h*≪*H*), and that the characteristic wavelength *L* is large compared with the water depth (*H*≪*L*). In such a regime the vertical components of the acceleration and velocity of water particles are neglected while the horizontal components are independent of water depth.

With these assumptions and a scale analysis the nonlinear set of shallow-water equations can be derived, given by
3.1where *h* is the sea surface elevation (deviation from mean sea surface), *H* the mean depth of the ocean, **v** the two-dimensional depth-averaged horizontal (particle) velocity, *g* the acceleration due to gravity and *S* a source term comprising Coriolis coefficient, bottom friction, eddy viscosity, etc.

Since the wave propagation in the deep ocean is basically represented by a gravity wave behaviour, particle velocities **v** are negligible, and thus the advective (nonlinear) term in the momentum equation (second line in (3.1)) vanishes. Furthermore, as *h*≪*H* the continuity equation (first line in (3.1)) can be written as ∂*h*/∂*t*=−*H*∇⋅**v**. Together, these two modifications form the linear shallow-water equations. Many implementations allow for solving either the linear or the nonlinear equation set by just a parameter switch.

Several established operational codes implement equation set (3.1), such as MOST [21], TUNAMI and NAMI-DANCE both originating from the same source [22], COMCOT the Cornell Multi-Grid Coupled Tsunami Model [23], and the CEA code used at CENALT [24–26], among others. These methods are based on finite-difference mesh discretization of the shallow-water equations (uniform or nested). More flexible mesh alternatives are implemented by finite-element discretizations, e.g. SELFE [27] and TsunAWI [28].

Finite-volume (FV) discretization schemes prove their robustness, flexibility and applicability to the underlying hyperbolic systems of partial differential equations in the VOLNA code [29] and in the adaptive mesh implementation GeoClaw [30]. More specifically, the VOLNA code combines modern numerical techniques for hyperbolic systems with real-world application-oriented design. It can be run efficiently in realistic environments and uses unstructured meshes. For efficiency, the code has recently been ported to graphics processing unit (GPU) architectures. LeVeque *et al.* [31] wrote a useful review on the FV method and adaptive mesh refinement.

A further research code based on Galerkin methods (continuous and discontinuous) is the code TsunaFlash or TAM [32]. It comprises triangular tree-based mesh refinement and is currently further developed in the framework of the project ASCETE [33].

### (b) Boussinesq equations

In recent years, the question whether one needs dispersive wave models or not has been tackled by various groups. A comprehensive paper about this topic is the paper by Glimsdal *et al.* [34]. When the propagation of waves is influenced by frequency dispersion, the long-wave assumption becomes questionable. The simplest model consists in averaging over the water column and keeping only first-order effects of nonlinearity and frequency dispersion. These classical Boussinesq equations are valid for wavelengths down to roughly a few water depths. Several formulations are available. Historically, in the early 1990s, several authors used numerical dispersion to replace physical dispersion [35] in finite-difference schemes. When landslides occur near the water surface, some of the generated waves are relatively short (of the order of the water depth). The classical Boussinesq equations can lead to some numerical oscillations. It is then better to use the generalized Boussinesq equations [36], for which the group velocity of the short waves does not go to zero. There are so many versions of the generalized Boussinesq equations that there is no point trying to write them down in this short review. There are also the fully nonlinear shallow-water generalized Serre equations [37]. They have been rediscovered several times and are in fact used in FUNWAVE [38,39]. Implementations of these dispersive wave models are in-house models or available as standard codes, open-source or commercial. The majority of codes use finite-difference schemes. There is, for example, the Pedersen and Løvholt model, which is designed for long-distance propagation of dispersive tsunamis [40]. With this model, the Japan tsunami of 11 March 2011 can be taken across the Pacific Ocean on a standard desktop during some hours of CPU time. The standard models, such as COULWAVE [41,42] and FUNWAVE [38,39], use more demanding numerical schemes and incorporate a number of effects that are not relevant for oceanic propagation. Hence, simulations of oceanic propagation on single CPUs using these models may therefore be too time-consuming, unless the recent parallel versions (e.g. [43,39]) are used. Recently, Kazolea *et al.* [44] developed a high-order well-balanced unstructured FV scheme on triangular meshes for modelling weakly nonlinear and weakly dispersive water waves over slowly varying bathymetries. The FV scheme numerically solves the conservative form of the equations following the median dual node-centred approach, for both the advective and dispersive parts of the equations.

As emphasized by Glimsdal *et al.* [34], there is no consensus in the tsunami community about the significance of the extra physical features such codes inherit. For seismic tsunamis, the extra cost in terms of both CPU time and memory is not worth it. For tsunamis generated by landslides, and subaerial landslides in particular, more sophisticated wave models are needed, as demonstrated for the potential La Palma tsunami [45,46] and the 1998 Papua New Guinea tsunami [47–49]. For the case of La Palma, the waves are also dispersive in the far field. Long-term propagation of dispersive waves may be approximated by ray (optical) methods [50,51]. If the waves are moderately dispersive, Boussinesq models without the optical approximation are now capable of simulating the far-field tsunami propagation over transoceanic distances (e.g. [52–54]). A review on landslide tsunamis can be found in the present theme issue [55].

### (c) Navier–Stokes equations

The three-dimensional Navier–Stokes equations are the most complete description of fluid wave behaviour and are applied in situations where complex wave interactions are to be studied.

The basic equations can be given by 3.2Additional conservation laws and equations of state complete the system (e.g. [56]). An example for an implementation and application of a Navier–Stokes equation solver in tsunami modelling is given in [46,56], where a three-phase model (air–water–slide) is used to simulate landslide-generated tsunamis.

Another example is the SAGE hydrocode, a multi-material adaptive grid Eulerian code, using a Godunov scheme. This is a commercial code, originally developed by M. Gittings for Science Applications International (SAIC) and Los Alamos National Laboratory (LANL). SAGE has been applied to study asteroid impact-generated tsunamis [57] and landslide-generated tsunamis [45].

## 4. Wave breaking, induced currents, transport of debris

A lot of papers devoted to tsunamis focus on wave run-up. But tsunamis also lead to enhanced currents and damage in harbours and ports. These currents are still poorly understood. As described in Lynett *et al.* [58] and also in [59], observations and modelling of the currents induced by the 2011 Tohoku and 2004 Indian Ocean tsunamis have shown that the strongest flows in harbour basins are governed by horizontally sheared and rotational shallow features, such as jets and large eddies. When examining currents in harbours, one must use modelling approaches that both include the relevant physical processes in the governing equations and use a numerical scheme that does not artificially damp these features. A lack of proper representation of the physics associated with these phenomena may provide drag force estimates that are off by one order of magnitude or more. The implementation of this type of analysis into tsunami hazard studies may then lead to forecasts ranging from an unaffected port to one in which 300 m long container vessels can be detached from their moorings.

Lynett *et al.* [58] showed that the accurate prediction of water surface elevations, which is the present goal of most tsunami hazard mitigation studies, does not necessarily correlate with an accurate prediction of the localized and possibly damaging currents in a port. The model must be able to reproduce the rotational flow patterns as well as the local and turbulence-driven currents in a port or harbour in order to predict damage potential. In other words, quantifying vertical water level rise, buoyancy forces and straight-line current estimates is not sufficient for understanding the hydrodynamic phenomena that can lead to damage in ports and harbours. While the current speeds associated with an unbounded rotational structure are significant, the potential for enhanced currents exists as an eddy approaches a solid boundary or another strong circulation [60]. In the three incidents described during the 2004 Indian Ocean tsunami where the mooring lines of a large ship were broken and the ship drifted freely with the currents in the harbours [61–63], each was located immediately adjacent to a 90° corner in the layout of the container terminal.

Also of interest is the evaluation of tsunami risks offshore in order to formulate guidelines for how far offshore ships must go in order to be safe from strong tsunami currents. Of more recent interest is the interaction of tsunamis with marine renewable energy installations. There were no wave energy converters or current turbines installed at the time of the 11 March 2011 Japan tsunami, but there is some information available about offshore wind turbines, which seem to have survived the tsunami. Until now, there was no interest in the behaviour of tsunamis a few hundreds of metres from the shoreline [64]. This will change in the future.

## 5. Numerical methods

In the previous section, we presented a number of physical description approaches for the deterministic mathematical modelling of tsunami wave behaviour. Except for the study of specialized phenomena, the given equation sets represent simplifications of the full nonlinear hydrodynamic Navier–Stokes equations. Traditionally, the method of finite differences had dominated the field of tsunami simulation [21,22]. However, current numerical methods for solving hyperbolic or weakly parabolic equation sets have been used in recent developments of tsunami simulation software. This section will present a number of current numerical techniques as well as technical developments for improving the efficiency of the solution process.

A number of numerical methods rely on the conservative form of the governing equations, written generically as a flux-form balance equation:
5.1where **U** is the vector of balanced quantities, is the time interval and is the *d*-dimensional spatial model domain, **F** is a flux function and *S* is a function of sources. An example for the typically used nonlinear shallow-water equations with **U**=(*ϕ*,*ϕ***u**)^{⊤}, *ϕ*=*gh* the geopotential height, *g* the acceleration due to gravity, *h* the fluid depth and **u**=(*u*,*v*)^{⊤} the two-dimensional velocity, is
5.2with **I**_{2} the 2×2 identity, **f**_{c} the Coriolis term, *b*=*b*(**x**) the time-independent bottom topography, and *τ*_{b} the bottom friction term. We assume the coordinate **X**=(*t*,**x**)^{⊤}, with the two-dimensional spatial coordinate and the time. Compatible boundary and initial conditions need to be given in order to solve the system.

### (a) Galerkin-type numerical methods

Galerkin-type methods have been introduced to ocean modelling applications for their robustness with respect to complex computational domains [28]. The key idea of Galerkin-type discretization schemes is to replace the continuous function space, in which the solution is sought, by a finite-dimensional (discrete) approximation to that space—usually a piecewise polynomial representation. We briefly outline the process in the following.

The discretization proceeds as follows. We replace the function **U** in (5.1) by linear combinations of piecewise polynomials
5.3where are coefficients and form a polynomial basis for the discrete polynomial function space . Usually, test functions are recruited from as well and we can reformulate (5.1) as
5.4We decompose *Ω* into *M* non-overlapping elements (e.g. triangles) and define *B*_{i} as piecewise orthogonal polynomials:
5.5Integration by parts of the flux term on the elements and replacing the face flux integrals by their numerical approximation **F*** yield the weak form
5.6Repeating the integration by parts leads to the strong form (see [65])
5.7Note that we have derived a method that communicates element information purely by face integrals. If we assume continuous conforming finite-element basis functions, the interior face integrals vanish, and we retain a continuous Galerkin (CG) or finite-element method [66]. If we allow for discontinuous interfaces between elements, we obtain non-zero interior face integrals and a discontinuous Galerkin (DG) method. DG methods enjoy several useful properties: they preserve the local conservation properties of the continuous equations discretely; they are well suited for high-order local basis functions; and they are very well scalable on parallel computing environments, as only nearest-neighbour communication over edges is necessary.

Replacing **U** by its discrete approximation , rearranging the terms, and inverting the mass matrix **M** with leads to a semi-discretized system of ordinary differential equations (ODEs):
5.8This system can then be solved by appropriate ODE solvers, e.g. strong stability-preserving Runge–Kutta schemes. DG methods—and in particular those capable of simulating inundation processes—usually employ slope limiting for stability and positivity preservation [33].

### (b) Finite-volume methods

The FV methods rely again on the flux form of the equations. Their basic assumption is the Eulerian perspective of a fluid dynamics description. In order to discretize the governing equations (e.g. (5.1)) the computational domain *Ω* is approximated by a partition into disjoint FVs *τ*_{e}, *e*=1:*M*, with
5.9The fluid flow is then described by fluxes over the faces of each *τ*_{e}. In particular, the cell integrated conserved quantities (mass, energy, etc.) are balanced by the fluxes over cell faces:
5.10where is a discrete approximation to the cell average of **U**; is an approximation to the time-integrated flux over face *f* and |*τ*_{e}|, *Δt*; |*f*| are the volume/area of cell *τ*_{e}, the time-step length and face area/length respectively; and *ν* is a factor accommodating for the correct scaling due to the definitions of cell mean values. Note that *F*_{f}(*t*) might not be well defined at cell interfaces, as the cell averages of neighbouring cells may jump at the faces. Therefore, *F*_{f}(*t*) is usually replaced by a numerical flux, which is an approximate solution of the Riemann problem at the cell interface (for details of FV method design consult [67] or [68]).

Recent FV implementations for tsunami modelling include the VOLNA code [29]. The adaptive mesh refinement supporting implementation GeoClaw is another example of a modern FV method applied to tsunami modelling [30,31]. The GeoClaw software is described in [69].

### (c) Other methods

The simulation of the violent impact of tsunami-like bores with structures can be performed using traditional grid-based methods, including standard software packages such as Delft3D or OpenFOAM, or particle methods, such as the single-phase, weakly compressible three-dimensional smoothed particle hydrodynamics (SPH) model of [70,71]. In order to avoid large fluctuations in the pressure field and to obtain accurate simulations of the hydrodynamic forces, a Riemann solver-based formulation of the SPH method is used in [70,71]. Time histories of the water surface elevation as well as time histories of the pressure distribution and net total force acting on the structure can be computed. Results show that the magnitude and duration of the impulsive force at initial bore impact depend on the degree of entrapped air in the bore front. Although ensuring a stable pressure field, the Riemann solver-based SPH scheme is believed to induce excessive numerical diffusion, as sudden and large water surface deformations, such as splashing at initial bore impact, are marginally reproduced.

It is well documented in the literature related to breaking wave impacts that entrapped air bubbles can significantly inhibit impulsive shock pressures at the initial impingement of the water with the structure [72].

## 6. Accelerating model performance

Since forecast time is critical in tsunami modelling, several approaches for model performance improvement have been proposed and will be presented briefly here.

### (a) Porting on GPU

One target of current modelling approaches is the real-time assessment of tsunami hazard, where the accurate simulation of a tsunami event requires large numbers of unknowns (corresponding to high-resolution meshes) combined with low times to solution. One powerful computing paradigm consists of general-purpose graphics processing units (GPGPUs) accessible by special programming tools (CUDA or OpenCL). It turns out that some numerical methods can benefit greatly from porting them to GPGPUs, while others are limited by their memory access and computation patterns.

One approach using a high-order ADER (arbitrary high-order derivative) discontinuous Galerkin method for accurate tsunami wave simulation is described in [73]. This code benefits moderately from the GPGPU architecture with a maximum speed-up of approx. 30 over the original CPU code. It turns out that the large number of computations for higher-order polynomial approximations in the ADER-DG method is beneficial for the low memory bandwidth of this type of architecture. On the other hand, the necessary communication of larger numbers of degrees of freedom per element inhibits very large performance gains.

The VOLNA code has also been ported to GPGPUs and demonstrated an increased potential for accelerated performance [74].

### (b) Optimization of single- and multiple-core performance

A more conventional optimization strategy consists of single-core optimization by common software engineering techniques (data locality, restructuring of data layout, kernels, etc.). An impressive performance on single cores is demonstrated for shallow-water equation solvers in [75]. A Sierpinski space-filling curve is used to order the unknowns such that local access patterns can be derived. By this technique, the floating-point units of the processor can almost always be filled with the required data. This is a particularly beneficial optimization strategy for current memory bandwidth-limited hardware architectures.

Parallelization is usually achieved by message passing (MPI) or multi-threading shared memory (OpenMP) programming paradigms. Examples of codes being parallelized by these techniques are the earlier mentioned GeoClaw [30], the operational finite-element code TsunAWI [76], or the Boussinesq solver described in [43]. TsunAWI combines single-core performance optimization with space-filling curves and multi-threading to achieve considerable performance gains.

### (c) Adaptive mesh refinement

While code optimization and parallelization can lead to considerable performance improvements, the large number of unknowns resulting from a scale gap of five orders of magnitude in each dimension can still limit the effectiveness of such measures. Reducing the number of unknowns or even the order of the problem is possible with adaptive mesh refinement techniques [77]. Adaptively refined meshes resolve the tsunami wave interaction locally, in areas of interest, but leave the domain unresolved where no wave interaction is expected. An example is given in figure 2.

These techniques require an automatic and reliable assessment of the mesh refinement criterion. Additionally, an efficient mesh management and data access strategy needs to be employed. Efficient meshing is demonstrated in [78]. An early and successful application of an adaptive mesh refinement code is given in [57], where the commercial SAGE code uses adaptive meshes for studying tsunami hazard generated by asteroid impact. A versatile and well-validated adaptively refined software for tsunami simulations is also given by GeoClaw [30].

## 7. Statistical emulators

A very efficient and viable solution to the problem of generating many computationally expensive simulations is to use a statistical surrogate model, the so-called statistical emulator, in place of the computer model. The statistical emulator gives perfect predictions at the input points that are used in its generation process (it interpolates) and it enables a nonlinear interpolation to predict the model output at any other input point that is within the computational domain. Statistical emulation does not accelerate the model itself.

The significant advantage of using the emulator is that it is much less computationally demanding to be evaluated and, therefore, it can be employed to carry out fast predictions and inexpensive analyses, such as sensitivity and uncertainty analyses. The sensitivity analysis describes how significantly the model outputs are affected by changes in inputs and the uncertainty analysis estimates the uncertainties in the predictions due to uncertain inputs. The use of the emulator in place of the simulator adds an extra source of uncertainty, since it is an approximation. However, the uncertainty induced by using it can be estimated and included into the overall assessment of uncertainty.

Geist [79] studied the effect that rupture complexity has on the local tsunami wavefield. He used a stochastic source model. Løvholt *et al.* [80] studied the variability in coseismic slip along a fault plane and the subsequent tsunami run-up that occurs along a coastline in the near field. Simple linear models for the tsunami generation, propagation and run-up estimation were used. Sarri *et al.* [81] presented an application of Gaussian process statistical emulation for the accurate approximation of a landslide-generated tsunami model. The emulator was built from the combination of selections for regression and residuals covariance functions and parameters with a significantly small number of model evaluations. The resulting emulator was validated using leave-one-out diagnostics, concluding that the emulator's predictions highly agree with the computer model's evaluations. Additionally, the emulator was used further for the computationally demanding sensitivity and uncertainty analyses. Sarri *et al.* [81] demonstrated the potential of reducing significantly the computational cost using the emulator in place of the computer model.

A very recent study on the quantification of uncertainties for a tsunami computer model using another statistical surrogate is performed by Sraj *et al.* [82]. The authors investigate the uncertainties in the resulting wave elevation predictions due to the uncertainty in a specific parameter, which is the Manning's friction coefficient. This coefficient is a parametrization of the effect of bottom friction. Sraj *et al.* [82] use polynomial chaos to build a surrogate model that is a computationally cheap approximation of the computer model. The polynomial chaos surrogate model is used for efficient quantification of uncertainty in the predicted wave elevations, as well as for sensitivity analysis. Following that, the authors apply Bayesian inverse modelling to estimate the posterior distribution of the Manning's coefficient using data collected during a particular tsunami event.

Various sampling techniques have been developed which aim to reduce the number of numerical experiments by finding a small but representative sample in the input space. These techniques are commonly referred to as ‘experimental design’, e.g. [83], 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. In this procedure, all such designed points are queried. This has been a substantial advancement, compared with the regular grids approach, but still involves massive computations.

More recently, adaptive design [84] and machine (or statistical) learning algorithms have been developed for the ‘active experimental design’, also known as sequential design. They use already computed results as a guide to select future query points [85,86]. To achieve this, a statistical model *f*(*x*) of the experiment is built and is constantly updated as new results arrive. Using the predictions of this statistical model (emulator), future query points are selected according to the objective of the experiment, until the objective has been achieved. This dynamic approach further reduces the computational costs and allows investigation of larger parameter ranges. It was recently applied to the problem of run-up amplification behind an island by Stefanakis *et al.* [87].

## 8. Conclusion

Almost 10 years ago Bernard *et al.* [1] wrote a future-looking section, entitled ‘Frontiers for the future’. In that section, there were three items: (i) forecasting tsunami impacts, (ii) human response simulations and (iii) real-time remote sensing. Our review shows that the forecast of tsunami impacts has improved. New unstructured and adaptive grid techniques allow one to bridge the gap of scales involved in tsunami propagation and inundation simulation. Robust numerical methods for accurately representing the wave behaviour in diverse modelling tasks have been successfully applied. Computational acceleration by adaptive mesh refinement, hardware acceleration by GPU utilization and parallelization make real-time tsunami forecast technically feasible. First approaches to handle uncertainty quantitatively have been published. Thus, the vision of tsunami forecasting—(i) above—has become almost a reality. However, since the much more expensive deployment of sensors and measurement capabilities mentioned under (iii) has not kept up, a real-time tsunami forecast is only possible in the far field [88]. In the near field, an accurate assessment of the source is still not feasible. This means that mitigation measures for this case still have to rely on more traditional approaches and a well-educated population.

## Authors' contributions

J.B. and F.D. planned the review paper. They contributed equally to the writing and review of the final manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

The authors gratefully acknowledge support through project ASTARTE (Assessment, Strategy And Risk Reduction for Tsunamis in Europe) FP7-ENV2013 6.4-3, grant no. 603839. J.B.'s work was further supported by Volkswagen Foundation through project ASCETE (Advanced Simulation of Coupled Earthquake Tsunami Events) grant no. AZ 85434.

## Acknowledgements

The authors thank Costas Synolakis, Finn Løvholt, Audrey Gailler, Serge Guillas and Randall LeVeque for their careful reading of the manuscript.

## Footnotes

One contribution of 14 to a theme issue ‘Tsunamis: bridging science, engineering and society’.

- Accepted July 29, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.