Adaptive modelling of long-distance wave propagation and fine-scale flooding during the Tohoku tsunami

,


Introduction
On 11th March 2011 an unexpected very large earthquake struck the northeastern shore of Honshu: the Tohoku region.While the earthquake itself caused remarkably little damage -a testament to the high quality of earthquake preparedness in Japan -the subsequent tsunami caused the loss of more than 15,000 lives and extensively damaged infrastructure, overwhelming the extended network of sea defences.This disaster dramatically underlined the vulnerability of coastal communities to hazards such as tsunamis, even in countries as well prepared as Japan.A positive note is that thanks to the extensive, long-term investment in monitoring infrastructure in Japan, this event is certainly the best ever documented large earthquake and tsunami in history.Extensive datasets exist documenting ground deformation, offshore and inshore tsunami waves, global atmospheric waves, extent of inundation and damage etc...This wealth of data will hopefully lead to improved understanding of these extreme events, and of our capability to mitigate their consequences on society.One of the scientific challenges is to be able to construct a consistent model of the sequence of events leading to the observed timeseries for very different but tightly coupled processes.In this context, the main goal of this article is to demonstrate how adaptive tsunami modelling can effectively be used to obtain a consistent model of both long-distance tsunami wave propagation and local flooding.
In adaptive models the spatial resolution is not fixed in time but varies according to the local complexity of the solution.This is different from both constant-resolution models (which use a resolution constant in space and time) and variable-resolution models (which use a resolution variable in space but constant in time).The increased flexibility of adaptive methods is particularly valuable for processes which involve a wide and time-varying range of scales.Many fundamental physical processes, such as turbulence, belong to this category and tsunami waves are a clear geophysical example of such scalings.As demonstrated in [Popinet (2011c)] exploiting the "gaps" in the scale distribution of typical tsunami waves leads to a computational cost (i.e.number of grid points) which scales roughly like ∆ −1.4 , with ∆ the mesh size, rather than ∆ −2 for constant-or variable-resolution models.This reduction in the effective dimension (or "fractal dimension") of the problem leads to increasingly large gains in computational efficiency with spatial resolution.
While adaptive methods are now relatively common for industrial computational fluid dynamics applications [Yiu et al. (1996), Popinet (2003), Popinet (2009)], they are comparatively less developed for geophysical fluid dynamics [Liang et al. (2004), Popinet and Rickard (2007), Liang and Borthwick (2009), Popinet et al. (2010)].The quadtreeadaptive tsunami model used in this article, implemented within the Gerris flow solver framework is one of the few adaptive tsunami models published [George and LeVeque (2008), Harig et al. (2008), Popinet (2011c)].A summary of the model is given in the next section.Section 3 describes a classical validation test case and section 4 gives a detailed comparison between field data and model results for the Tohoku tsunami.

Numerical method
This section gives an overall summary of the numerical method.A detailed description can be found in [Popinet (2011c)] and references therein.
Tsunamis are most commonly modelled using a long-wave approximation of the mass and momentum conservation equations for a fluid with a free-surface.In this approximation, the slopes of both the free-surface and the bathymetry are assumed to be vanishingly small.Although these assumptions are often violated in the case of tsunamis, particularly close to shore, the long-wave approximation still gives reasonable results in practice.Better approximations include the Boussinesq wave equations which tend to be significantly more expensive to solve [Grilli et al. (2007)].
The long-wave approximation can be obtained either by formal vertical averaging of the Navier-Stokes equations with a free-surface or by considering mass and momentum conservation through vertical slices in the water column, as first done by Saint-Venant [de Saint-Venant (1871)].The resulting equations of motion are often written in differential form as the set of non-linear advection equations with h the free surface elevation, u = (u, v) the horizontal velocity vector, g the acceleration of gravity and z the depth of the bathymetry.This corresponds to a hyperbolic system of conservation laws which, using Stokes' theorem, can be expressed in integral form as where Ω is a given subset of space, ∂Ω its boundary and n the unit normal vector on this boundary.For conservation of mass and momentum in the shallow-water context, Ω is a subset of bidimensional space and q and f are written Dissipative terms (viscosity and/or friction) are also often added but they do not change the fundamental structure of the equations.This system of conservation laws is hyperbolic and admits an entropy inequality, which guarantees that the total energy h u 2 /2 + g h 2 /2 behaves physically.Another important property is that the steady state is a solution of this system.This simply reflects the expected property that a lake at rest remains so.Numerical methods for hyperbolic systems of conservation laws have been extensively developed in the context of compressible gas dynamics with applications in aeronautics in particular.
Many of these schemes are suitable for solving the Saint-Venant equations provided two important specificities are taken into account: the solver needs to properly account for vanishing fluid depths (equivalent to a vanishing density i.e. a rarefaction wave) and needs to verify exactly the lake-at-rest equilibrium solution (i.e.be a "well-balanced" scheme).The Saint-Venant solver used in this article is based on the numerical scheme analysed in detail by [Audusse et al. (2004)] which verifies both properties.The scheme uses a second-order, slopelimited Godunov discretisation.The Riemann problem necessary to obtain Godunov fluxes is solved using an approximate HLLC (Harten-Lax-van Leer Contact) solver which guarantees positivity of the water depth.The choice of slope limiter is important as it controls the amount of numerical dissipation required to stabilise the solution close to discontinuities (i.e.hydraulic jumps and contact discontinuities).Based on the study in [Popinet (2011c)], the Sweby limiter seems to be a good compromise between stability of short waves and low dissipation of long waves.A MUSCL-type discretisation is used to generalise the method to second-order timestepping.The timestep is constrained by a classical CFL condition with ∆ the mesh size.To summarise, the overall scheme is second-order accurate in space and time, preserves the positivity of the water depth and the lake-at-rest condition and is volumeand momentum-conserving.This scheme is implemented within the Gerris Flow Solver framework [Popinet (2003), Popinet (2003Popinet ( -2011))] which provides a number of additional capabilities including: quadtreebased adaptive discretisation, generalised orthogonal coordinates [Popinet (2011c)] (used for longitude-latitude discretisation in this article) and parallelism [Agbaglah et al. (2011)].An example of the quadtree structure is represented in Figure 1 together with its logical (tree) representation.This tree can be conveniently described as a "family tree" where each parent cell can have zero or four children cells.An important parameter is the level of a given cell in the tree.The root cell has level zero by convention and the level increases by one for each successive generation in the tree.The level of the cells in the tree are also given.
Quadtree discretisations are a good compromise between the simplicity of (inflexible) regular Cartesian meshes and the flexibility of (complex) unstructured meshes.New cells are easily added or removed from the hierarchy and fine-grained control of the spatial resolution is possible (in contrast to block-structured AMR for example [George and LeVeque (2008)]).Interpolation on newly-refined or coarsened cells is also designed to ensure conservation of quantities.
In the context of tsunami modelling, using adaptive meshes also requires a technique for efficient reconstruction of the detailed bathymetry.With spatial resolutions varying over several orders of magnitude (e.g. from 250 kilometres down to 250 metres for the results in this article) it is generally not possible or efficient to reconstruct and store the bathymetry at the finest resolution in main memory.A terrain database system has been developed specifically for this task within Gerris.It also relies on a tree-like hierarchy for fast access to very large bathymetric datasets.In constrast to the quadtree structure stored in main memory, the bathymetry database relies on a 2d-tree data structure [Bentley and Friedman (1979)] mostly stored on-disk so that the size of the database is disk-and not memory-limited.With this data structure the cost of constructing the average depth of a discretisation cell scales like O( N √ ) + O( m √ ) where N is the total number of samples in the database and m is the number of samples contained in the cell [Popinet (2011c)].In practice the typical cost of adaptive bathymetry reconstruction is roughly 10% of the total cost.

Tsunami runup onto a plane beach
The implementation of the Saint-Venant solver within Gerris has been validated both with simple test cases and complex tsunami models [Popinet (2011c), Popinet (2010a), Popinet (2010b), Popinet (2010c)].For the sake of completeness, this section describes a further test case which has not been previously published.The test case was proposed at the Third International Workshop on Long-Wave Runup Models [Liu (2004)] and seeks to validate the Saint-Venant solution of a plane wave running up a sloping beach.The geometry is chosen to replicate the spatial scales of the vertical cross-section of a typical (large) tsunami wave.The problem is thus one-dimensional for the (depth-averaged) Saint-Venant equations.The domain is 60 kilometres long, the slope of the beach is 1/10.The initial wave profile is given numerically [IWLWRM ( 2004)] and has maximum and minimum amplitudes of + 3 and − 8.8 metres respectively, located approximately 20 km and 8 km offshore respectively.The evolution with time of the wave profile close to the shore is illustrated in Figure 2. The (semi)-analytical solution is obtained using the initial-value-problem (IVP) technique introduced by Carrier et al [Carrier et al. (2003)] and is given numerically on the web site [IWLWRM (2004)].The initial wave minimum is located closer to the shore than the maximum so that the first wave is negative with a maximum dryout of the beach at around t = 175 seconds, the positive wave then runs up the beach and reaches a maximum elevation of approximately 15 metres at t = 220 seconds.For this test case only static mesh refinement is used.The spatial resolution is increased closer inshore according to where l is the level of refinement, x is the distance to the shoreline, L = 50 km, l min = 7 and l max is varied between 11 and 14 to study the convergence of the solution with spatial resolution.It is clear from Figure 2 that the numerical solution at the maximum resolution is a good approximation of the analytical solution.A quantitative assessment of the quality of the numerical solution is given in Figure 3.The average (L1-norm) and maximum errors at t = 220 seconds are computed as where (x i , y i ) is the discrete numerical solution for the wave profile, y(x i ) is the (interpolated) analytical solution at position x i and ∆ i is the (variable) grid size.The average norm shows first-order convergence and the maximum error slightly less than first-order convergence.This is not unexpected given that errors are dominated by the contact discontinuities at the wet-dry transition, where the Godunov scheme (or any other shock-capturing scheme) reduces to firstorder accuracy.From a tsunami modelling perspective, this test case demonstrates that, even for very simple coastline geometries, spatial resolutions of less than a few tens of metres are necessary to accurately capture wave runup and inundation on the shoreline.All the data necessary to reproduce this test case are given on the Gerris web site [Popinet (2011a)].

Results for the Tohoku tsunami
The adaptive model was used to simulate the 11th March 2011 Tohoku tsunami.The simulation domain was chosen so as to validate both long distance tsunami propagation and high-resolution local inundation.The Saint-Venant equations are solved on a longitude-latitude grid centered on the epicentre of the earthquake (142.597• E, 38.486 • N) and spanning 73 • in both latitude and longitude (Figure 4).The topography is obtained from two datasets: the GEBCO_08 30 arcsecond dataset [IOC (2011a)] and the SRTM 3 arc-second dataset [NASA ( 2011)] (for the Japanese islands only).Both datasets are stored using the 2d-tree (KDT) hierarchical terrain database system of Gerris with on-disk storage sizes of: • GEBCO_08: global dataset, 933 million points, 21 GBytes, KDT index: 65 MBytes • SRTM, 3 arc-second: Japan only, 58 million points, 1.4 GBytes, KDT index: 4.1 MBytes These on-disk datasets are accessed during the simulation to reconstruct the bathymetry as required as the spatial resolution evolves.The size of the KDT index gives an order of magnitude estimate of the amount of dynamic memory (RAM) required to access each dataset.A critical part of any tsunami model is the choice of an appropriate initial surface elevation (the "source model").In this study, I have chosen not to use a source model based on inversion of DART buoys and tidal gauges timeseries but rather based directly on seismic inversions.The rationale for this choice is that this allows independent validation of the source model using the DART buoys and tide gauges data, and also that, in a realtime forecasting scenario, using only seismic information may allow a quicker estimate of the source (because seismic waves travel faster than gravity waves).
Several deformation models for the Tohoku earthquake have been made available by the seismic research community, and some of them were released only a few hours after the event.A selection of such models is available on the USGS Tohoku web site [ USGS (2011)].Different methods can be used for the seismic inversion (combined with different weighting for the various seismic datasets) and the scatter in the resulting predicted deformations can be large.For this study I have used one of the deformation models obtained by Shao, Li and Ji, UCSB, more precisely their preferred "Preliminary model III" released on March 14th 2011 [Shao et al. (2011)].This model predicts vertical deformations which are significantly larger than any of the other models mentioned on the USGS Tohoku web site and was the only source model I tested which resulted in satisfying agreement with tsunami wave measurements.The initial vertical wave ele-vation was assumed to be equal to the total vertical bottom displacement predicted by the source model.This displacement was obtained by superposition of the 190 Okada subfaults displacements used in the seismic inversion.The fault model also includes the relative rupturing times of each subfault, but I chose to assume that the total rupturing time (3 minutes) was small enough to be negligible.The resulting initial vertical displacement is illustrated in Figure 5.For the adaptive model, the choice of spatial resolution is very flexible and can be based on several simultaneous criteria, depending on which specific information one would like to extract from the model.In this study I chose to illustrate how adaptivity can help obtain simultaneous solutions for two aspects of tsunami modelling which are generally considered as conflicting in term of resolution requirements: long-distance propagation and detailed inundation/damage predictions.In this study, both regimes are captured using a single criterion but different maximum spatial resolutions: the mesh size ∆ is adapted so that where ǫ is a control parameter set to 2.5 cm.ǫ can be interpreted as the maximum truncation error in the free surface elevation of a spatially first-order discretisation scheme.This corresponds to the maximum error expected near sharp discontinuities (i.e.shocks).In smooth regions, the scheme is second-order and errors will be significantly smaller.There are primarily two regions where one can expect larger discretisation errors: wave fronts and "breaking" waves in shallow water i.e. the two regimes of long-distance wave propagation and inundation we are interested in.I have shown in a previous study that adequate propagation of long-wave energy over large distances with the non-linear, shock-capturing numerical scheme requires spatial resolutions of the order of one arc-minute in order to minimise numerical wave dissipation [Popinet (2011c)].Following this, the maximum resolution for deep water wave propagation is set to 12 levels of quadtree refinement, giving a maximum deep water resolution of close to 1 arc-minute (which is also comparable to the spatial resolution of the deep water GEBCO bathymetry).Detailed inundation modelling, on the other hand often requires higher spatial resolution, in order to capture fine-scale topographic effects as well as wave runup (see the previous section).In this study the maximum resolution of "inundable areas" (any part of the domain above sea level) is set to 15 levels of refinement (or 8 arc-seconds, or ≈ 250 metres) which is comparable to the spatial resolution of the SRTM topographic dataset.Note that this choice of adaptive refinement does not rely on any a priori assumptions on the location of inundated areas and/or propagation of the tsunami wave and could thus be applied to any other event.The variable resolution is also used to implement "sponge layers" on the open boundaries of the domain.
Within bands approximately 3.65 degrees wide along all four boundaries of the domain, the spatial resolution is coarsened to only 5 levels (approximately 2.3 • resolution).The resulting numerical dissipation dampens waves before they exit the domain, thus minimising spurious wave reflections at boundaries.
With this particular choice of spatial resolutions, the global timestep is constrained by the CFL condition applied to gravity waves propagating in deep water (not by the highly-resolved waves propagating in inundated areas) and is of the order of 2 seconds.The initial evolution of the wave and mesh is illustrated in Figures 6 and 7. Figure 7 details the evolution of inundation and the associated adaptive mesh refinement in the Miyagi prefecture (Sendai) area which was highly impacted by the tsunami.The detail shown covers an area of roughly 220 × 180 km.  Figure 8 gives the corresponding evolution in the total number of grid points over 10 hours of modelled physical time.The initial number of grid points is approximately 50,000 which corresponds to the number of elements necessary to resolve the source of Figure 5 with the accuracy specified by (4).After a rapid initial increase, the number of grid points reaches a max-imum of approximately one million at t = 5 hours, the time at which waves start exiting the computational domain.This can be compared with 2 24 , the total number of grid points which would be necessary if the whole domain was resolved at a constant resolution of ≈ 1 arc-minute (the resolution necessary to model deep-water wave propagation correctly), as well as an upper bound of 2 30 ≈ one billion, if the whole domain was resolved at the maximum inundation resolution of 8 arc-seconds.A somewhat more realistic upper bound can be estimated by assuming that a static "nested mesh" approach could be used to resolve only low-lying (altitude less than e.g.50 metres) emerged land areas at the maximum resolution (8 arc-seconds) and the rest of the domain at a resolution of one arc-minute.This still gives about 50 million grid points, so that a fair estimate of the gain in number of grid points obtained with the adaptive method for this particular simulation is a factor of order 15 to 50.Note also that this factor is scale-dependent, so that as pointed out in [Popinet (2011c)] it will increase with increasing maximum spatial resolution.The gain in runtime is more difficult to estimate as it depends on the details of the respective implementations and optimisations of different methods.Previous comparisons show that the current quadtree implementation in Gerris is roughly five times slower than a comparable numerical scheme implemented on a regular Cartesian grid so that the gain in runtime would be a factor from 3 to 10.Note however that this also assumes that the speed of the regular Cartesian grid scheme would not be affected by the large memory requirements, which is unlikely in practice.The results presented here were obtained on a single-processor standard desktop PC with an average runtime of roughly two hours for one hour of physical time (this ratio is much lower in the initial phases of the simulation where few grid points are used).This runtime could be decreased by using the parallel capabilities of Gerris (including dynamic parallel load-balancing) [Agbaglah et al. (2011)].A general overview of the maximum wave elevation reached over 10 hours is given in Figures 9, 10 and 11.Maximum wave elevations in excess of 40 metres are predicted along the coast of the Oshika peninsula .These extreme values extend to the north of Ofunato, roughly between latitudes 38 • N and 39.5 • N. Large values are otherwise predicted along the Japanese coast from 35 • N to 42 • N (Figure 10).Offshore, the maximum amplitude is reached along an east-south-east direction with strong modulations due to the many small submerged and emerged topographic features in the Pacific (Figure 9).A first validation of the model is obtained by comparing the modelled wave elevation timeseries with observed sea level at DART buoys locations.Figure 12 illustrates these timeseries for all the DART buoys deployed in the simulation domain (see Figure 9) which were active during the event and which measured deviations larger than 10 cm.The DART timeseries were obtained directly from the raw data available in real-time on the DART website [NOAA (2011)] and were low-pass filtered via Fourier decomposition to remove the tidal components.The model correctly predicts both the arrival times and maximum wave elevations at all locations (keeping in mind that DART observations were not used to constrain the source model).The initial wave front is steeper than observed for most stations however.This was also observed in a previous study [Popinet (2011c)] and may be due to excessive anti-diffusion of the limiter in the shock-capturing scheme.It is also clear that the model reproduces many of the high-frequency, secondary wave reflections.For example, the clear modelled and measured signal around t = 4 hours for DART 21419 is most likely due to the reflection of the primary wave off the Emperor seamount chain (the underwater range extending south-east from DART 21416 in Figure 4).While the DART buoys data validate the long-distance, deep-water propagation, tidal gauges can give information on the inshore wave propagation.In this study, I chose to use only data available within the Global Sea Level Observing System (GLOSS) [IOC (2011b)].The GLOSS database system provides a unified way to access hundreds of sea level recording systems, in near real-time, which greatly simplifies the implementation of an automated tsunami forecasting system.Unfortunately, many existing tide gauges have not been integrated within GLOSS yet.For the area I consider only six sites are available for the Japanese islands.Another limitation of the GLOSS database is that the location of tide gauges is only recorded with an accuracy of one arc-minute.This can easily lead to tide gauges being incorrectly located in "dry" areas of the numerical model.For this study the locations of gauges have been adjusted to their most probable locations based on bathymetry and satellite images.
Figure 13 and 14 give a comparison between the numerical model and measurements for all the GLOSS stations in the simulation domain of Figure 4 which were operating at the time of the event and which measured deviations from tides larger than 10 cm.Predictions of tsunami timeseries at tidal gauge locations is notoriously more difficult than for open ocean buoys, as it is well known that measured tidal timeseries are often dominated by the resonant modes (seiches) of the small-scale bays and harbours where the gauge is located [Takahashi andAida (1963), Miller (1972)].Keeping this in mind, it is clear that the arrival time and amplitude of the first wave are very well predicted for all locations of Figure 13.This is particularly true of the stations at Hanasaki, Naha, Saipan and Pago bay and less so for Ofunato, Omaezaki, Ishigakijima and Wake island.Interestingly, the Ofunato tide gauge measured a positive first wave before ceasing transmission at t ≈ 30 minutes.The model predicts a negative first wave which is clearly related to the negative initial amplitude distribution of the source model (Figure 5).The characteristics of the following waves are very site-specific but can be roughly classified according to their dominant frequency.Some sites such as Hanasaki or Ishigakijima have dominant waves at relatively long periods ( ≈ 1 hour) which are characteristic of large-scale resonances from the coast to the edge of the continental shelf.These sites are located on open coastlines and are thus well-represented by the model and associated bathymetry.Other sites such as Naha and Omaezaki are located within harbours protected by breakwaters which are not represented at the maximum model resolution of 250 m (and in the topographic data used).As demonstrated by [Abe (2011)], for sites such as Omaezaki the detailed geometry of the harbour is crucial in determining the amplification of short-period components.This could explain the under-estimation of wave amplitudes for these sites.Finally, the sites at Wake island and Pago bay are characterised by very high frequencies (periods smaller than 15 minutes).Both the characteristic frequencies and amplitudes of the measured waves are well captured by the model.Both sites are located on small, isolated islands surrounded by deep ocean.The origin of the observed and modelled high-frequency waves is readily apparent on Figure 6 (bottom row).They are due to multiple scattering of the primary wave off the many seamounts and atolls of the north-west Pacific.These high-frequency signatures are also present (at lower amplitudes) in the modelled and measured DART buoys timeseries (Figure 12).Prediction of the extent of tsunami inundation as a function of detailed topography is an important aspect of tsunami modelling, particularly when used for urban planning and risk assessment.Extensive survey work after the Tohoku tsunami done in very difficult conditions by the Tohoku Earthquake Tsunami Joint Survey Group [TTJT (2011b), TTJT (2011a), Mori et al. (2011)] as yielded an impressive amount of data on this event: more than 5000 individual GPS records of wave height have been collected along the entire eastern Japanese coastline.In addition to this field data, various satellites have also recorded different aspects of the tsunami, including detailed maps of the extent of flooding one day after the event.Figure 15 (left) displays a detail of the topography contours (red), survey data points (blue) and extent of flooding (green) estimated using analysis of high-resolution Synthetic Aperture Radar (SAR) imaging [DLR (2011)] for the Miyagi prefecture area.
In the southern half of the map, the flooding extent is clearly controlled by the relatively steep topography, and less so for the northern half.Figure 15 (right) shows a comparison between the satellite estimation of the flooding extent on 12th March with both the modelled maximum extent of the tsunami (green) and the modelled inundated area (depth larger than 10 cm) at t = 10 hours (black).Both sets of curves are a satisfactory approximation of the satellite data, with the maximum extent curve (green) lying beyond the satellite data as expected, and the 10 cm flooding curve (black) generally underestimating slightly the extent of flooding compared with the satellite estimate.These curves are available as Google Earth KML files for the entire simulation domain of Figure 4 [ Popinet (2011b)] and can thus be compared visually with satellite images of the damaged areas.This also confirms the good model predictions of the extent of tsunami damage.Finally Figure 16 illustrates a statistical comparison of the surveyed and modelled wave elevation for all locations in the survey dataset.The survey data covers an area extending from 32 • N to 44 • N. The comparison in Figure 16 assumes that the wave elevation is mostly correlated with latitude which is broadly confirmed by the reasonably narrow range between the lower and upper quartiles (shaded areas).The model predictions are in broad agreement with the survey data, both for the qualitative wave distribution along the coastline and the elevations reached.Very good agreement is obtained between 36 and 37.5 • N as well as north of 42 • N (on the coast of Hokkaido where the Hanasaki tide gauge is located, Figure 13).The extreme wave elevations (peak values in excess of 40 metres) observed between 38.5 and 40.5 • N are also captured by the model with the correct amplitudes, but the modelled extremes occur about one degree south of the observed values: this corresponds to the Ofunato and Miyako regions respectively (see Figure 11).Both regions are characterised by very steep and indented coastlines which, together with proximity to the source, explains the high wave elevations reached there.The sources of error for both datasets (surveyed and modelled) are of course different: the model will suffer from a representation bias due for example to the inaccurate representation of steep terrain (even at the maximum resolution of 250 metres); while the surveyed data is likely to suffer from a sampling bias due to the difficulty of surveying every location and/or due to particular choices in the sampling procedure (e.g.focusing on tsunami damage to infrastructure etc...).Such a sampling bias could explain the discrepancy between modelled and surveyed data for the Miyagi area (corresponding to the inundation maps of Figure 15).The surveyed data for this area seems to be biased low (with many samples even reporting negative wave elevations) while the modelled elevation distribution seem to be more consistent with the observed extent of inundation.
Model or survey biases alone cannot however explain the large latitudinal shift in maximum wave elevation between the Ofunato and Miyako regions.The modelled distribution is consistent with the source model: the Oshika peninsula and Ofunato where the maximum wave elevations are predicted lie on a line perpendicular to the fault direction and going through the maximum of initial wave amplitude (Figure 5).In contrast the Miyako area lies at the edge of the area of maximum deformation.A possible explanation for the discrepancy is thus that the source model is locally inaccurate and that a significant fraction of the rupture energy travelled in a more northerly direction toward the Miyako coastline.Assuming that the total energy released should be unchanged, this would also divert some of the energy away from the Ofunato area.The short record from the Ofunato tide gauge (top graph of Figure 13) also suggests that the source model may be locally inaccurate: the observed first wave is positive whereas the modelled first wave is negative (with an amplitude of − 5 metres) which is consistent with the large negative amplitude area of the source model west of the epicentre.To try to clarify this issue, the model output is compared in Figures 17 and 18 with wave height records from a further six wave buoys and two bottom pressure gauges located a few tens of kilometres offshore the Sanriku coastline (locations indicated on Figure 5).An obvious discrepancy is that all the modelled timeseries show a large negative first (or second) wave which is absent in the field records.As noted before for the Ofunato tide gauge, this is consistent with the large negative amplitude in the source model (south of TM1 in Figure 5).For most sites, the peak amplitude of the largest positive wave (and subsequent waves) is correctly captured by the model with the notable exceptions of buoys 804, 803 and 801.The peak wave amplitude is under-estimated by a factor of two for buoy 804 and over-estimated by a similar factor at buoys 803 and 801.This is consistent with the bias in predicted wave elevations of Figure 16: buoys 803 and 801 are located offshore the Ofunato region (positive bias) whereas buoy 804 is located offshore the Miyako region (negative bias).This confirms that the biases in predicted wave elevations on the coastline are likely due to local errors in the source model rather than an inaccurate representation of the dynamics of flooding on this steep terrain.
Note also that several other source models proposed for this event suffer from a similar north-south bias.Using for example the source model proposed by Fuji et al [Fuji et al. (2011)], based on inversion of tide and pressure gauges and wave buoys (including all the sites of Figures 17 and 18) still leads to a large under-prediction of the wave amplitude and runup in the Miyako region (but does reduce the amplitude in the Ofunato region).Petukhin et al also report a very similar bias for their source model based purely on seismic inversion (see Figure 3 of [Petukhin et al. (2011)]) and suggest that this may be due to an inaccurate representation of the bathymetry for this area.Based on the results presented here, I do no think that this is a sufficient explanation since both the northern and southern parts of the Sanriku coast are described with similar datasets and do not differ markedly in term of complexity.The discrepancy is more likely due to the interaction of the main wave -which looks correctly described given the agreement with far-field measurements -with the initial wave field between the epicentre and the coast.

Conclusion
The quadtree-adaptive Saint-Venant solver implemented within Gerris was shown to provide robust and accurate solutions for the full range of wave phenomena associated with the Tohoku tsunami.The adaptive spatial resolution ranged from approximately 250 metres to 250 km for a domain spanning 73 • of latitude and longitude.The maximum spatial resolution proved adequate to resolve inundation processes in most areas while adaptivity allowed to track deep-water wave fronts at a typical resolution of 1 arc-minute ( ≈ 1.8 km) which was sufficient to ensure good long-distance wave energy conservation.A conservative estimate of the gain in computational efficiency obtained with the adaptive model is a factor of ten compared to a discretisation using a spatially-variable but constant-in-time resolution.This estimate increases to several orders of magnitude when compared to simple, spatially-constant (e.g.Cartesian) discretisations.Furthermore, as pointed out in [Popinet (2011c)] this gain is scale-dependent and will increase with maximum spatial resolution.The criteria used for adaptive mesh refinement are simple and do not rely on a priori assumptions about the solution.They are thus applicable to any tsunami event.
As in most tsunami models, the quality of the initial vertical displacement field is crucial to the accurate prediction of the subsequent wave evolution.The source model of Shao, Li and Ji was shown to lead to remarkably accurate predictions of initial wave amplitudes and arrival times as well as subsequent multiple reflections from bathymetry for all the DART buoys located within the simulation domain.Shao, Li and Ji only used information on seismic waves for their source inversion, so that comparisons between the modelled and observed surface waves are a true independent validation of both source and tsunami model consistency.This study demonstrates that accurate and consistent source models for tsunami wave propagation can be constructed using seismic waves only.Coupled with the high level of automation readily possible with the adaptive model, this opens the door to realtime forecasting of tsunami propagation and inundation, where the main time constraint is related to acquisition and processing of the fast (a few km/s) seismic waves signals rather than the slower (a few hundreds of m/s) surface gravity waves signals.
The results for deep-water waves are confirmed by comparisons with in-shore GLOSS stations timeseries.All stations show satisfactory agreement when allowing for the well-known dependencies of in-shore signals on (unresolved) topographic details such as breakwaters.Model results also compare well with satellite estimates of the extent of flooding as well as with extensive point survey data along the eastern Japanese coastline.However there is a significant discrepancy between surveyed and modelled wave elevations on the Sanriku coastline, where extreme wave elevations (up to 40 metres) were observed.These extremes values are reproduced in the model but further south than observed.Analysis of further inshore wave buoys and pressure gauge records indicate that this discrepancy is most probably due to an incorrect incoming wave rather than inaccurate modelling of inundation on the complex and steep shoreline.A similar discrepancy has been noted by other authors, using different source models [Fuji et al. (2011), Petukhin et al. (2011)].This may indicate that our understanding of the local seafloor deformations during the Tohoku earthquake is incomplete.A potential mechanism which could explain the discrepancy is the occurence of local submarine landslides (i.e.plastic rather than elastic deformations) triggered by the earthquake.

Figure 1 .
Figure 1.An example of quadtree discretisation (left) together with its logical representation (right).The level of the cells in the tree are also given.

Figure 2 .
Figure2.Topography, analytical (lines) and numerical (symbols) wave profiles for the times indicated in the legend.The numerical results are for l max = 14 i.e. a maximum spatial resolution of ≈ 3.6 metres.The numerical solution is represented using only one in every six discretisation points for clarity.

Figure 3 .
Figure 3. Convergence of the average and maximum errors with spatial resolution.The errors are for the predicted wave profile at t = 220 sec.

Figure 4 .
Figure 4. Extent of the simulation domain and bathymetry (altitude in metres).DART buoys (black labels) and GLOSS stations (red labels).

Figure 5 .
Figure 5.Initial vertical displacement (metres) derived from the subfaults model of Shao, Li and Ji, UCSB (Preliminary model III).The location of wave buoys (801-807) and bottom pressure gauges records (TM1 and TM2 located between TM1 and 802) used in this study are also indicated.

Figure 7 .
Figure 7. Detail of the wave evolution in the Miyagi prefecture area.The black line in all figures is the coastline.Rows from top to bottom, t = 1, 2, 4 hours.Left column: evolution of wave elevation (colorbar in metres).Right column: corresponding evolution of spatial resolution, yellow 12 levels ≈ 1 arc-minute, dark red 15 levels ≈ 8 arc-seconds ≈ 250 metres.

Figure 8 .
Figure 8. Evolution with physical time of the total number of grid points.

Figure 9 .
Figure 9. Maximum wave elevation reached over 10 hours.The (logarithmic) colorbar gives the maximum wave elevation in metres.The low values close to the boundaries of the domain are due to the mesh coarsening used to limit spurious wave reflections.

Figure 10 .Figure 11 .
Figure 10.Detail of the maximum wave elevation reached over 10 hours.The (logarithmic) colorbar gives the maximum wave elevation in metres.The model predict maximum wave elevations of more than 40 metres on the coast of the Oshika peninsula (south of Ofunato).

Figure 12 .
Figure 12.Timeseries of observed and modelled wave height at DART buoys locations, ordered by increasing arrival time from top to bottom.Horizontal axis is time in hours.

Figure 13 .Figure 14 .
Figure 13.Measured and modelled timeseries at GLOSS stations locations.Horizontal axis is time in hours.Distances from the epicentre are indicated in the legends.

Figure 15 .
Figure 15.Maps of inundation extent for the Miyagi prefecture area.Left: 0 and 10 metres contours (red), SAR estimate of inundation extent on 12 March 2011 (green), field survey data points (blue).Right: SAR inundation estimate (green), maximum modelled inundation extent (red) and modelled inundated areas at t = 10 hours (black).

Figure 16 .
Figure 16.Wave elevation reached as a function of latitude.For each dataset (modelled and surveyed) the curve is the median elevation in each 0.2 • -wide bin, and the shaded areas delimit the lower and upper quartiles of the wave elevation distribution.The vertical lines indicate the different regions displayed in Figure 11.

Figure 17 .
Figure 17.Comparison between recorded and modelled wave amplitudes for the inshore locations indicated in Figure 5.The stations are ordered by latitude from top to bottom.The horizontal axis is time in minutes.

Figure 18 .
Figure 18.Comparison between recorded and modelled wave amplitudes for the inshore locations indicated in Figure 5 (continued).The stations are ordered by latitude from top to bottom.The horizontal axis is time in minutes.