The UBO-TSUFD tsunami inundation model : validation and application to a tsunami case study focused on the city of Catania , Italy

Nowadays numerical models are a powerful tool in tsunami research since they can be used (i) to reconstruct modern and historical events, (ii) to cast new light on tsunami sources by inverting tsunami data and observations, (iii) to build scenarios in the frame of tsunami mitigation plans, and (iv) to produce forecasts of tsunami impact and inundation in systems of early warning. In parallel with the general recognition of the importance of numerical tsunami simulations, the demand has grown for reliable tsunami codes, validated through tests agreed upon by the tsunami community. This paper presents the tsunami code UBO-TSUFD that has been developed at the University of Bologna, Italy, and that solves the non-linear shallow water (NSW) equations in a Cartesian frame, with inclusion of bottom friction and exclusion of the Coriolis force, by means of a leapfrog (LF) finite-difference scheme on a staggered grid and that accounts for moving boundaries to compute sea inundation and withdrawal at the coast. Results of UBO-TSUFD applied to four classical benchmark problems are shown: two benchmarks are based on analytical solutions, one on a plane wave propagating on a flat channel with a constant slope beach; and one on a laboratory experiment. The code is proven to perform very satisfactorily since it reproduces quite well the benchmark theoretical and experimental data. Further, the code is applied to a realistic tsunami case: a scenario of a tsunami threatening the coasts of eastern Sicily, Italy, is defined and discussed based on the historical tsunami of 11 January 1693, i.e. one of the most severe events in the Italian history.


Introduction
Tsunamis are long ocean waves with horizontal displacement of the fluid particles predominant on the vertical displacement and with motion nearly uniform along the vertical column from the sea surface down to the sea floor.Most of the current tsunami models, such as the non-linear shallow water (NSW) and the Boussinesq approximations, make use of depth-averaged water-wave equations to compute tsunami propagation across the ocean basins and tsunami evolution in the near-and on-shore areas.In this paper we restrict our attention only to NSW equations that have been long and widely investigated by the tsunami scientific community.One of the first and significant efforts to study the impact of long waves against the coasts was the analytical study by Carrier and Greenspan (1958) who explored the dynamics and run-up of long ocean waves climbing up a constantslope beach.Their paper can be considered as a milestone in tsunami analytical studies and inspired several works in the following years (e.g.Peregrine, 1967;Thacker, 1981;Synolakis, 1987;Tadepalli and Synolakis, 1994;Li and Raichlen, 2001;Carrier et al., 2003;Kanoglu et al., 2004;Tinti and Tonini, 2005) that had also the purpose to provide analytical reference solutions for numerical models.An interesting and exhaustive review of analytical methods and solutions that have been developed in the frame of tsunami science was prepared some years ago by Synolakis and Bernard (2006).In parallel with theoretical studies, laboratory experiments were carried out with the aim to validate theories and to confirm analytical achievements, but also to enlighten new S. Tinti and R. Tonini: The UBO-TSUFD tsunami inundation model phenomena and to stimulate new analytical and numerical analyses (e.g.Hammack, 1972;Synolakis, 1987;Grilli et al., 1994;Liu et al., 1995;Fritz et al., 2003;Sue et al., 2010).
Despite the fundamental role of analytical studies, they have the limitation to provide explicit solutions only for ideal waveforms and very simple basin geometries.The need of considering more realistic scenarios, together with the impressive and rapid technological advancement of computational tools, has led to a fast development of numerical methods to solve tsunami equations and to produce tsunami simulations.Among the codes that are more extensively used today, one can mention COMCOT (Liu, 1994), TUNAMI-N2 (Imamura, 1996) and MOST (Titov and Synolakis, 1998) that all solve NSW equations through a finite-difference technique.Other more recent NSW solvers are GeoClaw (George and LeVeque, 2006), formerly known as TsunamiClaw, that applies an adaptive finite-volume method and is under continuous development (see e.g.LeVeque et al., 2011), and the parallelized finite-difference model developed by Nicolsky et al. (2010).
Numerical models need testing and validation, which seems a trivial statement against which nobody could raise sensible objections.However, in case of codes that have large application and relevant social impact such as tsunami simulation packages that might be used to elaborate scenarios in the context of tsunami early warning systems, testing and validation become a rather sensitive issue, which requires the specification of agreed procedures of accreditation.One of the first attempts of addressing this topic was held at the second workshop on long wave run-up organised at Friday Harbor, San Juan Island, Washington, in 1995 under the sponsorship of the Geo-Hazard program of the US National Science Foundation and reported in Yeh et al. (1996).Here a series of benchmark cases was proposed to validate the flooding prediction of tsunami numerical codes.Years later, after the idea of benchmarking models against common problems gained consensus, in the third workshop of the series held in Catalina Island, California, in 2004, four benchmark cases were selected and proposed (see Liu et al., 2008) that since then have been considered standard validation tests in the tsunami modellers community.In later papers, Synolakis et al. (2007Synolakis et al. ( , 2008) ) stressed once more the importance of validating/verifying a numerical model before using it for tsunami mitigation plans or for tsunami warning purposes to avoid unrealistic predictions that could cost a high death toll or create unjustified panic.It is, however, important to point out that good results of a model in a series of benchmark tests do not ensure obtaining correct inundation results in a realistic case, since the prediction success depends on a number of additional factors (for instance, the availability of adequate topo-bathymetric datasets), but they no doubt give more reliability to the model predictions.In other terms, success in benchmark is a necessary, but not sufficient, condition for a model to provide credible computational results.
The goal of this paper is to present the NSW finitedifference tsunami model that was developed and is maintained by the Tsunami Research Team of the University of Bologna and that is called UBO-TSUFD.The code has been already extensively applied to several cases of tsunami propagation and inundation, especially in the frame of the studies on European tsunamis carried out in the FP6 EU-funded projects TRANSFER (http://www.transferproject.eu/)and SCHEMA (http://www.schemaproject.org).Furthermore, it has been used to study the case of the 2009 Samoa tsunami in the Pacific Ocean and the implication on tsunami forecasting strategies in the perspective of tsunami early warning systems (Tonini et al., 2011a).This is the first time the code is presented in a scientific journal.The work is mainly focused on the validation of the propagation and the run-up calculation of the code by illustrating its performance in four benchmark cases; two of which are taken from the set of the Catalina Island benchmark problems quoted above.In the next section a brief introduction of the NSW theory will be presented, followed in Sect. 3 by the description of the numerical scheme as well as the run-up algorithm used in the UBO-TSUFD code.Section 4 is dedicated to the benchmark exercises.In Sect. 5 the code is applied to study one of the most important tsunamis of the Italian history (see the Italian Tsunami Catalogue, Tinti et al., 2004), namely the tsunami following the 11 January 1693 earthquake that occurred in eastern Sicily, Italy, and that flooded several towns.In the present work, we will concentrate our attention on the area of the Catania harbour and of the beach south of the town, called La Plaia, which seems to be particularly exposed to tsunami attacks.Finally, conclusions regarding the efficiency and future improvements of the model will be drawn.

Background theory
The NSW theory is an approximation of the more general set of Navier-Stokes equations, which govern the motion of all fluids.The NSW equations are non-dispersive (since the phase speed of the waves, being proportional to the square root of the ocean depth, does not depend on the wave frequency or wavelength) and hold under the assumptions of fluid incompressibility and of pressure hydrostaticity.Hence, the vertical component of fluid particles acceleration is negligible compared to the gravity acceleration and the horizontal components of the velocities are constant along the water column (see Stoker, 1957).This reduces the problem from three to two dimensions in space and one dimension in time.In the Cartesian plane (x, y), the NSW equations can be given the following form: (1) where x and y are the horizontal coordinates, u and v are the depth-averaged velocity components, respectively, along the x and y directions, g is the gravity acceleration, h is the undisturbed water depth, η is the water surface elevation measured from the still sea surface, f x and f y are the x and y components of the bottom friction (see Fig. 1).The above formulation does not take into account the Coriolis acceleration and therefore, strictly speaking, is adequate only for tsunami propagation over short distances and with short duration, more specifically for simulations covering a few hours of tsunami travel times over a range of hundreds or at most a few thousands of kilometres.We neglect Coriolis terms in the UBO-TSUFD version of the code presented in this paper, since they are not required in any of the treated benchmark problems and not even for the simulation of the 1693 tsunami off eastern Sicily.Indeed, in the Mediterranean all the main tsunami effects occur within a time interval too short for the Coriolis terms to be effective.
In the Eqs.( 1)-( 3), it is very common to introduce the discharge fluxes M and N (Goto et al., 1997) that are related to the horizontal velocities u and v according to the definitions: where D = h + η is the total water column.It follows that Eqs. ( 1)-( 3) can be re-formulated as In tsunami models the bottom friction components f x and f y are usually expressed in terms of the fluxes M, N and of a friction coefficient C f , which in turn can be given in terms of the Manning's roughness coefficient n (see Satake and Tanioka, 1997;Dao and Tkalich, 2007), the final result being Here typically the value of n ranges between 0.01 and 0.06 m −1/3 s, depending on the nature of the bottom.It is worth stressing that three of the benchmark problems that have been used to validate the code and that will be treated in Sects.4.1, 4.2 and 4.3 are inviscid, since they are based on analytical solutions or ideal conditions.Therefore, they are no use to test the part of the code handling the bottom friction.This was tested only through a benchmark case that is based on a laboratory experiment (see Sect. 4.4).

Description and validation of the code
The code UBO-TSUFD solves both linear and non-linear shallow water equations by means of a leapfrog (LF) finitedifference scheme on staggered grids and calculates tsunami inundation caused by both seismic and landslide sources.It also implements a nested grid algorithm, which allows one to increase resolution and to calculate detailed inundation maps only in specific target areas.Neither of such capabilities, i.e. treatment of landslide-induced tsunamis and grid nesting, will be illustrated in the present paper, since their validation and application will be shown in subsequent papers.All benchmarks considered here are specific to validate the propagation of tsunamis starting from given initial fields of elevations and discharge fluxes (Sects. 4.1,4.2 and 4.3) or from a forcing prescribed on the boundary (Sect.4.4).All cases have been performed on a single computational grid.The code has been developed in standard FORTRAN 90/95 and can be compiled by open source compilers, i.e. gfortran (http://gcc.gnu.org/wiki/GFortran), which is freely available for the most popular operating systems.

Numerical scheme
The set of Eqs. ( 6)-( 8) with the additional expressions ( 9)-(10) describes a partial differential equations problem of hyperbolic type.The code UBO-TSUFD uses a finite difference method, where the domain is discretised in a regular grid with a finite number of nodes spaced x and y, and likewise the time axis is discretised in regular steps t, and the derivatives are replaced by differences over small intervals.The accuracy of the method depends on the density of points considered and on the truncation error of the computation.The scheme adopted is a LF method over staggered grids in space and time, meaning that the fluxes M and N are calculated over a grid that is shifted with respect to the grid used to compute η by a grid half-step in all the variables, i.e. along x, y and t, respectively (see Fig. 2).Restricting to the non-advective and inviscid terms, Eqs. ( 6)-( 8) are expressed as follows: where These formulas are explicit and all the quantities can be calculated using the information of the previous time step, identified by the index k, while i and j refer to the x and y coordinates.Half steps, both temporal and spatial, are de-fined by semi-integer indexes.Advective terms are presented and described in Appendix A.
The numerical computations are performed in the order indicated by the previous equations.At the time step k, firstly one calculates the discharge fluxes M k+1/2 i,j and N k+1/2 i,j and then the sea surface elevation η k+1 i,j by using the values available at the previous time step.In case of a tsunami induced by an earthquake, the initial conditions assigned for t = 0 are expressed by η(x, y, 0) = η 0 (x, y) and M(x, y, 0) = N (x, y, 0) = 0, where η 0 (x, y) is related to the co-seismic vertical displacement of the sea floor.In terms of discretised variables, the initial conditions can be written as . Further details on this numerical scheme can be found in Goto et al. (1997).

Stability condition
The LF scheme is known to be stable according to the von Neumann stability criterion and the Courant-Friedrichs-Lewy (CFL) condition if x/ t > c max and if x/ t > c max , where c max is the maximum phase velocity of the tsunami wave on the grid (see Smith, 1985).This ensures that tsunami speed is always smaller than the speed associated with the grid, which is the speed required to cross one cell in the x or y direction within a single time step.Since c max = √ gD max , in principle c max depends on the solution itself and cannot be determined a priori.In practice, however, the value D max = (h + η) max is dominated by the term h in deep ocean and does not change sensibly during the tsunami computations.In case of grids with variable steps the above inequalities may be better written as min( x, y)/ t > c max .Usually it is quite easy to select a value of t small enough that the stability condition is fulfilled over the whole space domain and over the whole time interval used in the calculations.

Boundary conditions
At the boundaries of the computational domain there is the need to set conditions prescribing specific behaviour for the wave fields.In tsunami propagation problems one can pose a condition of full reflectivity, which is equivalent to limiting the domain by means of an infinitely high vertical wall and is useful to handle waves travelling in closed basins and in laboratory tanks.Alternatively, a condition of free transmission can be posed, which allows waves to go out of the domain through a totally transparent boundary.In the code UBO-TSUFD boundary conditions are applied to the discharged fluxes in the nodes that are placed in the right (east) and the upper (north) sides of the boundary cells, implying a geometric asymmetry of the computational grid which must be accounted for.For the sake of clarity, let us suppose that a domain is covered by a rectangular mesh with N x × N y cells, which are numbered with the indexes i = 1, 2, . . ., N x and j = 1, 2, . . ., N y .Let us further assume that the cell (i, j ) has the central node with coordinates x i and y j .According to Fig. 2 and the concept of grids staggered in space, the sea surface elevation η(i, j ) is computed in the centre of the cell, while the x component of the flux M is calculated in the midpoints of the left and right sides of the cell, i.e. in the points (x i ± 0.5 x, y j ).Taking this into account, all the cells of the first column of the mesh, i.e. cells (1, j ; j = 1, 2, . . ., N y ), are considered to be outside of the computational domain, except their right sides, with coordinates (x 1 + 0.5 x, y j ; j = 1, 2, . . ., N y ), that altogether form the left boundary of the grid.Instead, all the cells of the last column, i.e. cells (N x , j ; j = 1, 2, . . ., N y ), are an integral part of the domain, and their right sides constitute the right side of the grid.Bearing this in mind, the boundary conditions for the considered grid can be given through the following formulas: 1+1/2,j = 0, vertical wall on the left side ( 16) N x +1/2,j = 0, vertical wall on the right side (17) , open on the right side, ( where is the local phase velocity depending on the sea surface elevation and therefore changing with time.In the linear approximation, expressions (20) simplify to c 2,j = gh 2,j and c N x ,j = gh N x ,j .
Analogous conditions can be imposed on the lower and upper boundary of the mesh involving the y component of the flux N.
The above scheme can be easily generalised to domains of arbitrary shape.If the grid considered before contains a basin covered only by a subset of the grid cells, and if the right side of cell (i, j ) happens to form part of the left boundary of the basin (i.e.cell (i, j ) does not belong to the basin, but cell (i +1, j ) does), hence Eqs. ( 16) and ( 18) can be re-written as M k+1/2 i+1/2,j = 0, vertical wall on the left side (21) Likewise, also formulas ( 17) and ( 19) can be adapted to pose conditions on the right boundary of an irregular basin.
All expressions given before are adequate to describe tsunami propagation where the tsunami energy is initially already inside the basin, and either prevent the energy to flow out of the boundaries (total reflection) or permit it to pass completely across (free transmission).The case of inputting energy into the basin by forcing a wave through the boundary can be treated quite easily by prescribing the flux discharge as a function of time on the boundary cells.Let us revert to the case of the regular rectangular basin introduced before, and let us consider a forcing on its left boundary.Forcing can be originally expressed in terms of sea surface elevation, or of horizontal velocity or of discharge flux, but must be converted to discharge flux to fit the discretisation scheme.Usually physical forcing is applied to open, rather than closed, boundaries, so that if the forcing discharge flux on the cell (1, j ) at the time t k+1/2 is known and is denoted by F k+1/2 M 1+1/2,j , then the boundary condition on the left side transforms into In case that the sea surface elevation H k 1+1/2,j is prescribed, then the forcing flux can be easily deduced by the aid of the definition (4) under the assumption that the wave enters the basin from the left, namely Forcing from the right boundary as well as forcing from the lower and upper boundaries of the grid can be easily treated by properly adapting the previous formulas.

Run-up
Accurate run-up calculations are a very important objective for any tsunami model to compute inundation maps for hazard and risk assessment purposes, and testing this capability was one of the main reasons leading to the introduction and application of benchmarks (Liu et al., 2008).Indeed tsunami models of the first generation until 1990s were able to produce quite consistent results with regards to tsunami propagation, but often were found to be discrepant and unsatisfactorily inconsistent when computing flooding and related quantities like flow depth, run-up heights, penetration distances, etc.
The algorithm used to solve this issue is based on the dynamic identification of the dry and wet regions within the grid that is on the identification of the instantaneous shoreline position at any time step, and on moving the shoreline according to the local values of the sea surface elevation and of the discharge fluxes.Specifically, the discriminating condition for a cell to be dry or wet is the depth of the cell water column D. If D > D TH , where D TH is a pre-assigned S. Tinti and R. Tonini: The UBO-TSUFD tsunami inundation model threshold value, the cell is assumed to be wet, otherwise it is dry.In the tsunami computation cycle this decision is taken at any time step after calculating the sea water elevation η k+1 i,j by means of the Eq. ( 13).If the cell is found to be dry, then D k+1 i,j is posed equal to zero, and consequently η k+1 i,j = −h i,j .A further decision is taken according to the values of the discharge fluxes in the neighbouring cells.Let us consider the motion along the x direction and let us assume that the cell (i, j ) is found to be dry.Then let us take into account the neighbour two cells (i + 1, j ) and (i + 2, j ) on the right side.If also the cell (i + 2, j ) is dry, hence the cell (i + 1, j ) is forced to be dry, irrespective of its calculated value, since isolated wet cells are assumed to be the result of numerical noise and are not significant from a physical point of view.If the cell (i + 1, j ) is dry, either independent or as the effect of the previous forcing, then the flux on the cell side that is shared with the cell (i, j ) is equalised to zero, i.e.M k+1/2 i,j = 0.If both cells (i +1, j ) and (i +2, j ) are wet, then the boundary between the cell (i, j ) and (i + 1, j ) is assumed to belong to the instantaneous shoreline and is given the same horizontal velocity u as the velocity found on the boundary between the cells (i + 1, j ) and (i + 2, j ), i.e. u(x i+1/2 , y j , t k+1/2 ) = u(x i+3/2 , y j , t k+1/2 ).In terms of the discharged fluxes, this is expressed by Analogous computations are carried out for the cells on the left of the cell (i, j ), namely for cells (i − 1, j ) and (i − 2, j ), affecting the values of M k+1/2 i−1/2,j , as well for the direction y, in such a way to also find proper values for the fluxes N k+1/2 i,j −1/2 and N k+1/2 i,j +1/2 .

Validation of the code through benchmarks
The four benchmark problems selected to validate the code UBO-TSUFD are two analytical problems (Sects.4.1 and 4.2), a specific case developed to analyse the mass conservation (Sect.4.3) and one example taken from a laboratory experiment (Sect.4.4).Benchmarks 1 and 3 are taken from the set proposed in the Catalina Island 2004 workshop (Liu et al., 2008), though the former has been modified to make it suitable for comparison with the outcome of numerical models, as will be explained later on.This benchmark is a one-dimensional NSW analytical solution of an initial value problem, where a wave of given initial profile and initial null velocity propagates over a constant-slope beach.The second analytical solution is a very special two-dimensional problem, where the sea surface of a basin of paraboloidal shape is a plane that oscillates in a periodic way along one of the axes of the basin, entailing that one of the two fluxes, say N(x, y), is identically equal to zero.Here the initial con-ditions are given by specifying both the sea surface elevation η(x, y, 0) and the non-zero flux discharge M(x, y, 0).The third benchmark regards a volume conservation analysis that is applied to the simple case of a plane wave propagating in a flat ocean, as well as to the two benchmark cases mentioned above.The fourth benchmark is a laboratory experiment of a three-dimensional replication of the 1993 Hokkaido-Nansei-Oki tsunami attacking the coast of Monai, Okushiri, Japan.The basin was a tank delimited by two vertical walls where the bathymetry off Monai and the nearby coastal topography including the islet of Muen, located very close to the coast, were reproduced.A wavemaker placed on the tank side facing the Monai coast was used to produce a wave with specified sea-elevation time history.In the following subsections, the results of the benchmark exercises will be illustrated in detail.

Benchmark n. 1: tsunami run-up onto a sloping plane beach
This benchmark problem is given in the paper by Liu et al. (2008) where one can find all the links to the datasets relevant to conduct the test.The case is representative of a leading-depression N wave with initial profile depicted in the top panel of Fig. 3 that travels on a uniform incline with beach slope fixed at 1/10.The solution to the problem is given through analytical integral formulas (see Carrier et al., 2003) and can be calculated to a very high accuracy by means of numerical codes of integral computation, so that it can be considered as known exactly.For the benchmark, free surface and velocity profiles are given at t = 160 s, 175 s and 220 s up to a distance of 50 km from the shore, and shoreline position and velocity are given in the range from 0 to 360 s (one can find all the available data at http://isec.nacse.org/workshop/2004cornell/bmark1.html).
The benchmark treated here is not precisely the same as the original Catalina Island case, since in benchmarking conducted in the frame of the TRANSFER project, it was found that this case leads to a meaningless physical wave (see Fig. 4) with the appearance of a multi-valued solution in the stage when the sea retreat phase ends and shoreline velocity reverts from seaward to shoreward.In order to use a fully meaningful case and to keep a strong resemblance with a benchmark that in the meantime has become a classical one, a slightly modified version of the Catalina Island case was proposed in the TRANSFER project by reducing the initial wave amplitude by 10 % and by re-computing all benchmark curves.The proposed case exhibits a very similar dynamics, with a sharp passage from the lowest sea withdrawal to sea flooding, but curves remain always single-valued.This is the benchmark case adopted here.The code UBO-TSUFD is intrinsically 2-D, and the 1-D case can be solved by imposing initial conditions with all variables (η, h, M) uniform along the y direction, implying that N is identically equal to zero.Three simulations have been carried out by us- We observe that the coarsest resolution ( x = 50 m) is the only one suggested in the Catalina Island benchmark exercise, but we have used also higher resolution to test the convergence of the code.Offshore propagation is very much the same for all cases (Fig. 3), but relevant differences can be found for the shoreline position and velocity prediction (see Figs. 5  and 6), where non-linear factors become relevant.In particular, the 50 m grid evidently has a worse performance with respect to the other two numerical simulations that, on the other hand, result in being very similar to each other.This shows clearly that the numerical solution converges as the discretisation step is made smaller and smaller.This issue can be addressed more quantitatively by numerically estimating the difference between analytical and numerical coastline evolution (both position and velocity).The analytical curve has been re-sampled using the time steps listed above corresponding to the three different grid resolutions.The difference between the analytical and the numerical solutions at each time step i is where θ T i and θ C i are the theoretical and the calculated value, respectively.In Fig. 6 (bottom row) the mobile average of R i is plotted for each grid step, and the mean ( R) and the root mean square (RMS) of R i are listed in Table 1.
The code UBO-TSUFD can determine the position of the shoreline not better than the resolution permitted by the grid discretisation.In the case of this 1-D benchmark, the shoreline can change in time only by fixed steps x.In order to better compare the numerical and the analytical position of the shoreline, a specific post-processing routine has been developed that computes coastline position by extrapolation, as the intersection of the sea surface profile close to the coast and the sea floor profile.If at a given time t k , cell i is the leftmost wet cell of the grid, implying that the shoreline is in a position to the left of the cell centre x i , then the coordinate of the coastline x c (t) can be obtained simply as This formula is meaningful and consistent only if i , a condition that happened to be always fulfilled for the specific case of this benchmark, and therefore no more sophisticated algorithm has been elaborated for this purpose.The velocity of the coastline is computed by numerical differentiation of the shoreline curve x c (t) after smoothing.Figure 6 depicts the coastline position (left panel) and the shoreline velocity (right panel) vs. time.It can be seen that the fit is quite good with the exception of the phase of transition between the retreat and the advancement of the sea when numerical velocity curves tend to be less sharp than the analytical solution, with fit improving as the space step decreases.The case studied here is interesting not only for benchmark purposes, but also as an instructive example showing how a wave with a leading depression and close to breaking attacks the coast.Figure 7 serves this purpose, since it displays the computed sea surface elevations in a sequence of nodes placed close to the area affected by sea regression and coastal inundation.Since the slope inclination β is 1/10, the position of each tide-gauge given in the figure caption can be scaled easily to the elevation of the sea-floor.All computed curves (referring to the case x = 10 m) are quite regular proving that the code UBO-TSUFD performs well even in the critical area of changing from wet to dry and/or vice versa (when the curve is at a constant level, this means that the corresponding node is dry).Notice that the minimum water level, i.e. induced by the arrival of the leading trough, occurs earlier in nodes that are more offshore, which is expected from a trough travelling shoreward.Notice further that the transit from sea retreat to sea inundation is quite quick and that, surprisingly, records attain the maximum water level in reverse order, i.e. the highest values are reached first in the nodes most inland.This shows that when the crest arrives at the coast, the shoreward propagation is not the only key to interpret the wave behaviour, but that also reflection against the coast enters into play and seems to be dominant.

Benchmark n. 2: tsunami run-up of a planar surface oscillating in a paraboloidal basin
Some exact solutions to NSW equations were found by Thacker (1981) for a three-dimensional paraboloidal basin and for particular classes of initial conditions.Here we consider solutions for a paraboloid with a still sea surface being an ellipse of semi-axes L and l.With no loss of generality, the ellipse axes can be taken to lie along the axis x and y.
The sea surface is a plane that oscillates around one of the two horizontal axes and implies that, if the oscillation occurs along one axis, then the component of the velocity along the other axis is null.The exact solutions can be expressed by the following analytical expressions in terms of the sea elevation η and the velocity components u and v: The set of Eqs. ( 27) and ( 28) represent the plane oscillation along the x-axis and along y-axis, respectively.
The basin topo-bathymetry is given by the equation: where h is positive downward.Here H 0 is the depth of the vertex of the elliptical paraboloid, i.e. the deepest point of the basin, A is a constant which determines the amplitude of the motion, and ω 1 and ω 2 are wave frequencies in the x and y directions, respectively.The range of validity of the solutions ( 27) and ( 28) is given by The coastline is the line satisfying the condition D(x, y, t) = 0, and, when ω 1 t = (2p +1)π/2 (p being an integer number), the sea surface elevation η is identically zero, so that the coastline lies on the horizontal plane h(x, y) = 0 and coincides with the ellipse x considered the surface of the basin in conditions of no motion.Notice that the solution is harmonic with frequency ω 1 (ω 2 ) in the discharge flux M (N) and that all the particles move synchronously, but it is not harmonic in the elevation η, that, due to the non linearity of the problem, shows the additional frequencies 2ω 1 (2ω 2 ) and zero, through the term cos 2 (ω 1 t) (cos 2 (ω 2 t)).
For the UBO-TSUFD simulations, the parameters above have been set as follows: A = 235 m, L = 4700 m, l = 1300 m and H 0 = 201.42m, and the resulting basic periods turn out to be T 1 = 470 s and T 2 = 130 s.The grid has been built with x = 10 m, y = 10 m, and the selected time step, respecting the CFL conditions, is t = 0.1 s.The initial con- ditions in input to the UBO-TSUFD model have been derived from formulas ( 27) and ( 28) by computing the values of η 0 i,j , M −1/2 i+1/2,j and N −1/2 i,j +1/2 .Both cases have been simulated, but only the results of the one that is described by the Eqs.( 27) and (29) will be shown.In Fig. 8 the geometry of the basin and the positions of a number of tide gauges are given.The comparison between the numerical and analytical solutions is made by means of Figs.9-11.The results obtained are in very good agreement with the analytical solution.In Fig. 11 26).This point of the coastline moves on the x-axis around the point of the still-water shoreline with coordinates (x = −L = −4700 m, y = 0 m).It is important to note that in the numerical simulation the velocity along the y-axis does not remain equal to zero, but is affected by a numerical noise that grows together with time.This is also illustrated in Fig. 12 where the water elevation (left column) at time instants t = T 1 /4 and t = 3T 1 /4 is not perfectly planar (η = 0) and the N flux is not equal to zero after t = 0.

Benchmark n. 3: volume conservation analysis
One of the basic steps in tsunami modelling validation is to verify how the model deals with the mass conservation.This can be done by checking the variation of the total water volume V (t) with respect to its initial value V (0) at t = 0 (the density of the water is considered uniform and equal to 1).The volume variation index can be defined as follows: Synolakis et al. ( 2007) and Synolakis et al. (2008) stated that the volume variation should be calculated only for the displaced water depth η, because using the total water depth η + h could mask errors.Moreover, the authors stated that a reasonable difference between initial and final displaced volume V D (t) should be less than the threshold of 5 %.However, this criterion could be of hard applicability in cases where the initial displaced volume is equal or very close to zero, because the effect of any small volume variation would lead to instabilities in formula (31).
In order to avoid this issue, we define the following index: where the difference between V (t) and V (0) is normalized to the maximum value of the displaced volume V D (t).To study the behaviour of the code UBO-TSUFD with respect to mass conservation, a specific benchmark was developed that consists of a positive plane wave (maximum height of 2 m) which propagates in a long flat channel that ends on a constant-slope beach.The geometrical sketch of the channel and the offshore propagation are shown in Fig. 13.The wave splits into two fronts of half amplitude which propagate in opposite direction, one moving offshore and one moving onshore.Two different spatial steps (5.0 m and 2.5 m) have been used.The difference between the results of two grids are shown in Fig. 14, where inundation and retreat phases are shown more in details.This is a 1-D case, so formulas ( 31) and ( 32) have been calculated by using section areas A(t) instead of volumes.The balance is considered during the time span needed by the wave travelling shoreward to be completely reflected.The variations of the displaced and total area have been calculated for both grid step resolutions and are plotted vs. time in the bottom panels of Fig. 15.It can be noted that almost no area changes are found while waves travel along the flat section of the channel (for the first 10 s), whereas they are found later during the inundation/withdrawal phase.It can be further noted that the finer grid gives area discrepancies much better than the coarser one (though the computed wave profiles of the two cases portrayed in Fig. 14 seem very alike); however, these remain well below 1 %.
The same analysis is proposed for the simulations conducted for benchmarks 1 and 2 with results given in the panels, respectively, of the first and second row of Fig. 15.As for benchmark 1, it is seen that passing from the coarsest (50 m) to the finest grid (5 m), one gets not only a significant improvement of the solution (as discussed in Sect.4.1), but also a remarkable reduction of the area error that remains well below 1 ‰.It is confirmed that our model conserves area/volume while handling pure offshore propagation (the error curves remain close to 0 in the first 100 s of simulation for all the grids used) and that errors occurs while waves interact with the coast.Indeed in benchmark 2 where the water oscillation involves all the basin and where waves interact all the time with the boundary, the volume error grows steadily with a rate that can be estimated to be less than 2 ‰ per period.

0
T 1 / 2  All records of nodes with negative x coordinates are in phase, while the ones corresponding positive x coordinates are phase-shifted by half period (see Fig. 8).The oscillation is dominated by the period T 1 .From records of nodes that are always wet, it is easy to also see the influence of the frequency zero that introduces a non-zero average elevation, and makes the modules of maximum and minimum values quite different from each other.Gauges 11, 12, 13 and 14 are on land at the same altitude, and are intermittently dry and wet.

Benchmark n. 4: tsunami run-up onto a complex three-dimensional beach
The M w = 7.8 Hokkaido-Nansei-Oki earthquake occurred on 12 July 1993, at 13:17 UTC, southwest of Hokkaido in the Sea of Japan.The strong quake was followed by a large tsunami that hit the nearby island of Okushiri with devastating power, killing more than 230 persons (Shuto and Matsutomi, 1995;Titov and Synolakis, 1997).Tsunami run-up heights up to 30 m were observed, with an extreme run-up height of 31.7 m that was measured at the head of a very narrow gulley near the village of Monai in Okushiri Island (Hokkaido Tsunami Survey Group, 1993).Here we have simulated the wave propagation and run-up of a laboratory experiment instead of the real event, consisting of a 1/400 scale physical model and conducted by Matsuyama and Tanaka (2001)  (one can find all the available data at http://isec.nacse.org/workshop/2004 cornell/bmark2.html).
The laboratory benchmark has been prepared by using a part of the flume and its precise dimensions are 5.446 m in length (x-axis) and 3.406 m in width (y-axis), with a still water largest depth of 0.135 m.The bathymetry and the coastal topography used in the laboratory experiment are given on a regular square-cells grid with a side length of 0.014 m.On the midpoint A of the entrance side of the tank (left side), the measured record of the water elevation of the incoming wave is known.The sidewalls as well as the topography reconstructed on the right side of the tank are fully reflective.The data available for a comparison are three aligned probes, numbered as channel 5, 7 and 9 that are located off, but close to, the coastline in the respective positions (4.521, 1.196), (4.521, 1.696), (4.521, 2.196) given in meters.Figure 16 is a plot of the basin geometry reconstructed in the laboratory tank facility where the Monai pocket beach and the on-land steep topography are reproduced as well as the little island of Muen sitting just in front of the beach.This Figure also shows the positions of all the probes.Furthermore, a movie zooming on the area inundated with the highest run-up, with frames separated by 0.05 s, is also available.The computational grid has been taken to be equal to the basin geometry mesh ( x = y = 0.014 m) and the selected time step t is 0.001 s.The laboratory tank has a frictional resistance, so frictional terms of Eqs. ( 9) and ( 10) cannot be neglected in this benchmark problem.The Manning's roughness coefficient n has been set equal to 0.025 s m −1/3 , which is adequate for the near-shore, rivers and coasts with no buildings and even for flumes like the experimental tanks (see Linsley and Franzini, 1979;Muhari et al., 2011).The boundary value problem has been solved by using the Eqs.( 23) and ( 24) for all the nodes of the left side of the grid, which means that the measured excitation functions have been assumed to be a pure incoming wave and not the total wave.All the other boundary conditions have been set as fully reflective.
The forcing function is plotted in Fig. 17 and is known until about 22 s.For subsequent times it has been assumed to be identically null.In Fig. 18 the comparison between numerical and experimental probe records is shown.Indeed, since the probes are not located exactly in correspondence of grid nodes, the nearest neighbours are selected for the graphs.The laboratory signals exhibit a quite relevant noise during the first 10 s of the experiment that is very likely to persist even later affecting the main wave recording, as also noted by Zhang and Baptista (2008).The agreement is quite good, especially in the rising side of the curves where the arrival of the steep (almost discontinuous) front is well reproduced both in time and amplitude.The descending phase is slightly delayed in the computations, with the effect of overestimating the water elevation.Consider that the water is quite shallow at the probe position (around or less than 1 cm) and that the highest water elevation is around 3-4 cm, which means that the regime is strongly non-linear.In Fig. 19 inundation snapshots zoomed on the narrow gulley where the highest run-up was observed are given and comparison is made between the frames taken from a recorded movie Matsuyama and Tanaka ( 2001) and the computed fields of water surface elevation reproduced side by side.The agreement is very good.Both the physical model and the numerical model are able to reproduce the tsunami run-up height observed at Monai.This occurs in the central panels corresponding to 16.65 s of elapsed time.Here one sees a tongue of water climbing along the gulley up to a level of about 8 cm that, taking into account the scale factor, is equivalent to 32 m.

Tsunami scenario based on the 1693 eastern Sicily event: application to the coasts of Catania, Italy
Eastern Sicily was struck by a strong seismic sequence starting on 9 January 1693.The main shock, of magnitude M w = 7.4 (CPTI Working Group, 2004), occurred on 11 January and it has been considered, most probably, to be responsible for the subsequent tsunami whose source determination is still the subject of an open debate.Tsunami observations (Baratta, 1901;Tinti et al., 2004) and tsunami modelling suggest that the source should be located offshore (Tinti et al., 2001) and a number of hypotheses have been proposed in literature on a possible offshore fault rupture (see e.g.Bianca et al., 1999;Sirovich and Pettenati, 1999;Monaco and Tortorici, 2000;Argnani et al., 2005;Gutscher et al., 2006).The most promising rupture area seems to be located along the Hyblaean-Malta escarpment since it provides the best fit with tsunami observations in terms of first wave polarity and wave amplitude (Tinti and Armigliato, 2003).Furthermore, along the Hyblaean-Malta escarpment some landslide bodies have been identified, giving support to the hypothesis of a tsunami induced by a landslide or by a combination of an earthquake and a landslide (Armigliato et al., 2007).All the studies on this tsunami quoted above made use of a tsunami propagation model based on the finite element technique (called UBO-TSUFE) that solves the NSW equations on a fixed domain, i.e. a domain with a fixed coastline consisting of a vertical wall following the contour line of a pre-selected sea depth.Therefore, technically the UBO-TSUFE model is not able to compute run-up heights, but only the maximum sea surface elevation.For the field application of the code UBO-TSUFD, we have selected this case since there is still uncertainty on the most plausible seismic fault that ruptured.In a previous study, the code UBO-TSUFD was already applied to this area (see Tonini et al., 2011b) to explore scenarios of tsunamigenic faults, both locals and remote, and to explain the technique of the worst-case credible tsunami scenario analysis (WCTSA) that was introduced in the frame of the already mentioned project SCHEMA.In the example treated here, we use a tsunamigenic source not analysed in the paper by Tonini et al. (2011b) that was however proposed in the literature as a possible source of the 1693 tsunami.The source model is the one named HM4 in the paper by Tinti et al. (2001), where the fault is represented by two disjointed segments located along the Hyblean-Malta escarpment, but the co-seismic on-fault slip has been increased from 1.7 m to 3.0 m, in order to get a magnitude M w = 7.17, closer to the latest CPTI4 estimate.The focal parameters are listed in Table 2.This source model was given great importance in Italy in studies carried out one decade ago, since it was taken as the reference earthquake to evaluate seismic damage in the town of Catania in a project involving a large part of the national community of geologists, geophysicists and earthquake engineers (see the GNDT monograph by Faccioli and Pessina (2000), and more precisely chapter 3 by Pessina).
The computational grid is shown in the left panel of metric sources of data as explained in the paper by Tonini et al. (2011b) and summarised in Table 6 of the same paper.Notice that in that paper the grids used for the tsunami simulations differ from the one used in the present application.A zoom on the area of the town of Catania and of its commercial harbour is displayed in the right panel of Fig. 20, where one can also see the position of a series of virtual tide gauges placed all in shallow waters (sea depth < 100 m) that enable us to follow the tsunami signal evolution along its way to the shore and into the harbour.
Benchmark 3 Step 5.0 m Step 2.5 m 0 0.00000 0.00005 0.00010 0.00015 0.00020 0.0000 0.0005 0.0010 0.0015 0.0020 Benchmark 2 0 1 2 3 4 5 6 t (min) -1e-7 0 1e-7 2e-7 3e-7 4e-7 5e-7 6e-7 Benchmark 1 Step = 50 Step = 10 Step = 5 0 1 2 3 4 5 6 t (min) Benchmark 1 Step = 50 Step = 10 Step = 5  The initial sea surface elevation of the tsunami, i.e. taken to be equal to the vertical displacement of the sea floor computed through Okada's formulas (Okada, 1992), is shown for the entire grid in the panel on the upper left corner of Fig. 21 that also gives the portrait of the tsunami wave fields taken after 2, 4 and 6 min in the other panels of the same upper row.The lower row displays the same fields but zoomed on a (10 km × 10 km) box around the area of Catania.The tsunami reaches the coast very soon.Within only 2 min the tsunami arrives at the Sicily coast first hitting the promontory of Monte Tauro.In about 4 min the positive leading front hits the offshore structures of the harbour of Catania and of  Augusta.Restricting the attention to the zoomed boxes of the lower row of Fig. 21, one sees that the tsunami attacks the coast with a long linear front that is almost parallel to the coastline and is mainly the effect of the northern fault segment rupture.At around 4 min, the front starts interacting with the external jetty of the Catania harbour and at about 6 min the following deep trough also reaches the jetty, while the leading wave starts flooding the La Plaia beach (i.e. the linear flat beach located south of the town in the homonymous Gulf).
Figure 22 shows the records of the sea surface elevation calculated in positions displayed in Fig. 20 (left panel).The main tsunami signature is the crest-trough system, with trough (about 2 m) much larger than the crest (about 1 m), and with around 3 min duration, i.e. visible in all records from 1-4 and in record 6 that are located offshore along the tsunami propagation path.This system is recorded at different times with increasing amplitude (trough passing from 2 to circa 3 m) as it approaches the coast, and is followed by a second system that is the front reflected from the coast.The separation of the two signals in the records is expectedly larger for nodes more distant from the coast.After these two main groups of oscillations, records tend to die out and, after one hour, sea level changes are negligible and irrelevant.In positions inside the harbour (records from 7-9), however, waveforms are quite different.The offshore tsunami signature is weaker and tends to be substituted by a series of regular oscillations that form before in the inner harbour basin (record 9) and then also in the outer basin (records 7 and 8).These are synchronous oscillations involving the harbour of Catania, very persistent, with amplitude decreasing from the innermost (record 9) to the outermost (record 7) gauge and with a period of about 15-20 min.These can be probably interpreted as Helmholtz mode oscillations of the harbour that probably remain trapped in the coastal belt, and are able to also influence the record of the coastal gauge 5, located at the beginning of the La Plaia beach, where the tsunami main signal is followed by a long-series of long-period oscillations.
Figure 23 shows maximum elevation fields, with maxima attained in the two-hour interval after the earthquake.One can distinguish clearly the contributions of the two fault segments, with the northern fault affecting chiefly the Gulf of Catania and southern fault affecting predominantly the Gulf of Augusta.Two almost parallel energy beams go ENE direction and are slightly deflected by the bathymetry.Flooding is visible in the Gulf of Catania in correspondence of the La Plaia beach (central panel) where sea penetration can reach 100 m (corresponding to two on-land cell being inundated), but it is more significant in the southern sector of the gulf around the Simeto river mouth, located in a natural park called Oasi del Simeto (right panel) where sea can penetrate up to 500 m.
The main purpose of this application in the context of this paper is to show that the code UBO-TSUFD is able to treat realistic tsunami cases that are of interest for Italy and for the Mediterranean, providing reasonable results for tsunami propagation and inundation.The code is able to capture the most interesting features of tsunami propagation, including the interaction of the tsunami with the coast and the excitation of resonance in the harbour as well as the flooding of flat beaches and especially of river mouth.The detailed comparison of the numerical results to the historical observations of the tsunami especially in the town of Catania would have required the reconstruction of the shoreline and bathymetry at the event occurrence time, i.e. in 1693.Changes are expected to be minor, with the only important exception of the harbour area.The 1693 harbour was indeed restricted only to what is today the inner west basin.The inner NS dock was the outer jetty in 1693, while all other structures were built later (see the right panel of Fig. 20).Hence the town was much more exposed to tsunami attacks than today and was severely inundated (Boccone, 1697;Tinti et al., 2004), which is not obtained in the present simulation because the computational grid accounts for the harbour geometry as it is today.

Conclusions
This paper has presented the main features of the numerical code UBO-TSUFD that solves the NSW equations by using a LF method on staggered structured grids and accounts for moving shorelines, which enables one to compute wave propagation and coastal flooding.The code solves both initial-value problems as well as boundary-condition problems, and it is therefore adequate for computing tsunamis generated by earthquakes inside the computation grid and tsunamis generated outside the grid and entering the basin Benchmark n. 1 is a 1-D case where non-linearity is quite strong, and is a modified version of one of the Catalina Island tsunami benchmarks (Liu et al., 2008), where the modification was needed to avoid triplication of the exact solution (see Fig. 4).The code appears to be capable of predicting correctly the benchmark functions and to converge to stable solutions when the grid size is made smaller and smaller.
Benchmark n. 2 is again 1-D, but in a 2-D geometry and with a coastline that is at any time curvilinear, more precisely an ellipse, which is the intersection of the planar water surface with a paraboloid representing the basin.The case treated here is mildly non-linear and a stronger non-linearity would have been obtained by taking higher values of the ratio A/L in the Eqs.( 27) and ( 28): this would have also implied the clear appearance of the higher frequency 2ω 1 (2ω 2 ) in the sea surface oscillations.However, the case is significant because run-up and run-down values are remarkably of the order of several meters and are realistic, and sea penetration and retreat involves very many cells of the grid that pass from wet to dry conditions and vice versa during the periodic motion.The code has shown to be able to reproduce correctly the sea level position and the shoreline position.
Benchmark n. 3, which consists in a positive plane wave moving over a flat channel and climbing a constant-slope beach, has been developed for successfully testing how the UBO-TSUFD model conserves the mass.
Benchmark n. 4 is an example of a wave with strong non-linearity and close to breaking that propagates basically in one direction but that, due to reflections, refractions and diffractions soon become a complex 2-D field especially after the first interaction with the coast.The incoming wave evolves quite soon in a sort of a bore with a quite steep front that goes around the islet of Muen by diffraction and climbs up the narrow gulley of Monai reaching runup heights much larger than in rest of the neighbour coast.This case has also served to test the inclusion of the bottom friction terms in the NSW equations set that are implemented via Eqs.( 9) and (10).All the available laboratory features have been reproduced correctly by the code, including the time and the modality of the inundation in Monai: time series of the probes off Monai and movie images show a very satisfactory agreement.
The last example treated is a realistic case of tsunami attacking the coast of eastern Sicily in southern Italy.The case considered is a scenario inspired by the 11 January 1693 earthquake and tsunami, which was one of the most destructive in the Italian history, and shows that the code can provide sensible results of tsunami propagation and inundation.The case is not a real reconstruction of the historical event since no attempts have been made to reconstruct the original geometry and/or geomorphology of the coast, for example, relevant is the discrepancy in the harbour area where the current geometry was used in the tsunami computations.Rather, the case can be seen as the scenario of a possible future earthquake from a source that deserves, however, more geotec-tonic validation and can be seen at the moment only as the fruit of an interesting but unproven speculation.
In conclusion, the code UBO-TSUFD can be considered as validated in its capability to compute non-linear long wave propagation also including energy dissipation through bottom friction expressed by the aid of the Manning's roughness coefficient, in 1-D and in 2-D, when the wave is excited by initial conditions (either sea level field or discharge fluxes or both) or by a boundary conditions.Excitation by a landslide and propagation in a nested grid system are left for future validation tests.

Fig. 1 .
Fig. 1.Benchmark n. 1: definition sketch, where η is the water elevation, h is the water depth, x c is the coastline position and β is the slope angle.

Fig. 2 .
Fig. 2. Sketch of the staggered grid technique.Upper panel: (left) definition of the discretised variables on a space cell and (right) on the time axis but restricted to sea surface elevation η and the flux in the x direction M. Elevation η is computed in the centre of the cell, while fluxes M and N are shifted by a half spatial step along their respective directions.Lower panel: sample of a grid with grey cells on the bottom and on the left side of the grid.These are ghost cells that are outside the computational domain where only the discharge flux component on the side in common with white cells is computed.

Fig. 3 .Fig. 4 .
Fig. 3. Benchmark n. 1: initial and later sea surface profiles computed through the code UBO-TSUFD with different grid steps compared to analytical solutions.At this scale, all curves overlap.A more detailed comparison of the shoreline area is given in Fig. 5.

Fig. 5 .
Fig. 5. Benchmark n. 1: cross sections zoomed in the shoreline area.Markers correspond to grid points at different resolutions.

Fig. 6 .Fig. 7 .
Fig. 6.Benchmark n. 1: coastline position (upper left) and velocity (upper right) predicted by UBO-TSUFD compared to the analytical solution.Residuals, i.e. differences between analytical and numerical values, are shown for coastline position (bottom left) and velocity (bottom right).Residuals have been calculated by means of a moving average over a time interval of 2.0 s.

Fig. 8 .
Fig. 8. Benchmark n. 2: contouring of the topo-bathymetry of the paraboloidal basin and positions of the calculated tide-gauge records.The black line is the coastline in still water conditions, i.e. an ellipse with semi-axes L = 4700 m and l = 1300 m.Tide gauges form couples symmetric with respect to the x-and y-axes.Examples of x symmetric couples are (5, 6), (7, 8), (9, 10), and examples of y symmetric couples are (1, 4), (2, 3) and (7, 9).Some of these nodes are in deep water and always remain wet, while others like the couples (11, 12) and (13, 14) belong to the wet-dry area.

Fig. 9 .
Fig.9.Benchmark n. 2: cross sections of the sea planar surface taken at each quarter of the period T 1 , i.e. equal to 470 s.After a full oscillation of the water surface, the agreement between analytical and numerical results is still very satisfactory, despite of the numerical noise shown in Fig.12.The values of R are the mean residuals calculated over the entire domain in m.They tend to increase along with the simulation time

Fig. 10 .
Fig.10.Benchmark n. 2: calculated tide gauge records (black dashed line) compared to the analytical solution (red solid line) for the planar oscillation along the x-axis.All records of nodes with negative x coordinates are in phase, while the ones corresponding positive x coordinates are phase-shifted by half period (see Fig.8).The oscillation is dominated by the period T 1 .From records of nodes that are always wet, it is easy to also see the influence of the frequency zero that introduces a non-zero average elevation, and makes the modules of maximum and minimum values quite different from each other.Gauges 11, 12, 13 and 14 are on land at the same altitude, and are intermittently dry and wet.

Fig. 12 .
Fig. 12. Benchmark 2: 2-D snapshots of water elevation (left column) and N flux computed by means of the code UBO-TSUFD.Differently from the theoretical problem, the velocity along the y-axis is not null and this leads to an error between numerical and analytical results which increases with time.

Fig. 13 .Fig. 14 .
Fig. 13.Benchmark 3: propagation of a positive plane wave over a flat channel ending with a constant-slope beach.

Fig. 15 .Fig. 16 .
Fig. 15.Mass conservation: time histories of the variation of the total (left column) and displaced (right column) volume/area are plotted for cases treated in benchmarks 1 and 2 as well for the specific case of benchmark 3.

Fig. 17 .
Fig. 17.Benchmark n. 4: forcing function as measured in the probe placed in position A of the tank (Fig. 16).

Fig. 19 .
Fig. 19.Benchmark n. 4: snapshot comparison.Images on the left column have been extracted by the movie recorded during the experiment(Matsuyama and Tanaka, 2001).Movie frames are separated by 0.05 s and here we have selected the numbers 45, 55 and 70.On the right column the corresponding calculated snapshot with the same time separation are given.We have assumed that frame 45 corresponds to the time 16.15 s of simulation.Space coordinates are the same as in the computational grid.Topography contour lines are at 1 cm intervals.

Table 1 .
Benchmark 1: mean ( R) and root mean square (RMS) of the residuals θ T i − θ C i for shoreline position and velocity calculated over 400 s simulation time.The residuals are plotted in Fig.6.
−1 km onshore to 50 km offshore.Lateral sidewalls which are parallel to the wave propagation are fully reflective, the right side boundary is open for free transmission and the left side is closed by the topography (100 m high).
Benchmark n. 4: comparison of time histories at the benchmark probes.Experimental records (black dashed line) are quite noisy and deflect from zero before the arrivals of the wave.The wave has a leading depression.The arrival of the positive front is simulated quite well, including the abrupt change of elevation close to the crest.The decreasing part of the record is slightly delayed in the simulations.The comparison is limited to 30 s because, beyond this time, the simulation becomes strongly non-linear (breaking waves) and NSW cannot reproduce this situation well.Moreover, the input forcing function does not consider reflecting waves which start to propagate backward and laterally before the forcing is finished.Hence this effect, leading unavoidably to misleading results, cannot be considered in the simulation.

Table 2 .
Benchmark 1: mean ( R) and root mean square (RMS) of the residuals θ T i − θ C i for shoreline position and velocity calculated over 400 s simulation time.The residuals are plotted in Fig. 4.1.