Air-burst Generated Tsunamis

TL;DR

This study investigates whether air-bursting small asteroids over water can generate tsunamis, finding that such events pose minimal tsunami threat based on numerical simulations and analytical models.

Contribution

The paper introduces a numerical and analytical analysis of air-burst generated tsunamis, highlighting their negligible threat compared to other tsunami sources.

Findings

Numerical simulations show very little tsunami generated by air bursts.

An explicit solution model explains the phenomena observed in simulations.

More sophisticated models suggest compressibility and dispersion effects are minor.

Abstract

This paper examines the questions of whether smaller asteroids that burst in the air over water can generate tsunamis that could pose a threat to distant locations. Such air burst-generated tsunamis are qualitatively different than the more frequently studied earthquake-generated tsunamis, and differ as well from impact asteroids. Numerical simulations are presented using the shallow water equations in several settings, demonstrating very little tsunami threat from this scenario. A model problem with an explicit solution that demonstrates and explains the same phenomena found in the computations is analyzed. We discuss the question of whether compressibility and dispersion are important effects that should be included, and show results from a more sophisticated model problem using the linearized Euler equations that begins to addresses this.

Equations

$$\begin{array}{rl}{h}_t+(hu{)}_x=& 0\\ (hu{)}_t+(h{u}^2+\frac{1}{2}g{h}^2{)}_x=& \frac{-h{p_e}_x}{{\rho }_w},\end{array}$$

$$m=x-st.$$

$$-s{h}_m+(hu{)}_m$$

$$-s(hu{)}_m+(h{u}^2+\frac{1}{2}g{h}^2{)}_m$$

$$u(m)=\frac{s(h(m)-h_0)}{h(m)}.$$

$$-s(s(h-h_0){)}_m+{\left(\frac{s^2(h-h_0{)}^2}{h}\right)}_m+{\left(\frac{1}{2}g{h}^2\right)}_m=\frac{-h{p_e}_m}{{\rho }_w}.$$

$$\frac{s^2}{2}{\left(\frac{h_0^2}{h^2}\right)}_m+g{h}_m=\frac{-{p_e}_m}{{\rho }_w}.$$

$$\frac{s^2}{2}\left(1-\frac{h_0^2}{h(m{)}^2}\right)+g{h}_0\left(1-\frac{h(m)}{h_0}\right)=\frac{p_e(m)}{{\rho }_w}.$$

$${\left(\frac{h_0}{h(m)}\right)}^2=\frac{h_0^2}{{\left(h_0+h_r\right)}^2}\approx 1-\frac{2h_r}{h_0},$$

$$h_r=\frac{h_0{p}_e}{{\rho }_w(s^2-c_w^2)}.$$

$$x(t)=\frac{A}{1-{\omega }^2}\cos (\omega t).$$

$$p_e=p_{ambient}\exp (-0.1(x-st{)}^2),$$

$$h(x,t)=h_r(x-st)-\left(\frac{s}{c}+1\right)\frac{h_r(x-ct)}{2}+\left(\frac{s}{c}-1\right)\frac{h_r(x+ct)}{2},$$

$$\begin{array}{rl}{\tilde{\rho }}_t+{\rho }_w{\tilde{u}}_x+{\rho }_w{\tilde{w}}_z& =0\\ {\rho }_w{\tilde{u}}_t+c_a^2{\tilde{\rho }}_x& =0\\ {\rho }_w{\tilde{w}}_t+c_a^2{\tilde{\rho }}_z& =-\tilde{\rho }g,\end{array}$$

$$h(x,t)=h_0+h_r(x,t).$$

$$\text{bottom:}\tilde{w}(x,z=0,t)$$

$$\text{top:}\frac{\partial h_r(x,t)}{\partial t}$$

$$\text{pressure bc:}{c}_a^2\tilde{\rho }(x,h_0,t)-{\rho }_{w}g{h}_r(x,t)$$

$$-s{\tilde{\rho }}_m+{\rho }_w{\tilde{u}}_m+{\rho }_w{\tilde{w}}_z$$

$$-s{\rho }_w{\tilde{u}}_m+c_a^2{\tilde{\rho }}_m$$

$$-s{\rho }_w{\tilde{w}}_m+c_a^2{\tilde{\rho }}_z$$

$$\tilde{w}(m,0)$$

$$\tilde{w}(m,h_0)$$

$$c_a^2\tilde{\rho }(m,h_0)$$

$$p_e(m)=A_k{e}^{ikm},$$

$$h_r(m)$$

$$\tilde{\rho }(m,z)$$

$$\tilde{u}(m,z)$$

$$\tilde{w}(m,z)$$

$$-s{\rho }_w\tilde{u}+c_a^2\tilde{\rho }=0.$$

$$\tilde{\rho }=\frac{s{\rho }_w\tilde{u}}{c_a^2},$$

$${\tilde{u}}_m\left(1-\frac{s^2}{c_a^2}\right)+{\tilde{w}}_z$$

$$-{\tilde{w}}_m+{\tilde{u}}_z$$

$${\left(\begin{array}{c}\hat{u}\\ \hat{w}\end{array}\right)}_z=\left(\begin{array}{c}-g/c_a^2\hat{u}+ik\hat{w}\\ -ik\hat{u}(1-s^2/c_a^2)\end{array}\right)=\left[\begin{array}{cc}-g/c_a^2& ik\\ -ik(1-s^2/c_a^2)& 0\end{array}\right]\left(\begin{array}{c}\hat{u}\\ \hat{w}\end{array}\right).$$

$$\left(\begin{array}{c}\hat{u}\\ \hat{w}\end{array}\right)=a_{+}{\mathbf{v}}_{+}{e}^{{\mu }_{+}z}+a_{-}{\mathbf{v}}_{-}{e}^{{\mu }_{-}z},$$

$${\mu }_{\pm }=\frac{\frac{-g}{c_a^2}\pm \sqrt{\frac{g^2}{c_a^4}+4k^2(1-s^2/c_a^2)}}{2}.$$

Full text

Keywords: tsunami; asteroid-generated air-burst; shallow water equations;

         linearized Euler equations.

1 Introduction

In Feb. 2013, an asteroid with a 20 meter diameter burst 30 km high in the atmosphere over Chelyabinsk, causing substantial local damage over a 20,000 $\text{km}^{2}$ region [19]. The question arises, what would be the effect of an asteroid that bursts over the ocean instead of land? The concern is that the atmospheric blast wave might generate a tsunami threatening populated coastlines far away.

There is little literature on air-burst-generated tsunamis. Most of the literature on asteroids study the more complicated case of water impacts, where the meteorite splashes into the ocean [23, 6, 5]. This involves much more complicated physics. The only reference we are aware of that relates to a blast-driven water wave is from the 1883 volcanic explosion of Krakatoa [8]. The authors report a tide gauge in San Francisco registered a wave that could not be explained by a tsunami. There is also some analytic work in [9], where they derive asymptotic formulas for water waves from explosions and from initial cavities. There is more literature on meteo-tsunamis. These are also driven by air-pressure events and have similarities to our case, but occur in a different regime of air speed and water depth.

This paper studies the behavior of air-burst generated tsunamis, to better understand the potential threat. In the first part of the paper, we present simulations under a range of conditions using the shallow water equations and the GeoClaw software package [2]. We compute the ocean’s response to an overpressure as calculated in [1]. The overpressure was found by simulating the blast wave in air, and extracting the ground footprint. Roughly speaking, the blast wave model corresponds to the largest meteor that deposits all of its energy in the atmosphere without actually reaching the water surface. With this forcing, if there is no sizeable response then we can conclude that that air-bursts do not effectively transfer energy to the ocean, and there is little threat of distant inundation.

Typically the shallow water equations are used for long-distance propagation, since they efficiently and affordably propagate waves over large trans-oceanic distances. Other alternatives, such as the Boussinesq equations, are much more expensive, and at least for the case of earthquake-generated tsunamis the difference seems to be small [12].

In general, our results using the shallow water equations suggest that air-burst generated tsunamis are too small to cause much coastal damage. Of course, depending on local bathymetry there could an unusual response that is significant. For example, Crescent City, California is well-known to be subject to inundation due to the configuration of its harbor and local bathymetry. However, we find that to generate a large enough response so that the water floods the coastline, the blast has to be so close that the blast itself is the more dangerous phenomenon. This is also the conclusion reached by Gisler et al. [6] and Melosh [14] for the case of asteroid water impacts.

In the second part of this paper, we study model problems to better understand and describe the phenomena we compute in the first part. The first model problem is based on the one-dimensional shallow water equations for which we can obtain an explicit closed form solution. It assumes a traveling wave form for the pressure forcing. Actual blast waves only approximately satisfy this hypothesis for a short time before their amplitudes decay. Nevertheless, the model explains several key features that we observe in the two-dimensional simulations. We observe a response wave that moves with the speed of the atmospheric forcing. There is also the gravity wave, or tsunami, moving at the shallow water wave speed, that is generated by the initial transient of the atmospheric forcing. We study in detail the response wave, or ‘forced’ wave, but the two are closely related. The analysis shows that the forced wave is proportional to the local depth of the water at each location, a phenomena clearly seen in our computations. The model problem also allows us to assess the importance of nonlinear modeling. For most physical situations related to air-burst tsunamis, the linear and nonlinear models give similar predictions.

In our final section we assess the effect of corrections to the shallow water equations arising from compressibility and dispersion using a second model problem - the linearized Euler equations. Air bursts have a much shorter time scale than earthquake-generated tsunamis, comparable to the acoustic travel time to the ocean floor. This leads to the question of whether compressibility of the ocean water could be a significant factor. In addition, air bursts have much shorter wavelengths, on the order of 10 to 20 kilometers, at least for meteors with diameter less than 200 meters or so. Recall that the shallow water model results from assuming long wavelengths and incompressibility of the water. Our results show that for air-burst generated tsunamis, dispersion can be significant but that compressibility is less so, suggesting interesting avenues for future work.

This work is an outgrowth of the 2016 NASA-NOAA Asteroid-generated Tsunami and Associated Risk Assessment Workshop. The workshop conclusions are summarized in [16]. Several other researchers also performed simulations, and videos of all talks are available on-line222All presentations are available at https://tsunami-workshop.arc.nasa.gov/workshop2016/sched.php.

2 Two-dimensional Simulations

In this section we present results from two sets of simulations. We use a 250MT blast, which roughly corresponds to a meteor with a 200 meter diameter entering the atmosphere with a speed of 20 km/sec. Generally speaking, this is the largest asteroid that would not splash into the water.333 Initially we used a blast wave corresponding to a 100MT blast, but since no significant response was found we do not include those results here. We also did simulations where we increased the pressure forcing by a factor of 1.2 with no change to the conclusions. For each location, we did several simulations varying the blast locations with no meaningful difference in results, so we only present one representative computation in each set of simulations.

In the first set of results, we locate the blast in the Pacific about 180 kilometers off the coast of near Westport, Washington. This spot was chosen since it is well studied by the earthquake-generated tsunami researchers due to its proximity to the M9 Cascadia fault [17, 4]. By the time the waves reach shore they have decayed and are under a meter high. Since they do not have the long length scales of earthquake tsunamis, we not see any inundation on shore. In the second set of results we move the location offshore to Long Beach, California, where there is significant coastal infracture, and has also been studied extensively in relation to earthquake tsunamis [21]. We place the blast approximately 30 kilometers from shore, so that there is less time for the waves to decay. In all simulations, bathymetry is available from the NOAA National Center for Environmental Information web site.

To perform these simulations, we use a model of the blast wave simulated in [1]. The ground footprint for the overpressure was extracted, a Friedlander profile was fit to the data, and its amplitude as a function of time was modeled by a sum of Gaussians. The 250MT model is shown in Fig. 1, with a few of the profiles drawn in black to illustrate their form. The profiles start with the rise in pressure from the incoming blast, and are followed by the expected rarefaction wave (underpressure) some distance behind. Note that the maximum amplitude is over 4 atmospheres, but decays rapidly from its initial peak.

In the model, the blast wave travels at a fixed speed of 391 m/sec. This may be less accurate at early times. If the asteroid enters at a low angle of incidence, the blast wave travels more quickly when it first hits the ground. This would also lead to a more anisotropic response when see from the ground. Here, however, we assume the blast wave is radially symmetric. We then use this model of the overpressure as a source term in our two-dimensional shallow water simulations using the software package GeoClaw.

GeoClaw is an open source software package developed since 1994 [11] for modeling geophysical flows with bathymetry using the shallow water equations. It is mostly used for simulations of tsunami generation, propagation and inundation. GeoClaw uses a well-balanced, second-order finite volume scheme for the numerics [10, 3]. Some of the strengths of GeoClaw include automatic tracking of coastal inundation, robustness in its handling of dry states, a local adaptive mesh refinement capability, and the automated setup that allows for multiple bathymetry input files with varying resolution. A bottom friction term is included using a constant Manning coefficient of 0.025. The results below do not include a Coriolis force, which we have found to be unimportant. There is no dispersion in the shallow water equations. In 2011 the code was approved by the U.S. National Tsunami Hazard Mitigation Program (NTHMP) after an extensive set of benchmarks used to verify and validate the code [7].

2.1 Westport Results

For this set of simulations, the 250MT blast is located at $-126.25^{\circ}$ longitude and $46.99^{\circ}$ latitude, about 30 kilometers from the continental slope. The ocean is 2575 meters deep at this spot. The blast location is about 180 kilometers from shore. Many simulations were performed with different mesh resolutions. The finest grids used in the adaptive simulations had a resolution of 1/3 arc second. Three bathymetric data sets were used – a 1 minute resolution covering the whole domain, a 3 second resolution nearer shore, and a 1/3 arc second bathymetry that included the shoreline itself.

In Fig. 2 we show a Hovmöller plot through the center of the blast location at a fixed latitude. On the left is the atmospheric overpressure for the first 300 seconds. This is the forcing that travels at 391 m/sec. On the right is the amplitude of the water’s response. Two waves traveling at different speeds are visible. A shallow water gravity wave travels with speed $\sqrt{gh}$, which at the blast location is 158 m/sec. It is evident that the blast wave travels approximately twice as fast as the gravity wave. The blast wave reaches the edge of the graph in just over 150 seconds instead of the 300 seconds of the main water wave (in blue, since it is a depression). Also visible in the wave height plot is a wave that starts off in red and travels at the same speed as the blast wave, and whose amplitude decays more rapidly. Here the color scale saturates below the maximuym the maximum value in ecah plot so that smaller waves are visible.

Fig. 3 shows the maximum amplitude found at any time in the simulation at that location. Note that the color bar is not linear in this plot, so that the different levels can more easily be seen. Nearest the blast location the maximum wave amplitude is over 10 meters, but it decay rapidly. As the waves approach shore, the waves are amplified in a non-uniform way by the bathymetry. The coastline is outlined in black. The light gray contour line represents the location of the waves after 30 minutes. We do not see any inundation of land, although admittedly at this resolution it would be hard to see.

Fig. 4 shows the time history of wave heights through several gauge plots. The left plot shows 7 gauges placed $0.1^{\circ}$ apart (about 10 kilometers at this latitude), starting about 1 kilometer from the blast. The gauge closest to the blast location has a maximum amplitude that reaches 5 meters. Subsequent gauges show a very rapid decay in maximum amplitude. These positive elevation waves are the water’s response to the blast wave overpressure, and travel at the same speed as the blast wave. Most of the ocean’s response at this location appears as a depression, not an elevation. The negative amplitude wave travels at the gravity wave speed, $\sqrt{gh}$, where the water has depth $h$. It shows much less decay in amplitude. For example, looking at gauge 3 and 5, the peak amplitude decays from 2.7 meters to 0.72 meters in about 50 seconds, whereas between 100 and 200 seconds, the trough decays from -5.3 meters to -4.1 meters in about 100 seconds.

Fig. 4 right shows gauges approaching the shore, starting about 100 kilometers away from the blast. These are not equally spaced but are placed from .25 to $.1^{\circ}$ apart (from 25 to 10 kilometers at this latitude), becoming closer as they approach shore and the bathymetry changes more rapidly. Shoaling is observed as the wave amplitudes increase, seen in gauges 17 and higher. The maximum elevation is between 0.5 and 1 meters, and its duration is short, at least compared to earthquake-generated tsunamis.

Finally, Fig. 5 shows several close-ups of the region near shore. The waves are of uneven strength due to focusing from the bathymetry. The maximum amplitude is around 1 meter. The sequence shows waves reflecting from the coastline but not flooding it. Some waves enter Grays Harbor, but they small amplitude and do not flood the inland area either.

2.2 Long Beach Results

For this set of experiments we move the simulations to Long Beach, California. We locate the blast very close to shore so that the waves do not have time to decay. Again we have detailed bathymetry at a resolution of 1/3 arc second between Catalina Island and Long Beach, and use a 1 minute dataset outside of this region. The blast is located at $-118.25^{\circ}$ longitude and $33.41^{\circ}$ latitude, where the ocean is 797 meters deep. This is about 30 kilometers from shore. Fig. 6 shows the region where the blast is located, and a zoom of the Long Beach harbor where we will look for flooding.

Fig. 7 shows the ocean response at several points in time. The black circle on the plots indicates the location of the air blast. In the plot marked at 25 seconds, note how the wave height in red that is closest to the blast location is not as circular as the air blast itself. It also decays faster than in the Westport computation. This will be explained by the model problem presented in the next section.

There is a breakwater that protects long Beach. It reflects most of the waves that reach it, with only a small portion getting through the opening. Waves that reach the harbor go around the breakwater, and are reflected from the shoreline back into this region.

Fig. 8 shows a plot of the maximum water amplitudes seen in the harbor area. We do see some overtopping of land, but it is very small. In several locations it reaches 0.5 meters, where the inlet exceeds its boundaries, and on the dock in the middle. The region with the largest accumulation is just outside the harbor before the breakwater, where the maximum amplitude seen is between 3 and 6 meters. There is a steep cliff here however and the water does not propagate inland. Paradoxically, in other experiments where the blast was located closer to shore by a factor of 2, there was no overtopping. This can also be explained by our model problem in the next section.

3 Shallow water model

In this section we present a one-dimensional model of the shallow water equations (SWE) that explains much of the behavior seen in the previous examples. In the SWE, the atmospheric overpressure appears as an external forcing $p_{e}$ in the momentum equation. In one space dimension it is

[TABLE]

where $h$ is the height of the water surface over the bottom, $u$ is the depth-averaged velocity of the water in the $x$ direction, and $g$ is gravity, and $\rho_{w}$ is the density of water. ${p_{e}}$ is the external pressure forcing, and it’s $x$ derivative is ${p_{e}}_{x}$. We assume constant bathymetry in the model. See [22] for these equations, or [13] for a complete derivation. In this section (and the next), the conclusions are in the last few paragraphs after the analysis.

3.1 Derivation and Analysis

As stated in the introduction, we simplify the pressure forcing by assuming it has the form of a traveling wave, and look for solutions $h$ and $u$ that are traveling waves too. This means they are functions only of the moving variable

[TABLE]

so that $\partial_{t}h=-s\partial_{m}h(x-st)=-sh_{m}$. The equations (1) become a pair of ordinary differential equations

[TABLE]

Equation (2) can be integrated to give $-sh+hu=\mbox{\em const}$. We evaluate the constant by taking $m\to\infty$, where $u\to 0$ and $h\to h_{0}$, with $h_{0}$ the undisturbed water height. (We assume the overpressure has localized support, and goes to zero as $m\rightarrow\infty$.) Therefore $-sh+hu=-sh_{0}$. This may be re-written as

[TABLE]

We use this to eliminate $u$ from (3), which gives

[TABLE]

After some algebra, this leads to

[TABLE]

As before, this may be integrated exactly. Again we use the boundary conditions $h\to h_{0}$, $u\to 0$, and $p_{e}\to 0$ as $m\to\infty$. The result is

[TABLE]

To summarize, eq. (5) is the water’s response according to shallow water theory. The solution of the differential equation system (2) and (3) is an algebraic relation between the overpressure and the response height. It tells us that the water height at a point $m=x-st$ is determined by the overpressure at the same point.

To get a better feel for the behavior of the solution (5) we linearize it, writing $h(m)=h_{0}+h_{r}(m)$ where $h_{r}(m)$ is the response height. The linearization uses the relation

[TABLE]

which is valid when $h_{r}\ll h_{0}$. This is our case, since the change in wave height $h_{r}$ is a number in meters where $h_{0}$ is typically measured in kilometers. The linearization of (5) is:

[TABLE]

Fig. 9 shows that the full response theory (5) and the linear approximation (6) are very close to each other for the parameters of interest. The plot uses a constant depth of $h_{0}=4$ km, and takes $\rho_{w}=1025\,\mbox{kg}/\mbox{m}^{3}$. The maximum difference between the nonlinear and linear wave heights in Fig. 9 is half a meter, when the overpressure is five atmospheres.

To enumerate the consequences of the response predicted by (5) and (6), we observe:

The response wave height $h_{r}$ is linearly proportional to the depth $h_{0}$. A pressure wave over deep ocean has a stronger effect than a pressure wave over a shallower continental shelf.

This explains why locating the blast in the Long Beach case closer to shore had less of an effect. If the distance offshore of the air blast from Long Beach is halved to 15 kilometers, the ocean is only 90 meter deep, resulting in approximately 1/10 the impact response. On the other hand, there is almost no difference in the decay rate of the shallow water waves before they reach shore. 2. 2.

If $s>c_{w}$, then $p_{e}$ and $h_{r}$ have the same sign. The response height is positive in regions of positive overpressure. This contradicts an intuition that positive overpressure would depress the water surface. This response is similar to the case of a forced oscillator in vibrational analysis. Consider for example $\ddot{x}=-x+A\cos(\omega t)$. The steady solution is

[TABLE]

For $\omega>1$, the response $x(t)$ has the opposite sign from the forcing $A\cos(\omega t)$. For pressure forcings with speeds slower than the water speed, the water response would be a depression, with $h_{r}$ negative.

This response is clearly seen in the all the simulations. The wave that travels at the speed of the blast wave is an elevation. However, in the Long Beach results, we can see that the response wave is not uniformly circular when the depth of the water changes rapidly. Note that since the speed of sound in air is 343 m/sec, the water would have to be more than 12 kilometers deep for the gravity wave speed to exceed the speed of the pressure forcing. Hence in all cases on earth we expect an elevation of the sea surface beneath the pressure wave. 3. 3.

The response is particularly strong when the forcing speed $s$ is close to the gravity wave speed $c_{w}\approx 200$ m/sec., (for $h_{0}=4$ km). In this case we have a Proudman resonance [20, 15]. This is the regime for meteo-tsunamis, in basins whose depth leads to gravity wave speeds that match the squall speeds. These speeds are much slower than the speed of sound in air.

3.2 Shallow Water Model Computations

This subsection illustrates the behavior of the analytic model described above with numerical simulations. We solve the equation set (1) again using GeoClaw, with $u=h=p_{e}=0$ for $t\leq 0$, and

[TABLE]

with $p_{ambient}$ = 1 atm. for $t>0$.

These initial conditions use an impulsive start for the air blast pressure wave at time $t=0.0$, so it is not a traveling wave. This generates gravity waves, one moving left and the other moving right, at speeds $c=\pm\sqrt{gh_{0}}$, in addition to the forced water wave traveling at speed $s$. The exact solution to the linearized shallow water equations can be found by matching the conditions at time 0. It is a combination of the forced wave solution to the inhomogeneous equation, presented above in equation (6), plus the solution to the homogeneous equation. Writing the full solution in terms of the response wave $h_{r}$ gives:

[TABLE]

consisting of a left-going and right-going gravity wave traveling with speed $c$, along with what we have been calling the response wave.

Eq.  (7) shows that the left-going tsunami wave will have a smaller amplitude in absolute value for $s/c>1$ than the right going wave. Also, the latter will be a depression, since it has amplitude $-0.5\cdot(1+s/c))$.

Numerical results illustrating this are shown in Fig. 10. The figure contains two curves for each time. Solid lines show the water wave heights, and dashed lines show the air overpressure profiles. In Fig. 10 left, the speed $s=.350$ km/sec, somewhat larger than the speed of sound in air. For this case, since $s>c_{w}$, (for $h_{0}$ = 4 km, $c_{w}=\sqrt{gh_{0}}\simeq$ 0.198 km/sec), the forced wave height is positive since the overpressure is. Note that the gravity wave at the same point in time trails the pressure wave. The right-moving gravity wave is a depression, the left moving wave is a smaller elevation. Since this calculation is in one space dimension, the waves do not decay. In the two dimensional shallow water equations, the gravity wave decays with the square root of distance. Also, the pressure blast wave, and therefore the leading water response would both decay too.

By contrast, Fig. 10 right shows the water’s response for an overpressure moving at $0.120$ km/sec, slower than the gravity wave ($s<c_{w}$). The tsunami waves travel at the same speed in both computations, but they have different amplitudes and signs. The tsunami wave is the opposite sign as the wave due to the pressure. This is consistent with conservation of mass.

To give a more complete picture, two more experiments with $s>c_{w}$ but different forcings are shown in Figure 11. The figures on the left use a Gaussian pressure forcing but their magnitude is ramped up for the first 100 seconds. This results in quite a different-looking gravity wave. The right figure uses a typical Friedlander blast profile described in section 2 for the overpressure, but keeping the amplitude constant at 1 atm. It looks similar to the Gaussian example above.

4 Linearized Euler model

In this section we analyze a more complete model of the ocean’s response to an airburst, to uncover possible shortcomings of the shallow water model of Section 3. We model the water using the Euler equations of a compressible fluid, which will bring in the effects of compressibility and dispersion. Another possibility would be to use one of the forms of the Boussinesq equations, but that also assumes incompressible flow, and would be more difficult to analyze. (See however a nice comparison of SWE and Serre-Green-Nagdhi Boussinesq results in [18].) We continue to neglect Coriolis forces, viscosity, friction, the Earth’s curvature, etc. We linearize the Euler equations and the boundary conditions, since Fig. 9 suggests that linear approximations are reasonably accurate for these parameters.

4.1 Derivation and Analysis

Our starting point for this section is the linearized Euler equations with linearized boundary conditions. A derivation is given in the Appendix. An explicit solution is not possible, and the results will depend instead on wave number. We will use wave number $k=\frac{2\pi}{L}$, where the length scale $L$ for the atmospheric pressure wave is on the order of $10-20$ kilometers. This is very short relative to earthquake-generated tsunamis, which can have length scales on the order of 100 kilometers or more. As before, those not interested in the analysis can skip to the end of the section for a summary of the main points.

The linearized Euler equations and boundary conditions that we use for analysis are

[TABLE]

where $c_{a}$ is the speed of sound in water. (We use $c_{a}$ for acoustic to distinguish it from the gravity wave speed $c_{w}=\sqrt{gh}$). Here, $\widetilde{\rho}$ is a small perturbation of $\rho$ (and the same for the other variables), except for

[TABLE]

where $h_{r}$ is again the water’s disturbance height for consistency with the previous section. The boundary conditions are:

[TABLE]

As in section 3, we will assume the atmospheric pressure forcing has the form $p_{e}(x-st)$, and look for solutions of the same form, functions of $m=x-st$ and $z$. The system (8) becomes

[TABLE]

The boundary conditions become

[TABLE]

This system now includes the effects of dispersion and water compressibility.

These equations cannot be solved in closed form for general $p_{e}$. Therefore, we study the response using Fourier analysis. We will take a Fourier mode of the overpressure

[TABLE]

with amplitude $A_{k}$ and compute the response as a function of $m$. The responses will have the form

[TABLE]

The hat variables are the Fourier multipliers. The partial differential equations (12a-c) become ordinary differential equations with wave number $k$ as a parameter.

Note that (12b) depends only on derivatives with respect to $m$. Integrating it gives

[TABLE]

The constant of integration is zero for each $z$ since as $m\rightarrow\infty$ we know $\widetilde{u}=0$ and $\widetilde{\rho}=0$. This gives an expression for $\widetilde{\rho}$ in terms of $\widetilde{u}$,

[TABLE]

which we can use in (12a) and (12c). After substituting for $\widetilde{\rho}$ and dividing by $\rho_{w}$, the remaining system of two equations is

[TABLE]

Substituting the Fourier modes (15c-d) into (17), and differentiating $\widetilde{u}$ and $\widetilde{w}$ with respect to $m$ gives an ordinary differential equation in $z$ for the velocities,

[TABLE]

The general solution to this 2-by-2 system is the linear combination

[TABLE]

where $\mu_{\pm}$ and $\bf v_{\pm}$ are the eigenvalues and eigenvectors of the matrix in (18), and the scalar coefficients $a_{\pm}$ are chosen to satisfy the boundary conditions. The eigenvalues are

[TABLE]

The eigenvectors (chosen to make the algebra easier so they are not normalized) are

[TABLE]

The boundary condition at $z=0$ is (13a). To apply it, note that $\widehat{w}$ corresponds to the second component of the eigenvectors ${\bf v}_{\pm}$. We find that $a_{+}=-a_{-}$. Henceforth we call this coefficient simply $a$.

Next, we substitute the Fourier modes (15c-d) into the remaining boundary conditions (13b) and (13c). We use the pressure forcing equation (14) in the form $\widehat{p}_{e}=A_{k}$. The result is

[TABLE]

Using equation (16) to substitute $\widehat{\rho}=\frac{s\rho_{w}}{c_{a}^{2}}\widehat{u}$ in (22b) gives an expression for $\widehat{h}_{r}$

[TABLE]

This can be used to replace $\widehat{h}_{r}$ in (22a) to get

[TABLE]

The final steps are using the form of the solution (19) in (24) to solve for the coefficient $a$. With this, everything is known, and $\widehat{u}$, $\widehat{w}$ and the response height $\widehat{h}_{r}$ can be evaluated.

Putting it all together we get

[TABLE]

Grouping terms, the final expression to solve for $a$ (using the definition (20) for $\mu_{\pm}$) is given by

[TABLE]

To summarize, given an overpressure amplitude $A_{k}$ with wavelength $k$, equation (26) gives the scalar coefficient $a$ in the velocity equations, then we solve for $\widehat{u}$ and $\widehat{w}$ using (19), and use (23) to get the Fourier multiplier for the wave height response.

4.2 Linearized Euler Model Computations

We evaluate these results using the following parameters: an ocean with depth $h_{0}=4\,\mbox{km}$, ocean sound speed $c_{a}=1500$ m/sec, $\rho_{w}=1025\,\mbox{kg/m}^{3}$, and atmospheric overpressure of $A_{k}=1\,\mbox{atm}$ with pressure wave speed $s=350$ m/sec, faster than the gravity wave speed of about $200$ m/sec. The responses are linear in the overpressure amplitude $A_{k}$, so we do not evaluate these curves for any other overpressures.

Fig. 12 (left) shows the surface wave height $\widehat{h}(k)$ as a function of length scale $L$, and (right) the amplitude of the surface velocities $\widehat{u}(h_{0},k)$ and $\widehat{w}(h_{0},k)$ are shown. There are two curves in each plot: one uses the physical acoustic water wave speed of $c_{a}=1500$, and the other uses a very large non-physical acoustic speed in the water of $c_{a}\times 10^{8}$. The latter corresponds to the intermediate model of finite depth but incompressible water. This should, and does, asymptote in the long wave ($k\rightarrow 0$) limit to the result of the shallow water equations. The difference between the blue and green curves shows approximately a 20% reduction in the amplitude of the longer length scales due to compressibility (but that this is not the amplitude of the total wave response yet). Note also that the $u$ velocity asymptotes to the shallow water limit, and the $w$ velocity approaches zero. The velocity curves show less of an effect due to compressibility.

For atmospheric forcing from asteroids with air bursts, the length scales of interest are closer to the short end, perhaps 10 or 20 kilometers. In this regime, the compressibility effects are around 10% or less. But at these wavelengths, dispersive effects reduce the response predicted by shallow water theory by nearly half!

This becomes more clear by comparing the forced wave response to a Gaussian pressure pulse instead of using just a single frequency. We use the pressure pulse

[TABLE]

take the Fourier transform, multiply by the Fourier multipliers shown in Fig. 12, and transform back. Figure 13 shows the results for two different water depths $h_{0}$: 4km and 1km. The blue curve uses the water wave speed $c_{a}$=1500 m/sec, and the red curve uses the limiting $c_{a}$. Compressibility changes the height by less than 10% in both figures. However, in the deeper water, the shallow water response is almost 70% larger, and has a narrower width since there is no dispersion. In the right figure, the water is shallower, and the linearized Euler results are closer to the shallow water results.

In Fig. 14 we fix the horizontal length scale at 15 km and instead vary the speed of the pressure wave $s$. This figure again uses $h_{0}=4000$ meters, and $c_{a}=1500$ m/sec. Three curves are shown: the linearized Euler, and the nonlinear and linearized shallow water responses. There is much more difference in this set of curves, particularly around the regions where resonance occurs. Here too we see that the wave height response to the linearized Euler forcing is negative for pressure forcing speeds $s\lesssim 150$ and again unintuitively, positive for larger $s$. There is also a section of the red curve that is missing, corrresponding to the regions where there is no smooth solution. Note also that the overpressure speed where the resonance occurs is significantly slower for the linearized Euler than for the SWE.

5 Conclusions

We have presented several numerical simulations using the shallow water equations over real bathymetry that demonstrate the ocean’s response to a 250MT air-burst. There is no significant wave response from the ocean, in either the forced wave or the gravity waves after a short distance. Our calculations show that the amplitude of the pressure wave response decreases much more rapidly than the gravity waves do. The blast had to be very close to shore to get a sizeable response. Thus the more serious danger from an air burst is not from the tsunami, but from the local effects of the blast wave itself.

Several unexpected features found in the simulations were explained using a one-dimensional model problem with a traveling wave solution for the SWE. One of the main results, that the wave response height is proportional to the depth of the water, explains why putting the blast on a continental shelf close to shore did not generate more inundation than putting it further away in deeper water.

We also looked at the water’s response to an air burst using the linearized Euler equations. In this case the traveling wave model problem shows that the amplitudes of the important wave numbers in the ocean’s response are greatly decreased. We do not yet know what this means for the gravity wave response. In addition, we expect the character of the water’s response to be different, since dispersive waves will generate a wave train characterized by multiple peaks and troughs. The effect of this on land, and whether it causes inundation when the SWE response does not, is something we plan to investigate in the future.

6 Appendix

In this appendix we start with the nonlinear Euler equations for a compressible inviscid fluid with nonlinear boundary conditions at the interface between ocean and air. The static unforced solution to these equations is determined by hydrostatic balance. The hydrostatic pressure is $p_{0}$, and the hydrostatic density is $\rho_{0}$. Since the static density variation is small (under 2%), we will end up neglecting it and proceed to linearize the equations, deriving eq. (8)–(11) in section 4.

This time there are two spatial coordinates, a horizontal coordinate $x$, and a vertical coordinate $z$. The (flat) bottom is $z=0$. The moving top surface is $z=h(x,t)$. The horizontal and vertical velocity components are $u$ and $w$ respectively, and the water density is denoted by $\rho$. The Euler equations are

[TABLE]

There is a “no flow” boundary condition at the bottom boundary,

[TABLE]

The kinematic condition at the top boundary [24] states that a particle that moves with the surface velocity stays on the surface,

[TABLE]

The dynamic boundary condition at the top is continuity of pressure,

[TABLE]

The left side of (30) is pressure in the water evaluated at the top boundary. The right side is the atmosphere’s ambient pressure, which is the sum of the static background atmospheric pressure $p_{\text{a}tm}$ and the dynamic blast wave overpressure $p_{e}(x,t)$.

For static solutions ($p_{e}=0$, $u=w=0$), the nonlinear equations (27) reduce to the hydrostatic balance condition

[TABLE]

Let the water density $\rho_{w}$ be the density at the water surface. If the density differences are small (as they turn out to be), we may use a linear approximation to the equation of state,

[TABLE]

where $c_{a}$ is the acoustic sound speed in water at density $\rho_{w}$, $c_{a}^{2}=\frac{dp}{d\rho}(\rho_{w})\;.$ The behavior of $\rho_{0}(z)$ is found by substituting this into (31): $\frac{dp_{0}}{dz}=c_{a}^{2}\frac{d\rho_{0}}{dz}=-g\rho_{0}\;.$ Therefore, for any two heights $z_{1}$ and $z_{2}$, we have

[TABLE]

If $z_{2}-z_{1}=4\,\mbox{km}$, and $c_{a}=1500\,\frac{\mbox{m}}{\mbox{sec}}$, then $\frac{g}{c_{a}^{2}}(z_{2}-z_{1})<.02\;.$ Therefore, the density varies by less than about $2\%$ between the water surface and bottom.

We denote small disturbance quantities with a tilde, except for the wave height response $h_{r}$, which we use for continuity with the previous sections. For example, the water density is $\rho_{0}(z)+\widetilde{\rho}(x,z,t)$. These disturbances are driven by the atmospheric overpressure $p_{e}(x,t)$. We substitute the expressions $\rho=\rho_{0}+\widetilde{\rho}$, $u=\widetilde{u}$, $w=\widetilde{w}$ (since the velocities are linearized around zero), and $p=p_{0}+c_{a}^{2}\widetilde{\rho}$ into the Euler equations (27) and calculate up to linear terms in the disturbance variables. Using the hydrostatic balance condition (31), this gives

[TABLE]

Finally, we replace the (slightly) variable $\rho_{0}(z)$ with the constant $\rho_{w}$. The resulting equations, which we use for analysis are

[TABLE]

The bottom boundary condition (28) is already linear. For the top boundary conditions, we express the water height as the sum of the background height $h_{0}$ and the disturbance height $h_{r}$:

[TABLE]

To leading order in $h_{r}$, $\widetilde{u}$ and $\widetilde{w}$, the linear approximation to the kinematic boundary condition (29) is

[TABLE]

For the dynamic boundary condition (30), which was $p(x,h(x,t),t)=p_{atm}+p_{e}$, we use the Taylor expansion and the perturbation approximation

[TABLE]

For the undisturbed quantities, the pressure at the top is $p(h_{0})=p_{\mbox{\scriptsize\em atm}}$. The hydrostatic balance relation (31) in the water (applied at the top) is $p_{0,z}(h_{0})=-g\rho_{w}$. Making these substitutions gives

[TABLE]

giving the result

[TABLE]

Summarizing, the linearized Euler equations are (32), with linearized boundary conditions (28), (33) and (35).

Acknowledgments

We are particularly grateful to Mike Aftosmis, Oliver Bühler, and Randy LeVeque for more in-depth discussions. It is a pleasure to thank our colleagues at the Courant Institute for several lively discussions. We thank Michael Aftosmis for providing the blast wave ground footprint model. This effort was partially supported through a subcontract with Science and Technology Corporation (STC) under NASA Contract NNA16BD60C.