Adaptive numerical modelling of tsunami wave generation and propagation with FreeFem++

,


Introduction
Tsunami waves represent undeniably a complex natural process.Moreover, they represent a major risk for exposed coastal areas including also the local populations, infrastructure, etc.The present work is devoted to the modelling tsunami generation and propagation processes.Moreover, this article is designed as a tutorial paper in order to show to the readers how easily these processes can be modelled in the framework of the FreeFem++ open source finite element software.Traditionally, tsunami waves are modelled using hydrostatic models [1][2][3][4].In the present manuscript we employ a non-hydrostatic BOUSSINESQ-type system to be specified below.This class of models is distinguished by the application of the so-called BOUSSINESQ approximation [5].They can be used to study a variety of water wave phenomena in harbors, coastal dynamics and, of course, tsunami generation and propagation problems [6][7][8][9][10].
In this study we consider a BBM-BBM system derived by D. MITSOTAKIS in 2D over a variable bottom in space h(x, y) and in time ζ(x, y, t) [11]: Submitted to Geosciences, pages 1 -24 www.mdpi.com/journal/geoscienceswhere and g is the acceleration due to gravity.System ( 1) is an asymptotic approximation to the three-dimensional full EULER equations describing the irrotational free surface flow of an ideal fluid Ω ⊂ R 3 [12,13], which is bounded below by −z b (x, y, t) = −h(x, y) − ζ(x, y, t) and above by the free surface elevation η(x, y, t) (cf. Figure 1).
The variables in (1) are X = (x, y) ∈ Ω and t > 0 are proportional to position along the channel and time, respectively.η = η(X, t) being proportional to the deviation of the free surface departing from its rest position and V = V(X, t) = u(X, t) v(X, t) = (u, v) = (u; v) being proportional to the horizontal velocity of the fluid at some height.In our study, we suppose that η = O(a), with the characteristic wave amplitude a (in other words, η is the difference between the water free surface and the still water level).Also we set λ = O( ) be the wave length.In addition, we limit ourselves to the case where η + z b > 0 (there are no dry zones in our computations).
This paper is organized as follows.In Section 2 we present the space and time discretization of Equations (1).In Section 3, we present the new domain adaptation technique.In Section 4, we establish the convergence of our numerical code, which validates the adequacy of the chosen finite element discretization.Then, with this code we simulate the propagation of a tsunami-like wave generated by the moving bottom (e.g. an earthquake).We present several test cases in various regions of the world.
First, we take a MEDITERRANEAN sea-shaped computational domain with flat bottom and we solve the sBBM system (1) in it.The mesh in this study is generated from a space image.Then, we consider the JAVA island region with real world bathymetry.Finally, we apply this solver to simulate a realistic example of a tsunami wave near the JAVA island which took place in 2006.We note that all numerical simulations were done using the FreeFem++ software [14], which is an open source platform to solve Partial Differential Equations (PDEs) numerically, based on Finite Element Methods (FEM).The main conclusions of this study are outlined in Section 5.

Discretization of the BBM-BBM system
In this section, we present the spatial discretization of (1) using Finite Element Method (FEM) with P 1 continuous piecewise linear elements.For the time marching scheme we use an explicit second order RUNGE-KUTTA method.

Spatial discretization
We let Ω be a convex, plane domain, and T h be a regular, quasi-uniform triangulation of Ω with triangles of maximum size h < 1. Setting V h = {v h ∈ C 0 ( Ω); v h | T ∈ P 1 (T), ∀T ∈ T h } be a finite-dimensional, where P 1 is the set of all polynomials of degree ≤ 1 with real coefficients and denoting by •; • the standard L 2 inner product on Ω, we consider the weak formulation of System (1): , we have: (2) For simplicity, we set , so that system (2) can be rewritten in the following way: (3) After integrating by parts, the left hand side of (3) becomes: and where Γ n is the boundary of the domain Ω.Dealing with the right-hand side F (E , U , V, Φ η ) of the first equation in System (3), we expand the two complex terms which are multiplied by A and Ã such as: On the other hand, we have: and, consequently, we deduce the final form of F (E , U , V, Φ η ) as follows: For the right-hand side G (E , U , V, Φ u ) of the second equation in System (3), we have: Finally, for the right-hand side H (E , U , V, Φ v ) of the third equation in System 3, we have: However, the model presented above contains some drawbacks.In particular, when the bathymetry function contains steep gradients, it causes instabilities in the numerical solution.We have to mention that this problem is well-known in the framework of BOUSSINESQ-type equations [15].In order to avoid this kind of problems and to have a robust numerical model, we take two measures.First of all, we perform the smoothing of the bathymetry data which is fed into the model.In this way, we avoid noise in the bathymetry gradient.As a second and more radical step, we neglect higher order derivatives of the bathymetry function as it was proposed earlier in [11].Thus, from now on we shall use the following system of equations: and

Time marching scheme
Our method is based on the explicit second order RUNGE-KUTTA scheme.For that, let us denote by (E n+1 , U n+1 , V n+1 ) and (E n , U n , V n ) the approximate values at time t = t n+1 and t = t n , respectively and by δt the time step size.Then, owing to (4), the unknown fields at time t = t n+1 are defined as the solution of the following system: where ) and

New domain adaptation, domains computation and initial data
We present here the new domain adaptation technique that will be compared in the sequel with the mesh adaptation used in FreeFem++.

New domain adaptation technique
Since some computation domains for many applications (here for Tsunami waves) may be huge and the initial data is concentrated in a small domain, before starting to propagate in the domain, we present here an idea to build a moving computation domain around the solution only, as when we use a mesh adaptation.The difference between these two methods is that the moving domain will be a cut from the initial one; i.e. all initials vertices, edges and boundary labels are conserved and a new label is defined for the new boundary; while the mesh adaptation technique don't conserve the initials vertices and edges, so when we make interpolation of solution from old to new mesh we will lose some information in the mesh adaptation technique but not with the moving domain.
Firstly, we cut from the initial mesh Thinit a circle or a rectangle zone Th where our initial solution lives (using trunc in FreeFem++), secondly we compute the initial solution u0 and we interpolate it to uadapt ∈ P 1 finite element (using interpolate in FreeFem++), and for each adaptation we follow this algorithm: • We deduce the limit min max of Th on x and y direction (using boundingbox in FreeFem++).
• We add epsadapt from each side in order to build the new rectangle Th1 (cutted from Thinit) that contains Th (using trunc in FreeFem++).• We interpolate uadapt ∈ Th to uadapt1 ∈ Th1 (using interpolate in FreeFem++).
• We smooth the function obtained from abs(uadapt1) ≥ erradapt using: where f = (|uadapt1| ≥ erradapt), with zero DIRICHLET BC only on the new boundary label of Th1 and a NEUMANN BC in the other boundary label.
Finally, we replace Th by the new mesh Thnew and we interpolate the solution from the old mesh to the new one (using interpolate in FreeFem++).We use a reflective Boundary Condition (BC) on the new boundary, i.e. zero NEUMANN BC for η and zero DIRICHLET BC for V, cause our BBM-BBM system gives artificial numerical explosion on the boundary if we do not use any BC or if we use only NEUMANN BC for η and V.
For the BBM-BBM system over a flat bottom, we use a mesh generated through a photo of the MEDITERRANEAN sea (a cut of the mesh around the CRETE island is shown in Figure 2 at left panel) and for the BBM-BBM system over a variable bottom in space and in time, we use a mesh generated using an imported bathymetry f xy for the sea near the JAVA island which can be downloaded from the NOAA1 website where in this case, the mesh generated is for the area, where the amplitude is zero.We can smooth the bathymetric data obtained from NOAA (cf. Figure 3, left panel) by solving (8) with f = f xy .For all simulations with realistic bathymetry, we use β = 20 in (8) to smooth the initial bathymetry after the generation of the mesh (cf. Figure 3, right panel) in order to ensure the stability of the numerical method, we also note that in order to be in a big deep water wave regime for the BBM-BBM system we change the depth close to the shoreline to 100 m.
The bathymetry data downloaded from the NOAA website are in degree coordinate and we need to convert them to meters.So, on the first hand, we must know the degree of Latitude (South and North) and of Longitude (West and East) of our domain where we can deduce the Latitude lat0 = .5(latSouth + lat North ) and the Longitude long0 = .5(longWest + long East ).On the other hand, we must take into account the spherical shape of the EARTH, even if it does not play significant role because of the small spatial scale of the experiments.So, we know that the radius of the EARTH near the equator is R equator = 6378, 137 km, and near to the pole R pole = 6356, 752 km, thus the radius of our domain equals to: So, we move the mesh of our domain using the following translation (coefl0 = πR/180):

Initial data 92
Tsunami waves considered in this study are generated by the co-seismic deformation of the Ocean's or sea's bottom due to an earthquake.The adopted modelling of the tsunami wave generation process is inspired by [8,11,16,17].The co-seismic displacement is computed according to the celebrated OKADA's solution [18,19].We assume the dip-slip dislocation process underlying the earthquake.The vertical component of displacement vector O(x, y) is given by the following formulas employing CHINNERY's notation, cf.[16,17]: and Here, W and L are the width and the length of the rectangular fault, (x, y) are the points where we computes displacements, (x 0 , y 0 ) is the epicenter, d = fault depth(x 0 , y 0 ) + W sin δ, δ is the dip angle, θ is the rake angle, D is the BURGERS's vector, U = |D| sin θ is the slip on the fault, φ is the strike angle which is measured conventionally in the counter-clockwise direction from the North (cf. Figure 4 (left)), µ, λ are the LAMÉ constants derived from elastic-wave velocities: where ρ c is the crust density, V P is the compressional-wave (P−wave) velocity, V S is the shear-wave (S−wave) velocity.The Matlab script to compute the OKADA solution can be downloaded at the following URL: https://mathworks.com/matlabcentral/fileexchange/39819-okada-solution/ We shall distinguish here the two types of tsunami wave generation mechanisms [20,21]: active and passive generation mechanisms.

Passive generation
We remind that the passive generation approach consists in transposing the bottom deformation on the free surface as an initial condition for tsunami propagation codes.In order to compute the initial data for η(x, y, 0) = O(x, y) in meters (cf. Figure 4 (right)), V(x, y, 0) = 0 which is referred to as a passive generation of a tsunami wave near the JAVA island, using our domain adaptive technique, we will use the fact that the solution is concentrated in the small rectangle [x 0 − 3.2W; where L = 100 km, W = 50 km, δ = 10.35In contrast to passive generation, the active generation approach consists in generating a tsunami waves by computing fluid layer interaction with moving bottom.For a more realistic case of the JAVA 2006 event, we use precisely this so-called active generation approach by following [8,22].In this case we consider zero initial conditions for both the free surface elevation and the velocity field, and assume that the bottom is moving in time.This case may be described by considering the bottom motion formula: where N x sub-faults along strike and N y sub-faults down the dip angle, H(t) is the Heaviside step function and α = log(3)/t r , where t r = 8 s is the rise time.We choose here an exponential scenario, but in practice, various scenarios could be used (instantaneous, linear, trigonometric, etc.) and could be found in [8,16,17,22,23].Parameters such as sub-fault location (x i , y i ), depth d i , slip U and rake angle θ for each segment are given in [8, Table 3].In this table, we notice that the fault's surface is conventionally divided into N x = 21 sub-faults along strike and N y = 7 sub-faults down the dip angle, leading to a total number of N x × N y = 147 equal segments.
For our special domain adaptivity technique, since the fault plane is considered to be the We show in Figure 6, the bottom displacement ζ(x, y, t) at time t = 100 s and t = 270 s using our domain adaptation technique.We note that after building the OKADA solution O(x, y) in the passive generation or O i (x, y) in the active generation, we can remark that this solution is non-local and decays slowly to zero, that is why in our domain adaptation technique we put 0 where the absolute value of the solution is less then min(|min We make the same thing without adaptive mesh in order to compare the solution using the same initial data.

Numerical simulations
In this section, we study first the rate of convergence of our schemes for the BBM-BBM System (4) with non-dimensional and unscaled variables i.e., with g = 1 over a variable bottom in space, which establishes the adequacy of the chosen finite element discretization and the used time marching scheme, for the flat bottom case, we refer to [24], where we use the same technique as in this paper.Then, we simulate the propagation of a wave, that is similar to a real-world tsunami wave generated by an earthquake, in the MEDITERRANEAN sea with the BBM-BBM model over a flat bottom.Then, we switch to the JAVA island region with real variable bottom in space.Finally, we study the active tsunami generation scenario which took place in 2006 near the JAVA island.In all numerical simulations we used P 1 continuous piecewise linear functions for η, u, v, h and ζ.

Rate of convergence
We present the evidence here, following the work done for the 1D case of the BBM-BBM system in [25], that the second order RUNGE-KUTTA time scheme considered for the BBM-BBM variable bottom in space is of order 2. We note that the function ζ(x, y, t) is only used for the generation of tsunami wave and, thus, will not be taken into account in the convergence rate test.In this example, we take bi-periodic Boundary Conditions (BC) for η h , u h and v h on the whole boundary of the square [0, 2L] × [0, 2L], where L = 50 and we consider the following exact solutions: adding an appropriate function to the right-hand side to make these solutions exact.We measure at time T = 1 and for and we end up with the results reported in Table 1.So, the L 2 rates for η and V is of order ∼ 2 and the H 1 rates for η and V is of order ∼ 1 as shown in the Figure 7 and which confirms the convergence of the second-order RUNGE-KUTTA scheme in time for the BBM-BBM system with variable bottom in space.
η, slope = 1.995V, slope = 1.992Order of accuracy in loglog scale for the H 1 norm of η and V.

Figure 7.
Rate of convergence for sBBM system with variable bottom in space.

Propagation of a tsunami wave in the Mediterranean sea with a flat bottom
We simulate here, the propagation of a wave that looks like a tsunami wave generated by an earthquake in the MEDITERRANEAN sea with the BBM-BBM System (4) with a flat bottom −h(x, y) = −1, 5 km which is the average depth of the MEDITERRANEAN sea.This wave was defined above in the passive generation part of the Section 3 where, in this case, the initial solution is concentrated in the small rectangle [x 0 − 5W; x 0 + 4W] × [y 0 − 1.5L; y 0 + 2.5L] and we take these following values: is the POISSON's ratio, U = 2.5 m, (x 0 ; y 0 ) = (2390.* scale, 590.* scale) and the fault depth 10 km.
In this example, we will take the fact that the LAMÉ constants µ and λ are given by the formulas We also use the following settings: for the step time δt = 1 s, a reflective BC for all the boundary, for the adaptmesh of FreeFem++: f e s p a c e V h i n i t ( T h i n i t , P0 ) ; V h i n i t hT= h T r i a n g l e ; r e a l Dx=hT [ ] .min ; Th=adaptmesh ( Th , uadapt , e r r = 1 .e −7, e r r g = 1 .e −2,hmin=Dx , i s o =tr ue , nbvx=1 e8 ) ; and for our domain adapt technique: isoadapt=5e-2, erradapt=1e-7, β=5e-3, epsadapt=2e-2.We note that, we adapt the mesh around the solution each 100 iterations i.e. each 10 s by using uadapt= In order to compare the results between adaptmesh of FreeFem++, our new domain adaptation technique and without using mesh adaptation, we plot in addition to the free surface elevation η in the Figures 8 → 9, the variation of η vs time in Figure 10   without mesh adaptation which is very promising method for the tsunami wave propagation.For the adaptmesh of FreeFem++ with err=1.e-7and errg=1.e-2,we also almost get the mass conservation with an error of order 9.5e − 4, but we obtain some difference in wave gauge with the full method which is due to the refinement mesh adaptation and the interpolation of the solution, although the computation time is almost the double of the new domain adaptation technique.Thus, we can go faster with our new domain adaptation technique if we can also deduce the mass matrix after cutting the mesh, of course, if the mass matrix does not change along the simulation of the full mesh.This is an outgoing project.

Propagation of a tsunami wave near the Java island: passive generation
We will take here the same initial data as defined above in the passive generation part of Section 3, we take δt = 1 s as the time step size and we note that, we adapt the mesh after computing the initial data for η and then every 50 s by using the following value for the domain adaptation uadapt= η + u + v, isoadapt=3e-2, erradapt=1e-4, β=5e-9, epsadapt=30e3.We compare here the results between our new domain adaptation technique and without using mesh adaptation.To this end, we plot the free surface elevation η in the Figures 13 and 14, the variation of η vs time (in Figure 15) at four numerical wave gauges placed at the following locations: (i) (107.345 (105.9 • , −10.35 • ) and (iv) (107.7 • , −11 • ) (see Figure 5 (left)) where (i) is the position of the epicenter.
However, because of the large variations of the bottom, shorter waves were generated, especially around CHRISTMAS Island (southwest of JAVA) and around the undersea canyon near the earthquake epicenter.
Finally, we present a comparison of the kinetic, potential and total energies with the full mesh (in Figure 16, top left panel) and with the domain adaptivity method (in Figure 16, top right panel) defined in [26] as follows: |V| 2 dz dxdy,  where ρ w = 1027 kg/m 3 is the ocean water density, the number of degrees of freedom (in Figure 16, down left panel) and the computation time of the simulation (in Figure 16, down right panel).We obtain here an error of order 2.6e − 4 between the total energy with domain adaptivity and without any adaptation.We present in Figure 17 the comparison of the maximum of the propagation of the solution between the full and domain adaptivity methods at t = 1750 s.

Propagation of a tsunami wave near the Java island: active generation
For a more realistic case as in the JAVA 2006 event, we use the active generation in order to model the generation of a tsunami wave as in [8,22].In this case we consider zero initial conditions for both the surface elevation and the velocity field, we take δt = 2 s as the time step size, we assume that the bottom described in the Section 3 is moving in time and we note that we adapt the mesh, before the end of the generation time t = 270 s, every three iterations i.e. every 6 s by using the following value for the domain adaptation uadapt = η + u + v, isoadapt=5e-2, erradapt=1e-4, β=5e-9, epsadapt=50e3 and then for t > 270 s every 25 iterations i.e. each 50 s.We compare here only the results between our new domain adaptation technique and without using mesh adaptation.To this end, we plot the free surface elevation η in the Figures 18 → 20.However, as in the passive case, because of the large variations of the bottom, shorter waves were generated, especially around the CHRISTMAS Island     (southwest of JAVA island) and around the undersea canyon near the earthquake epicenter.We plot the variation of η vs time (in Figure 21) at four numerical wave gauges placed at the following locations: (i) (107.345• , −9.295 • ), (ii) (106.5 • , −8 • ), (iii) (105.9 • , −10.35 • ) and (iv) (107.7 • , −11 • ) (see Figure 5 (left panel)) where (i) is the position of the epicenter.Finally, we present a comparison of the kinetic, potential and total energies with the full mesh (in Figure 22, top left panel) and with the domain adaptivity method (in Figure 22, top right) defined in (9), the number of the degrees of freedom (in Figure 22, lower left panel) and the computation time of the simulation (in Figure 22, lower right panel).We obtain here an error of order 2e − 5 between the total energy with domain adaptivity and without any adaptation.We present in Figure 23 the comparison of the maximum of the propagation of the solution between the full and domain adaptivity method at t = 1750 s.

Conclusion and Outlook
In this manuscript we demonstrated how to discretize the BBM-BBM system (1) using the FEM and dedicated open-source software FreeFem++.The use of this numerical technique was demonstrated in view of applications to tsunami wave modelling [17,27].The concrete cases of wave propagation in the MEDITERRANEAN sea and in JAVA island region (INDONESIA) were considered.The digital computing environment that we developed allows the integration of realistic data (bathymetry and geography) in a relatively simple software framework.The codes used in this study are made freely available for all our readers.Moreover, a novel mesh and domain adaptation technique was proposed to speed-up substantially the computations.The gain in terms of the CPU time after applying this technique can be clearly seen in Figure 22.The accuracy of the 'accelerated' solution is more than acceptable to make this technique useful in a variety of tsunami propagation problems.It goes without saying that this technique can be applied to other events and other regions of the world with minimal changes in the provided codes.
Regarding the perspectives of this study, we would like to develop also the parallel version of this code together with the domain adaptation technique to make computations practically faster than the real time tsunami wave propagation.However, we underline that even the current version can be efficiently run even on a modest laptop personal computer.

Figure 1 .
Figure 1.The sketch of the physical domain Ω .

Figure 2 .
Figure 2. Left (a): the mesh around CRETE island.Right (b): the place of : wave gauge and : epicenter.

Figure 4 .
Figure 4. Geometry of the source model (left) and the initial solution for η (right, min= −0.46 m, max= 0.71 m).
rectangle with vertices located at (109.20508 • (Lon),−10.37387• (Lat)), (106.50434• (Lon),−9.45925• (Lat)), (106.72382• (Lon),−8.82807• (Lat)) and (109.42455• (Lon),−9.74269• (Lat)), we will consider that our bottom displacement is concentrated on the big rectangle which is equidistant of 1 • from each side of the initial fault plane as in Figure5(left panel), then we compute each OKADA solution O i on a circle of center (x i − 10m, y i − 10m) and of radius 6 max(L, W) and at the end all the OKADA solution will be interpolated on the big rectangle before starting to compute the vertical displacement of the bottom ζ(x, y, t), in Figure5(right panel) we plot O 14 .For the computation of ζ(x, y, t), we start the mesh by a circle of center (x c − 5m, y c − 5m) and of radius 4 max(L, W) and we adapt the mesh each 3 iterations i.e. each 6 s by using the following value for the domain adaptation uadapt = ζ, isoadapt=5e-2, erradapt=1e-4, β=5e-9, epsadapt=50e3.

Figure 8 .
Figure 8.The mesh and the solution at t = 1000 s, with the full method at the left panel (a), the domain adaptive method at the center panel (b) and the FreeFem++ adaptation with err=1.e-7 at the right panel (c).

Figure 9 .Figure 10 .
Figure 9.The solution at t = 3000 s, with the full method at the left panel (a), the domain adaptive method at the center panel (b) and the FreeFem++ adaptation with err=1.e-7 at the right panel (c) (min= −5.5e −2 m, max= 4.7e −2 m, for the three case).

Figure 11 .
Figure 11.Comparison between the three method full (top panel), domain adaptation (middle panel) and FreeFem++ adaptation (down panel) of the maximum of the propagation of the solution of a tsunami wave in the MEDITERRANEAN sea for t = 6800 s (min= 0 m, max= .4m, for three cases).

Figure 12 .Figure 13 .
Figure 12.Comparison between the three methods: full, FreeFem++ adaptation and domain adaptation of the computation time of each adaptmesh, the number of degree of freedom and the computation time of the simulation.

Figure 16 .
Figure 16.Passive generation: comparison between the two methods the full one and domain adaptivity of the kinetic, potential and total energies, the number of degree of freedom and the computation time of the simulation.

Figure 17 .
Figure 17.Passive generation: comparison between the maximum of the solution at t = 1750 s, with the domain adaptivity method (left panel) and with the full one (right panel).

Figure 21 .Figure 22 .
Figure 21.Active generation: comparison between the two methods (the full one and domain adaptivity) of the free surface elevations (in meters) vs time (in seconds), computed numerically at four wave gauges where the gauge (i) correspond to the epicenter.

Figure 23 .
Figure 23.Active generation: comparison between the maximum of the solution at t = 1750 s, with the domain adaptivity method (left panel) and with the full one (right panel).
x 0 ; y 0 ) = (107.345• , −9.295 • ) and the fault depth 10 km.All these geophysical parameters can be downloaded from this file hosted by USGS: