1 Introduction
Many steps are required in modeling a tsunami arising from an asteroid impact in the ocean. The impact itself forms a crater that drives the eventual tsunami creation. Modeling this requires a complex threedimensional multiphysics hydrocode, since there are many physical processes and time scales. Once the tsunami has formed, it propagates hundreds or thousands of kilometers across the ocean. When the shoreline is reached, the ultimate goal is modeling the inundation risk to coastal populations and important infrastructure at a much smaller spatial scale (typically 10 meters or less).
This work addresses the last two steps, the longdistance propagation and coastal inundation. The goal is a highfidelity model that can accurately determine the inundation risk for particular sites using available bathymetric data sets. Since large scale ocean simulations are so computeintensive, for many tsunami modeling problems the twodimensional depthaveraged Shallow Water Equations (SWE) are used for the propagation step. These equations assume the wavelength is long relative to the depth of the ocean, as is typical for tsunamis generated by large earthquakes. However, these nondispersive equations are often insufficient for shortwavelength asteroidgenerated tsunamis, giving inaccurate results for both tsunami travel time and maximum shoreline runin, as noted in our own work [4] and in several other studies, e.g. [14, 25, 27]. This is also the case for landslidegenerated tsunamis and other short wavelength phenomena; see e.g. [9]. Shorter waves experience significant dispersion (waves with different periods propagate with different speeds), while the hyperbolic SWE are nondispersive. Dispersive depthaveraged equations can be obtained by retaining more terms when reducing from the threedimensional Euler equations to two space dimensions, giving some form of ”Boussinesq equations”. The additional terms involve higherorder derivatives (typically third order), and several different models have been proposed. When solved numerically, these equations generally require implicit methods in order to remain stable with physically reasonable time steps. By contrast, the hyperbolic SWE involve only firstorder derivatives and explicit methods are commonly used.
Our numerical model is based on the GeoClaw softwave (part of the open source Clawpack software project
[8]), which has been heavily used and wellvalidated for modeling earthquakegenerated tsunamis using the SWE [3, 10, 16]. The numerical methods used are highresolution, shockcapturing finite volume methods based on Riemann solvers, a standard approach for nonlinear hyperbolic problems [15]. In the case of GeoClaw, additional features are included to make the methods “well balanced” so that the steady state of an ocean at rest is preserved. Moreover, the shoreline is represented as an interface between wet cells and dry cells and robust Riemann solvers allow the determination of the fluxes at these interfaces. The wet/dry status of a cell can change dynamically as the tsunami advances onshore or retreats. This software implements adaptive mesh refinement (AMR), critical for solving problems with vastly different spatial scales from transocean propagation to communitylevel inundation modeling. However, the patchbased AMR algorithms are based on the use of explicit solvers, and the extension of the GeoClaw software to also work with implicit solvers for Boussinesq equations has been a major part of this project. The basic approach used can also be used more generally with the AMR version of Clawpack, which has many potential applications to other wave propagation problems when dispersive or dissipative (e.g. secondorder derivative) terms are included and implicit solvers are needed.2 Overview of Approach
We are building on the work of [12], in which the Boussinesq equations described in the next section were solved using an extension of GeoClaw, but only for the case of a single grid resolution, without the AMR capability. These equations have the form of the twodimensional SWE with the addition of “source terms” involving third order derivatives. Equations of this type can often be solved by fractional step or splitting methods: first advance the solution by solving the hyperbolic shallow water equations, and then advance the solution using terms associated with the higher order derivatives (or in the opposite order, as we have found to be advantageous). For the SWE part, we can use the standard GeoClaw solver. The third derivative terms require an implicit step as described below. We are currently using a sparse linear system solver called Pardiso [5]; several other groups have found that multigrid works nicely as well.
The difficulty in the solution algorithm comes from the combination of adaptive mesh refinement and implicit solvers. Adaptive mesh refinement is critical in bridging the scales of oceanic tsunamic propagation, which typically needs resolution on the order of kilometers, and inundation modeling, where a resolution on the order of 10 meters or less is required. The patchbased mesh refinement in GeoClaw refines in time as well as space, in order to satisfy the time step stability constraint of explicit methods. If a patch is refined in space by a factor of 4, then we also typically refine in time by the same factor (as required by the CFL condition for explicit methods) and so 4 time steps are taken for the fine patch to “catch up” to the coarse patch in time. The fine patch thus needs to interpolate ghost cells that fill out the stencil at each intermediate time step. Figure
1 indicates this schematically for refinement by 2.In our approach, the solution of the implicit equations is stored as additional elements in the solution vector, increasing the number of equations in two horizontal dimensions from 3 to 5. We also reverse the typical splitting order and perform the implicit solve first. At the initial and final times during this time step, the ghost cells are still interpolated in space from the coarse grid but do not need interpolation in time. The ghost cell values at intermediate times on the fine grid are interpolated from the coarse grid values at times
and .There are other issues to consider when combining Boussinesq and SWE in a single solver. The Boussinesq equations give a better model of dispersive waves over some regime, but lack any wave breaking mechanism. As large waves approach shore this can lead to very large magnitude solitary waves that should break. The nonlinear SWE perform better at this point; a shock wave develops that is a better representation of the turbulent bore formed by a breaking wave. Very large waves, such as those that might be formed by an asteroid impact, can undergo shoaling far out on the continental shelf and dissipate some of their energy in this manner (the van Dorn effect [13]). (Shoaling refers to the modification of wave heights when the wave enters shallower water, when the wave steepens and becomes higher frequency). So it is important to transition from using the Boussinesq equations in deep water to SWE closer to shore, based on some breaking criterion. Also the onshore flooding is well modeled by SWE, which is fortunate since the wettinganddrying algorithms of GeoClaw can then be used. In our initial work we have simply suppressed the higherorder derivative terms (switching to SWE) wherever the initial water depth was 10 m or less. A better wave breaking model that allows dynamic switching has not yet been implemented but will ultimately be incorporated.
Another issue to consider is the initial conditions for the simulation. For the asteroid impact problem, even the Boussinesq equations are not adequate to model the original generation or evolution of a deep crater in the ocean. We must start with the results of a threedimensional multiphysics hydrocode simulation and produce suitable initial conditions for the depthaveraged equations. We discuss this further in Section 4.
3 Equations and Numerical Methods
For simplicity we present the equations and algorithm primarily in one space dimension and time, since this is sufficient to illustrate the main ideas.
3.1 The Shallow Water and Boussinesq Equations
The shallow water equations can be written
(3.1) 
where = water depth, topography ( offshore), , with being sea level. Thus is the depth of water at rest. The depthaveraged horizontal velocity is , and so is the momentum, and finally is the gravitional constant. These equations are a long wave length approximation to the Euler equations, in the limit of small ocean depth relative to the wavelength of the disturbance. For earthquake generated tsunamis, a typical ocean depth might be 4 km, and a subduction zone can have a wavelength of 100 km or more, giving a small ratio.
The equations (3.1) have the form of a hyperbolic system of conservation laws with a source term in the momentum equation that is nonzero only on varying topography. The GeoClaw implementation incorporates the topography term into the Riemann solvers in order to obtain a wellbalanced method [16], which amounts to solving the equations (3.1) in the nonconservative form
(3.2) 
Peregrine [22] derived a Boussinesqtype extension on a flat bottom in the form
(3.3) 
These equations have some drawbacks however, and do not match the dispersion relation of the Euler equations as well as other models developed more recently. (For a historical review of Boussinesqtype models see [6])
Madsen and Sorenson [18] and Shaffer and Madsen [26] optimized the equations by adding a term with a parameter that could be chosen to match the water wave dispersion relation more closely. On general topography, they obtained
(3.4) 
where is the initial water depth, and matching the dispersion relation leads to an optimal .
The equations (3.4) appear to have the form of the SWE (3.2) together with source terms on the right hand side. However, the standard fractional step approach can’t be used for equations in this form because the source term involves derivatives.
These equations can be rewritten as
(3.5) 
where the differential operator is defined by
Now subtracting from both sides of the momentum equation from (3.5) gives
(3.6) 
By inverting we get
(3.7) 
where is computed by solving an elliptic system:
(3.8) 
The system (3.7) now looks like SWE plus a source term that involves only spatial derivatives.
3.2 Two dimensional versions
For completeness, we include the twodimensional version of these equations to show that they have a similar structure. In two dimensions let be the two horizontal (depthaveraged) velocities. Then the shallow water equations take the form
(3.9) 
where
(3.10) 
The Boussinesq equations of [18, 26], as used in [12], take the form
(3.11) 
where now the matrix consists of four linear differential operators:
(3.12) 
with
(3.13) 
As in 1D, in order to apply a fractional step method, we subtract from both sides of the momentum equation of (3.11) so that it becomes
(3.14) 
Inverting allows rewriting (3.11) as the SWE with a source term,
(3.15) 
where involves only spatial derivatives and is determined by solving the elliptic equation
(3.16) 
Discretizing this elliptic equation leads to a nonsymmetric sparse matrix with much wider bandwidth than the tridiagonal matrix that arises in one dimension due to the cross derivative terms. Other than the significant increase in computing time required, the fractional step and adaptive mesh refinement algorithms described below carry over directly to the two dimensional situation.
3.3 Numerics
We return to the onedimensional equations in order to describe the numerical algorithm and reformulation for patchbased adaptive refinement in space and time.
We solve the onedimensional Boussinesq equation (3.7), in which is determined as the solution to (3.8), by using a fractional step method with the following steps:

Solve the elliptic equation (3.8) for the source term . After this step the values are saved on each patch to use as boundary conditions for finer patches.

Update the momentum by solving over the time step (e.g. with Forward Euler or twostage RungeKutta method). The depth does not change in this step since there is no source term in the equation.

Take a step with the homogeneous SWE, using the results of step 2 as initial data. This step uses the regular GeoClaw software and Riemann solvers.
We solve the implicit system first and then take the shallow water step because this facilitates interpolating in time for values required on the edge of grid patches. In order to explain this in more detail, we introduce some notation for a simple case in one space dimension.
First suppose we only have a single grid at one resolution, with no AMR. We denote the numerical solution at some time by , the cellaveraged approximations to depth and momentum on the grid. We also use for the source term determined by solving the discrete elliptic system defined by on this grid. We also assume at the start of the time step that we have boundary conditions for and also for the Boussinesq correction , provided in the form of “ghost cell” values in a layer of cells surrounding the grid (or two layers in the case of since the highresolution explicit methods used have a stencil of width 5 because of slope limiters). On a single grid we assume that it is sufficient to use the Dirichlet condition in all ghost cells surrounding the grid, i.e. that there is no Boussinesq correction in these cells. This is reasonable for a large domain where the waves of interest are confined to the interior of the domain. We also use zeroorder extrapolation boundary condition for , which give a reasonable nonreflecting boundary condition for the SWE step [16].
A single time step of the fractional step algorithm on this grid then takes the following form in order to advance to at time :

Solve the elliptic system for . The righthand side depends on .

Advance the solution using the source terms (Boussinesq corrections):
(using Forward Euler, for example). 
Take a time step of length with the SWE solver, with initial data , to obtain . We denote this by .
In the software it is convenient to store the source term at each time as another component of the solution vector, so we also use to denote this full solution at time .
Now suppose we have two grid levels with refinement by a factor of 2 in time. We denote the coarse grid values at some time as above. We assume that the fine grid is at time , but that on the fine grid we must take two time steps of to reach time . We denote the fine grid values at time using lower case, and . We also need boundary conditions in the ghost cells of the fine grid patch. If a patch edge is coincident with a domain boundary, then we use the Dirichlet BC and extrapolation BCs for , as described above. For ghost cells that are interior to the coarse grid, we let represent a spatial interpolation operator that interpolates from coarse grid values to the ghost cells of a fine grid patch at time . This operator is applied to all three components of , i.e. to the source term as well as the depth and momentum, in order to obtain the necessary boundary conditions for .
Then one time step on the coarse grid, coupled with two time steps on the fine grid, is accomplished by the following steps:

Coarse grid step:

Take time step on the coarse grid as described above for the single grid algorithm, but denote the result by since these provisional values will later be updated.

Using the Dirichlet BCs on the domain boundary, solve for a provisional . This will be needed for interpolation in time when determining boundary conditions for on the fine grid, using and .


Fine grid steps:

Given and boundary conditions , solve the elliptic system for .

Update using the source terms, .

Take a shallow water step: .
Note that we use to denote the intermediate time. 
Obtain BCs at this intermediate time as .

Solve the elliptic system for .

Update using the source terms, .

Take a shallow water step: .


Update coarse grid:

Define by the provisional values by where there is no fine grid covering a grid cell, but replacing by the average of over fine grid cells that cover any coarse grid cell.

The final step is applied because the fine grid values are more accurate than the provisional coarse grid values.
We then proceed to the next coarse grid time step. Note that at the start of this step, the updated will be used to solve for . The provisional is discarded. Hence two elliptic solves are required on the coarse level each time step, rather than only one as in the single grid algorithm.
If the refinement factor is larger than 2, then the same approach outlined above works, but there will be additional time steps on level 2. For each time step the ghost cell BCs will be determined by linear interpolation in time between and .
If there are more than two levels, then this same idea is applied recursively: After each time step on level 2, any level 3 grids will be advanced by the necessary number of time steps to reach the advanced time on level 2. In this case there will also be two elliptic solves for every time step on level 2, once for the provisional values after advancing level 2, and once at the start of the next level 2 time step after on level 2 has been updated by averaging the more accurate level 3 values.
4 Computational Results
In this section we show an endtoend simulation of a hypothetical asteroid impact off the coast of Washington, from initial conditions to shoreline inundation. We present our initialization procedure in some detail. The section ends with a discussion of results.
4.1 Initialization Procedure
The computational results presented in this section use initial conditions of a static crater, illustrated in Figure 2(a). This is a standard test problem in the literature, from [27]. The crater is one kilometer deep, with a diameter of 3 kilometers, in an ocean of depth 4 kilometers. Depthaveraged equations are unsuitable for modeling the generation of the crater and the initial flow, since the ensuing large vertical velocity components are not modeled. The initial conditions for our depthaveraged simulations are taken from a threedimensional hydrocode simulation of the first 251 seconds after impact. The hydrocode ALE3D [20] was run by a collaborator [24] in a radially symmetric manner, and the surface displacement at time seconds was recorded, as shown in Figure 2(b). It proved too noisy to depthaverage the horizontal velocity from the hydrocode. Instead we set the velocity based on the surface displacement and assuming that the wave was a purely outgoing wave satisfying the SWE, for which the velocity then depends only on the ocean depth and surface displacement. We could then place the “initial” crater anywhere in the ocean.
This procedure for initialization of the velocity can be done because the wave speed for the SWE is independent of wave number, and the eigenvectors of the linearized Jacobian matrix give the relation between surface elevation and fluid velocity for unidirectional waves. However, in the Boussinesq equations the wave speed depends on wave number and using the initialization based on the SWE results in a small wave propagating inward as well. Better initialization procedures will be investigated in future research.
4.2 Adaptive simulation
As an illustration, we show a simulation of a hypothetical asteroid impact off the coast of Washington. We place the initial crater approximately 150 km west of Grays Harbor, to study the vulnerable area around Westport, WA, shown in Figure 3. This wellstudied area is in close proximity to the Cascadia Subduction Zone, which can generate Mw 9 earthquakes. The Ocosta elementary school in Westport was recently rebuilt to incorporate the first tsunami vertical evacuation structure in the U.S., due to the low topography of this region, with design work based in part on GeoClaw modeling [11]. Detailed bathymetry and topography data is available in this region at a resolution of 1/3 arcsecond [21], which is roughly 10 meters in latitude and 7 m in longitude at this location. For the ocean we use the etopo1 topography DEM [1] at a resolution of 1 arcminute.
A 7 level simulation was used, starting with a coarsest level over the ocean with arcminutes, and refining by factors 5, 3, 2, 2, 5, 6 at successive levels, with an overall refinement by 1800 for the level7 grids with arcsecond. On each grid so that the finestlevel computational grids are at a resolution of roughly 10 m in both and .
Figures 4 and 5 show snapshots of the simulation at the indicated times. The dispersion is clearly evident, with a much more oscillatory solution than would be obtained with the shallow water equations. The figures show how the fine grid patches move to follow the expanding wave. The refinement is guided to focus at later times only on waves approaching Grays Harbor. There are some reflections at the grid boundaries, but they are much smaller in magnitude than the waves we are tracking. The 6th level refined patch appears approximately 20 minutes after the impact, to track the waves as they approach Grays Harbor. Note how the bathymetry is refined along with the solution when the finer patches appear. The close up plots near Grays Harbor in Figure 5 shows “soliton fission”, a nonlinear dispersive wave phenomenon seen near the coast that can be captured with Boussinesq solvers [2, 19]. In this simulation, we switch from Boussinesq to SWE at a depth of 10 meters (based on the undisturbed water depth). Figure 5 also shows the waves sweeping over the Westport and Ocean Shores peninsulas. In this calculation, the finest level 7 grid was placed only over Westport, while the level 6 grid covers both peninsulas.
4.3 Discussion of Results
Since the wavelengths of asteroidgenerated tsunamis are shorter than those of earthquakegenerated tsunamis, dispersive effects may be very important. Dispersion gives rise to a highly oscillatory set of waves that propagate at different speeds, possibly affecting the arrival time of the first wave and leading to significant waves over a longer time period. These waves can undergo substantial amplification during the shoaling process on the continental shelf and break up into large solitary waves.
On the other hand, wave breaking can rapidly dissipate the energy in short wavelength waves. Moreover, a train of waves approaching shore results in nonlinear interactions in the swash zone, where waves run up the beach. In the swash zone, a large approaching wave may be largely negated by the rundown of the previous wave. Thus it is not clear a priori whether waves modeled with the Boussinesq equations will result in substantially different onshore inundation than would be observed is only using the SWE, for which computations are much less expensive. This question still remains to be answered.
5 Open Problems and Future Research
We have demonstrated that it is possible to develop a highfidelity modeling capability that includes propagation and inundation by combining the Boussinesq equations, shallow water equations, and patchbased adaptive mesh refinement in space and time. We are embedding this in the GeoClaw software framework, which has previously been wellvalidated for earthquakegenerated tsunamis. This will allow the efficient simulation of dispersive waves generated from asteroid impacts as they propagate across the ocean, combined with highresolution simulation of the resulting inundation on the coast. This capability could be very important in hazard assessment for an incipient impact.
Unfortunately the software is not yet robust enough for general use. While it often works well, stability issues sometimes arise at the edges of finer grid patches when refinement ratios greater than 2 are used from one level to the next. Larger refinement ratios are generally desirable, so that an overall refinement factor of several thousand between the coarsest and finest meshes can be achieved with 6 or 7 levels of mesh refinement. Other Boussinesq codes we are aware of use factor of 2 refinement [23, 7], so this may be an inherent instability. Moreover, numerical instabilities have also been observed by other researchers [17] even on a uniform grid when there are sharp changes in bathymetry rather than in the grid resolution. We are investigating this issue theoretically, and may find that a different discretization or even a different fomulation of the Boussinesq equations is required to obtain a sufficiently robust code.
We are also continuing to investigate the shoaling phenomena and the best way to incorporate wave breaking and the transition from Boussinesq to shallow water equations. This can have a significant impact on the resulting onshore runup and inundation.
Not discussed here but currently under investigation is the possibility of solving the implicit system of equations on multiple levels in a coupled manner, whenever the levels have been advanced to the same point in time. It is an open question whether a coupled system is more accurate and/or stable, and possibly less computationally expensive, in the context of patchbased adaptive mesh refinement that includes refinement in time.
Finally, as mentioned in section 4.1, additional research is needed on ways to initialize the depthaveraged model. Better initialization procedures will allow a more seamless transition from threedimentional hydrocode simulations of asteroid impacts in the ocean to our model of tsunami propagation and inundation.
Acknowledgments. We thank Darrel Robertson, a member of the ATAP team, for providing the hydrocode simulation results for our initial conditions.
This work was partially supported by the NASA Asteroid Threat Assessment Project (ATAP) through the Planetary Defense Coordination Office and BAERI contract AO9667.
References
 [1] C. Amante and B. W. Eakins. ETOPO1 1 ArcMinute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC24, http://www.ngdc.noaa.gov/mgg/global/global.html, 2009.
 [2] Toshitaka Baba, Narumi Takahashi, Yoshiyuki Kaneda, Kazuto Ando, Daisuke Matsuoka, and Toshihiro Kato. Parallel implementation of dispersive tsunami wave modeling with a nesting algorithm for the 2011 Tohoku tsunami. Pure and Applied Geophysics, 172(12):3455–3472, 2015.
 [3] M. J. Berger, D. L. George, R. J. LeVeque, and K. T. Mandli. The GeoClaw software for depthaveraged flows with adaptive refinement. Advances in Water Resources, 34(9):1195–1206, September 2011.
 [4] M. J. Berger and R. J. LeVeque. Modeling issues in asteroidgenerated tsunamis. NASA Technical Memorandum NASA/CR2018219786, ARCEDAATN53167, http://hdl.handle.net/2060/20180006617, 2018.
 [5] Matthias Bollhöfer, Olaf Schenk, Radim Janalik, Steve Hamm, and Kiran Gullapalli. Stateoftheart sparse direct solvers. In Ananth Grama and Ahmed H. Sameh, editors, Parallel Algorithms in Computational Science and Engineering, pages 3–33, Cham, 2020. Springer International Publishing.
 [6] M. Brocchini. A reasoned overview on Boussinesqtype models: the interplay between physics, mathematics and numerics. Phil. Trans. Royal Soc. A, 469(20130496), 2013.
 [7] D. Calhoun and C. Burstedde. Forestclaw: A parallel algorithm for patchbased adaptive mesh refinement on a forest of quadtrees, 2017. ArXiv:1703.03116.
 [8] Clawpack Development Team. Clawpack software. http://www.clawpack.org, 2021.
 [9] S. Glimsdal, G. K. Pedersen, C. B. Harbitz, and F. Løvholt. Dispersion of tsunamis: Does it really matter? Nat. Hazards Earth Syst. Sci., 13:1507–1526, 2013.
 [10] F. González, R. J. LeVeque, J. Varkovitzky, P. Chamberlain, B. Hirai, and D. L. George. GeoClaw Results for the NTHMP Tsunami Benchmark Problems. http://depts.washington.edu/clawpack/links/nthmpbenchmarks, 2011.
 [11] F. I. Gonzaléz, R. J. LeVeque, and L. M. Adams. Tsunami Hazard Assessment of the Ocosta School Site in Westport, WA. http://hdl.handle.net/1773/24054, 2013.
 [12] J. Kim, G. K. Pedersen, F. Løvholt, and R. J. LeVeque. A Boussinesq type extension of the GeoClaw model  a study of wave breaking phenomena applying dispersive long wave models. Coastal Engineering, 122:75 – 86, 2017.
 [13] D. G. Korycansky and P. J. Lynett. Offshore breaking of impact tsunami: The Van Dorn effect revisited. Geophysical Research Letters, 32(10), 2005.
 [14] D. G. Korycansky and Patrick J. Lynett. Runup from impact tsunami. Geophysical Journal International, 170:1076–1088, 2007.
 [15] R. J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, 2002.
 [16] R. J. LeVeque, D. L. George, and M. J. Berger. Tsunami modelling with adaptively refined finite volume methods. Acta Numerica, 20:211 – 289, May 2011.
 [17] F. Løvholt and G. Pedersen. Instabilities of Boussinesq models in nonuniform depth. International Journal for Numerical Methods in Fluids, 61:606–637, 2009.
 [18] P. A. Madsen and O. R. Sørensen. A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowlyvarying bathymetry. Coastal engineering, 18(34):183–204, 1992.
 [19] Masafumi Matsuyama, Masaaki Ikeno, Tsutomu Sakakiyama, and Tomoyoshi Takeda. A study of tsunami wave fission in an undistorted experiment. Pure and Applied Geophysics, 164(2):617–631, 2007.
 [20] A. L. Nichols and D. M. Dawson (Eds.). ALE3D User’s Manual: An Arbitrary Lagrange/Eulerian 2D and 3D Code System. Lawrence Livermore National Laboratory LLNLSM726137, 2017.
 [21] NOAA National Geophysical Data Center. Astoria, Oregon 1/3 Arcsecond MHW Coastal Digital Elevation Model. https://www.ncei.noaa.gov/access/metadata/landingpage/bin/iso?id=gov.noaa.ngdc.mgg.dem:5490, Accessed 2019.
 [22] D. H. Peregrine. Long waves on a beach. Journal of Fluid Mechanics, 27(4):815–827, 1967.
 [23] S. Popinet. A quadtreeadaptive multigrid solver for the SerreGreenNaghdi equations. J. Comp. Phys., 302:336–358, 2015.
 [24] D. Robertson. Personal Communication, 2019.
 [25] Darrel K. Robertson and Galen R. Gisler. Near and farfield hazards of asteroid impacts in oceans. Acta Astronautica, 156:262–277, 2019.
 [26] H. A. Schäffer and P. A. Madsen. Further enhancements of Boussinesqtype equations. Coastal Engineering, 26(12):1–14, 1995.
 [27] S. N. Ward and E. Asphaug. Asteroid impact tsunami: A probabilistic hazard assessment. Icarus, 145:64–78, 2000.
Comments
There are no comments yet.