A vertical eigenfunction equation is solved to examine the partial reflection and partial transmission of tsunami-generated gravity waves propagating through a height-dependent background atmosphere from the ocean surface into the lower thermosphere. There are multiple turning points for each vertical eigenfunction (at least eight in one example), yet the wave transmission into the thermosphere is significant. Examples are given for gravity wave propagation through an idealized wind jet centered near the mesopause and through a realistic vertical profile of wind and temperature relevant to the tsunami generated by the Sumatra earthquake on 26 December 2004.
A tsunami is a long ocean surface wave caused most often by an undersea or coastal earthquake. As a tsunami travels across the ocean, its movement generates gravity waves in the atmosphere. Numerical simulations indicate that these tsunami-generated atmospheric gravity waves can propagate rapidly upward and within a few hours reach peak amplitudes at altitudes of 250 km or higher (e.g., Hickey et al. 2009; Occhipinti et al. 2011). An improved understanding of tsunami-generated gravity waves could conceivably have practical applications for early tsunami detection and warning through near real-time atmospheric monitoring (Artru et al. 2005; Lovett 2010; Liu et al. 2006).
Like the tsunami itself, the tsunami-generated gravity waves have long horizontal wavelengths (hundreds of kilometers), fast horizontal phase speeds (~200 m s−1), and small amplitudes at the ocean surface (tens of centimeters). These properties allow the gravity waves to propagate upward through the lower and middle atmosphere without much geometrical spreading or critical-layer interaction, and without significant dissipation or wave breaking, even as the gravity wave amplitudes grow exponentially with height in response to the decreasing air density. One process that can affect tsunami-generated gravity waves in the lower and middle atmosphere is partial reflection, as noted in the early work of Peltier and Hines (1976).
Here we examine the partial reflection of tsunami-generated gravity waves. The altitude range of our model extends from the ocean surface into the lower thermosphere, to a maximum altitude of 180 km. We solve a vertical eigenfunction equation for the gravity waves and find that each vertical eigenfunction has multiple turning points—at least eight in one example. The gravity waves are thus converted between propagating and evanescent waves several times during upward propagation. The amplitude of the transmitted waves is significant and comparable to the amplitude of the reflected waves. The problem is related to wave tunneling through a series of thin evanescent layers (Fröman and Fröman 2002; Brown and Sutherland 2007).
Hickey et al. (2009) and Occhipinti et al. (2006, 2008, 2011) also calculate solutions related to the vertical eigenfunctions for tsunami-generated gravity waves. Their models extend to higher altitudes than ours and include processes relevant to the ionosphere that we omit, such as wave coupling to the plasma. Hickey et al. (2009) solve the Navier–Stokes equations for each Fourier component of the gravity wave spectrum. Occhipinti et al. (2006, 2008, 2011) solve a set of inviscid equations for each Fourier component. In the latter case, the model does not include downgoing reflected waves or background winds.
We supplement these studies with a more detailed analysis of partial reflection and partial transmission for altitudes up to the lower thermosphere. The present model is anelastic, inviscid, steady state in the reference frame moving with the tsunami, and two dimensional in a vertical plane with a height-dependent atmospheric background, including winds.
In one example, we use this model to study atmospheric gravity waves generated by the devastating tsunami that resulted from the Sumatra earthquake on 26 December 2004 (Titov et al. 2005). Gravity waves generated by this tsunami were also modeled by Hickey et al. (2009) and Occhipinti et al. (2006). As in these studies, we look at a location on the equator over the Indian Ocean west of Sumatra. Ionospheric measurements taken from as far away as Puerto Rico have identified wavelike disturbances interpreted as atmospheric gravity waves generated by the Sumatra tsunami (Lee et al. 2008).
The model formulation is given in section 2. The numerical technique for the vertical eigenfunction equation is described in section 3. An accompanying ray-tracing integration is introduced in section 4. The results are presented in section 5 for propagation through three different background atmospheres: a uniform background, a localized wind jet centered at 100-km altitude, and a realistic vertical profile of wind and temperature representative of atmospheric conditions at the time and location of the Sumatra earthquake. A summary of our main findings is given in section 6.
The model is linear, two dimensional, and steady state in a coordinate system (x, z) that moves horizontally at the phase speed c of the tsunami. This is a standard formulation for this application (e.g., Hickey et al. 2009; Artru et al. 2005). As in Hickey et al. (2009), we set c = −200 m s−1, corresponding to an ocean depth of about 4 km. For simplicity, we assume that the tsunami and gravity waves propagate in the east–west plane, so that x is a zonal coordinate and negative x is westward. The undisturbed ocean surface is at z = 0.
In the tsunami reference frame, the background wind is
where U′(z) is the background wind velocity in the reference frame fixed to Earth.
The tsunami displaces the ocean surface vertically by
where a is a constant, Ai is the Airy function, and x′ is x measured in units of 100 km. This function was used by Peltier and Hines (1976) and Hickey et al. (2009). The Airy function can be associated with a caustic near the front of the tsunami (Berry 2005).
The vertical velocity W(x) at the ocean surface produced by flow over the tsunami is expressed by the Fourier integral
Here U0 = U(z = 0), is the Fourier transform of h(x), and k is the horizontal wavenumber.
We mention a few characteristics of the tsunami, some of which can be found in Hickey et al. (2009). The wavenumber of maximum spectral amplitude is K0 = 1.55 × 10−5 m−1, with a wavelength 2π/K0 of about 404 km. For our choice of c = −200 m s−1, and for the wavenumber K0, the intrinsic wave period of the tsunami (in the reference frame fixed to Earth) is 33 min. For a maximum h of 0.5 m, as used in Hickey et al. (2009), and a ground wind speed of U′ = 10 m s−1, the maximum magnitude of the surface vertical velocity W(x) is W0 ≃ 1.2 × 10−3 m s−1.
The gravity waves are stationary in the tsunami reference frame. The vertical velocity of the gravity waves is (ρ0/ρ)1/2w(x, z), where ρ(z) is the background density with value ρ0 at z = 0. The density factor (ρ0/ρ)1/2 accounts for non-Boussinesq effects in the anelastic approximation [Smith 1979, his (2.22)] and increases exponentially with altitude. The plots in this paper show w(x, z) without the density factor in order to reveal the solution more clearly at all altitudes.
The main results of this paper are obtained by solving a vertical eigenfunction equation for , which is the horizontal Fourier transform of w(x, z). Spatial solutions are then calculated by the inverse Fourier transform
The vertical eigenfunction equation is
[Inverarity 2003, his Eq. (6.2)]. Here N(z) is the buoyancy frequency, S(z) = 1/H(z) = −(dρ/dz)/ρ is the reciprocal of the density scale height H(z), and Uzz is the second derivative of U with respect to z. Turning points occur where m2 = 0. The waves are propagating where m2 > 0 and evanescent where m2 < 0. The upper and lower boundary conditions for in (6) are given in the next section.
3. Numerical solution of the vertical eigenfunction equation
The numerical method for the vertical eigenfunction equation is standard and involves integrating (6) as two coupled first-order ordinary differential equations for and . The integration is performed for each Fourier component that is initially propagating (i.e., m2 > 0 at z = 0). The integration starts at the upper boundary at z = 180 km and proceeds down to z = 0 where the solution is scaled to satisfy the lower boundary condition (see below). At the upper boundary the solution is chosen to represent either an upgoing wave if m is real or a wave that decays exponentially with altitude if m is imaginary. Broutman et al. (2009) used this same method for mountain waves (of much shorter horizontal wavelength than here) and compared the numerical method with ray and uniformly valid approximations.
It is useful for interpretation to express the numerical solution wherever possible as the sum of slowly varying upgoing and downgoing waves. We write
where m = θz and the complex amplitudes A, B depend on z. When m is real, we choose the positive root for m in (7), so the terms with A and B represent upgoing and downgoing wave groups, respectively. Rearranging (8) and (9) gives the upward- and downward-propagating components as
[Fröman and Fröman 2002, their (2.2.1)]. This expression for ϵ is equivalent to that used in Einaudi and Hines (1970) but is expressed in terms of m2 (rather than m), which is real for all z, including evanescent regions.
We compute WKB solutions at altitudes where ϵ < 1. Broutman et al. (2009, p. 486) used and tested this criterion for an application with mountain waves and turning points. A smaller range of ϵ < 0.5 slightly decreases the extent of the region where WKB solutions are calculated (as plotted here in Figs. 5 and 10).
The lower boundary condition is
where A corresponds to the upgoing wave component in (8) and (10). This condition treats any downgoing gravity waves that reach the lower boundary as outgoing waves: there is no reflection back upward from the ocean surface. Occhipinti et al. (2008, p. 3) also use a lower boundary condition that involves only the upgoing wave component and a slowly varying approximation at the lower boundary.
For comparison, we also computed solutions with the lower boundary condition
where is given in (4). This treats the ocean as a reflecting surface for downward-propagating gravity waves. The gravity waves reflected upward from the ocean surface interact with the background and the turning points in the same way as the originally upward-propagating waves, and can result in constructive or destructive interference. Since our model is steady state, (14) corresponds to an infinite number of gravity wave reflections between the ground and the turning points.
We found in our examples that when (14) is used instead of (13), our main conclusions are the same: there is an approximate standing wave pattern below the turning points indicating significant partial reflection, and there is significant partial transmission through the turning points. The results presented below are obtained using (13) [except for the example of Fig. 7, which is calculated using (14) and is included for comparison].
4. Ray tracing
Some additional insight is provided by calculating ray paths for the gravity waves. The ray equation for the ray path x(z) is
where m is given by (7). In this context, x and k refer to values along the ray path.
The initial condition k(x, z = 0) for the ray tracing is obtained from the tsunami profile at the ocean surface. The tsunami is treated as a slowly varying wave train with a locally defined horizontal wavenumber K(x). The wavenumber K(x) is calculated from the vertical velocity W(x) in (2) and (3) using a Hilbert transform (Hahn 1996). The Hilbert transform gives an imaginary part for W(x), which can be added to the real-valued W(x) to form a complex function. The phase ϕ(x) of that complex function is then differentiated (using a finite-difference approximation) to obtain the local wavenumber K(x) = ϕx. The initial condition for the ray tracing is k(x, z = 0) = K(x).
Estimating the local wavenumber in this way is a well-known signal-processing application of the Hilbert transform. The Hilbert transform is numerically implemented here by the standard method involving a one-sided Fourier transform. Further information on the one-sided Fourier transform and the phase calculation from the Hilbert transform can be found in Hahn (1996, his sections 1.7 and 1.13).
Figure 1 shows W(x) in the top panel and the associated wavenumber K(x) in the bottom panel. For the ray tracing, we use the K values between the plotted circles in the bottom panel of Fig. 1. The circles have x values of −180 and 750 km and K/K0 values of −0.46 and −1.6. Over this range |Kx/K2| < 0.2, which is small compared to unity and thus an indication of a slowly varying wave profile and a well-defined local wavenumber K(x). The position at which K(x) = −K0 is x ≃ 360 km.
As the rays propagate upward, they encounter evanescent regions where m is imaginary and the ray path is not defined. In an evanescent region, we simply continue the ray straight upward (without horizontal displacement) until the waves are just above the evanescent region and become propagating again, with m real. This seems to be adequate for the thin evanescent regions in the results below. We ignore downgoing rays reflected from turning points.
We consider three model atmospheres: a uniform background (case 1), an idealized wind jet (case 2), and a background obtained from empirical atmospheric models for the Sumatra tsunami of 26 December 2004 (case 3).
The Fourier integral for w(x, z) in (5) is approximated by an FFT with 1024 points and a grid spacing of 15 km in x. In the following we will sometimes refer to the local complex magnitude of the vertical velocity. We denote this by |w|, a function of x and z that is obtained from w(x, z) using a numerical Hilbert transform (Hahn 1996).
Since m2 is real in our calculations, it follows that
is constant with height (Heading 1962, his section 4.4), where the asterisk indicates a complex conjugate and Im indicates the imaginary part. The quantity F is proportional to the wave momentum flux. In the examples below, the numerical solution achieved a constancy of F to at least three significant figures for all Fourier components.
a. Case 1: Uniform background
This is a reference case mainly for comparison with later examples. The background quantities in (7) are uniform, with U′ = 10 m s−1, N = 0.02 s−1, and H = 6 km. Figure 2 shows the vertical velocity w(x, z). Ray paths computed from (15) are superimposed on the plot with the thick line indicating the ray with k = −K0, the wavenumber of maximum |W(k)|.
Geometrical spreading is very weak in this example, as illustrated by Fig. 3. The local magnitude of the vertical velocity |w| is plotted as a function of x at two heights: z = 0 (dashed line) and z = 150 km. The latter is obtained from the Fourier solution [(5); solid line] and the ray-tracing solution (circles), which are in close agreement. The peak |w| near x = 300 km decreases with altitude owing to geometrical spreading, but only slightly, by about 3% from z = 0 to z = 150 km.
b. Case 2: Idealized wind jet
We next consider a simple case where most of the initially propagating vertical eigenfunctions have two turning points. As in case 1, we set N = 0.02 s−1 and H = 6 km. We add a wind jet of the form
with U0 = 10 m s−1, U1 = 60 m s−1, zm = 100 km, and L = 10 km. The form (17) is an idealization of the winds associated with migrating tides above the mesopause region [as in Hickey et al. (2009, their Fig. 3) and Fig. 8].
The top panel of Fig. 4 shows computed from (7). This plotted band of wavenumbers contains the k values that yield propagating waves at z = 0. The plotted curves correspond to every fifth initially propagating value of k on the Fourier grid used to calculate the wave solutions.
The bottom panel of Fig. 4 shows the turning-point heights as a function of k. Also shown for comparison is the shape of the tsunami spectrum |W(k)| (bottom of the plot) and the shape of the wind jet (dashed line at the left). Over most of the spectrum there are two turning points, which are located near 100- and 120-km altitudes. For k/K0 < −2.2 the number of turning points increases to four. For k/K0 < −2.95 the waves are evanescent at z = 0 and are not included in this calculation.
Figure 5 shows w(x, z)/W0 in the top panel. There is substantial partial reflection from the wind jet, revealed by the almost vertical orientation of the wave phases below the wind jet, representing a standing wave pattern. The upgoing and downgoing components of the wavefield, calculated from (8)–(11), are plotted respectively in the middle and bottom panels. These are calculated only in regions where ϵ < 1 in (12), leaving gaps in the plotted solution near the wind jet.
where A, B are given in (8)–(11). The subscripts 0, 1 refer, respectively, to values at z = 0 and at the top of the domain z = 180 km. Constancy of the momentum flux [see (16)] gives R2(k) + T2(k) = 1. For the present example m1 = m0.
At k = −K0, the reflection and transmission coefficients of Fig. 6 are both about 0.71. These values compare well with the spatial solution. At the top of the domain, there are only upgoing waves, and the maximum spatial solution |w|/W0 is about 0.67. At z = 0, the maximum |w|/W0 for the downgoing component is about 0.66.
It is useful to compare these results with the predictions of wave tunneling theory, with two turning points surrounding a single evanescent region. In the simplest case m2 is a quadratic function of z in the turning-point region and is vertically symmetric about the turning-point height (Fröman and Fröman 2002, their section 3.49).
The reflection and transmission coefficients for the two turning-point case can be derived theoretically and depend on the quantity
where z = a, b are the heights of the two turning points and mi is the imaginary part of m, assumed positive for this discussion.
The two turning-point theory gives the square of the reflection and transmission coefficients as
[Fröman and Fröman 2002, their (3.49.4)–(3.49.7)]. When the two turning points are well separated in the sense that μ is large, the transmission coefficient is exponentially small in μ, and the reflection coefficient is close to unity.
For a thin evanescent layer in the limit μ → 0, the reflection and transmission coefficients are equal, with value . The fact that this is so close to R(k) and T(k) for k = −K0 in this example is apparently fortuitous. Still, a basic result from two turning-point theory is relevant here: for closely spaced turning points the reflected and transmitted waves have similar amplitudes. In the present example, the Fourier components with two turning points (k > −2.2K0) have μ values that range from 0.74 to 1.05. This gives values of T2 that range from 0.50 to 0.57, somewhat smaller than the range for T of 0.50–0.74 in Fig. 6 for these wavenumbers.
Figure 7 shows the solution w(x, z)/W0 for the reflecting lower boundary condition (14). It has the same features as in the top panel of Fig. 5, which was computed with the nonreflecting lower boundary condition [(13)]. There is again a standing wave pattern below the turning points and transmitted waves above the turning points. The main difference is the overall amplitude of the solution. The peak magnitude for |w|/W0 is about 1.6 in the case of Fig. 5 and about 1.1 in the case of Fig. 7. This implies some destructive interference caused by the gravity waves reflected from the ocean surface.
As another test of the parameters for case 2, we increased the maximum wind jet speed to 100 m s−1, setting U1 = 90 m s−1 in (17), with all other parameters the same as before and with the lower boundary condition [(13)]. This increases the reflection coefficient, but the transmission coefficient is still relatively high. For k > −1.5K0 the range of T is 0.49–0.58. At this wind speed, the wind jet half-width in (17) has to be increased to L = 30 km to reduce the transmission coefficients to less than 0.1 for all k.
c. Case 3: The Sumatra tsunami
We now examine gravity waves from the tsunami following the Sumatra earthquake on 26 December 2004. Most of the tsunami energy propagated westward from Sumatra into the Indian Ocean (Titov et al. 2005). As in Hickey et al. (2009) and Occhipinti et al. (2006), we obtain background atmospheric profiles for 0300 UTC on the day of the earthquake, and for 0°, 85°E—a location in the Indian Ocean roughly 10° west of the earthquake epicenter.
The vertical profiles of density ρ(z) and wind velocity U′ are obtained, respectively, from the empirical models Mass Spectrometer and Incoherent Scatter Radar (MSIS) and Horizontal Wind Model (HWM), as described in Picone et al. (2002) and Drob et al. (2008). There are several versions of these models, dating back to the early 1990s. We use the latest versions NRLMSISE-00 and HWM07 (Drob et al. 2008). The wind model HWM07 is a major upgrade from earlier versions and has a completely revised specification of the migrating tides in the mesosphere and lower thermosphere. Therefore our wind profile U′(z) is different from that used by Hickey et al. (2009), which was derived from an earlier version known as HWM93.
The vertical profiles for U′, N, and H are shown in Fig. 8. The output from MSIS leads to large (and likely spurious) oscillations for the calculated N(z) in the lower troposphere; see also the oscillations in buoyancy period in Fig. 2 of Hickey et al. (2009). For our calculation, these oscillations were replaced by a representative tropospheric value of N = 0.01 s−1 for altitudes up to z = 15 km. The sudden change in the slope of N(z) near the equatorial tropopause at 15-km altitude generates some partial reflection, but the transmission is large, as indicated by the results below.
Figure 9 shows for the range of Fourier components that are propagating at z = 0. For k > −2.5K0 there are eight turning points per Fourier component and four intervening evanescent regions, all occurring at altitudes below 110 km. For these wavenumbers the transmission coefficient T(k) is greater than 0.55, as shown in the bottom panel of Fig. 9. For k < −2.5K0 there are as many as 15 turning points per Fourier component.
Figure 10 shows w(x, z)/W0 in the top panel. The standing wave pattern below ~100 km indicates significant partial reflection. There is also significant partial transmission. Above 100-km altitude the waves are well represented by the solution for upgoing waves, with phase lines tilting westward with increasing altitude. The upgoing (middle panel) and downgoing (bottom panel) components are computed from (10) and (11).
Figure 11 shows the vertical eigenfunctions for k = −K0 (thick line) and surrounding wavenumbers ranging from −0.74K0 to −1.26K0. The horizontal dashed lines show turning-point heights for k = −K0.
The vertical eigenfunctions pass through the evanescent regions without much decrease in amplitude. For example, note the curve farthest to the left near z = 80-km altitude. Above the turning point near z = 77 km, the curve seems to start to decay with height as an evanescent wave (m2 < 0), but this is shortly interrupted by the presence of another turning point near z = 83 km, above which the wave resumes vertical propagation (m2 > 0). Above 100-km altitude, there is a sharp reduction in the amplitude of some of the eigenfunctions owing to the absence of downgoing waves above the highest turning point, as shown in the bottom panel of Fig. 10.
Figure 12 shows the individual terms in m2 of (7). All terms are normalized by dividing by . The largest positive term is N2/U2. The largest negative term is, at most altitudes, −S2/4. The term with the sharpest vertical fluctuations is −Uzz/U, which has large local negative values near 50- and 100-km altitudes. The only term in m2 missing from the plot is k2, which is relatively small for this case. The most energetic Fourier components have of order unity.
The term N2/U2 in (7) is generally dominant in gravity wave applications. For tsunami-generated gravity waves this term is relatively small because U = U′ − c includes the large phase speed c of the tsunami. For realistic conditions the other terms of m2 can locally dominate N2/U2 and cause m2 to change sign repeatedly with altitude, explaining the abundance of turning points. While the exact altitudes of the turning points may be sensitive to some of the terms in m2, particularly −Uzz/U, the occurrence of multiple turning points is not sensitive to these details based on the experiments we have performed so far.
Hickey et al. (2009) modeled gravity waves generated by the Sumatra tsunami and found that the background winds can strongly limit the ability of the gravity waves to propagate into regions of the thermosphere above about 200-km altitude. This wind filtering effect was explained by the reflection of the gravity waves. In the case of their Fig. 6, it appears that most of the wind filtering occurs above about 180-km altitude, and that by 250-km altitude the wind filtering causes more than an order of magnitude reduction in the gravity wave amplitudes.
Our model does not cover these high altitudes, but as in Hickey et al. (2009) we did not find strong wind filtering below the thermosphere. We computed gravity wave solutions using background profiles from the latest and earlier versions of MSIS and HWM, for different times of the day during the Sumatra tsunami, and we also tested the effects of a strong wind jet (100 m s−1 peak speed) just above the mesopause region, as in Fig. 3 of Hickey et al. (2009). In all of these cases, the evanescent layers were sufficiently localized in the vertical to allow significant wave transmission. For altitudes below 150 km, partial reflection accounted for at most a factor of about 2 reduction in the upgoing wave amplitudes.
Above 200 km, the wind filtering as a result of turning points is likely to be stronger because the vertical scales of the background are much longer; see Fig. 3 of Hickey et al. (2009), particularly the east–west component of the wind, which is responsible for most of their wind filtering. The larger scales result in wider evanescent layers above turning points and a potentially greater exponential decay of the wave amplitudes with altitude.
We considered two cases with background atmospheres that vary with height: an idealized wind jet centered near the mesopause and a realistic background for the tsunami following the Sumatra earthquake. In both cases all of the computed vertical eigenfunctions were evanescent for part of the way up to the thermosphere, with alternating layers of propagating (m2 > 0) and evanescent (m2 < 0) waves. In the Sumatra case, there were at least eight turning points for each vertical eigenfunction, with four intervening evanescent layers.
We found comparable magnitudes for the reflection and transmission coefficients, defined in (18) and plotted in Figs. 6 and 9. Two turning-point theory for a thin evanescent layer also gives similar magnitudes for the reflection and transmission coefficients, as noted in section 5. We modeled tsunami-generated gravity waves up to 180-km altitude. At higher altitudes, the vertical scales of the background variations are larger, as in the profiles of buoyancy period and wind in Figs. 2 and 3 of Hickey et al. (2009). This can result in wider evanescent layers and stronger partial reflection.
We computed spatial ray paths for the gravity waves. The ray corresponding to the peak amplitude of the tsunami spectrum (indicated by the thick curves in Figs. 2, 5, and 10) tracked the largest gravity wave amplitudes reasonably well. It may be of interest to continue the ray paths to higher altitudes with a form of the dispersion relation that includes dissipative terms for the thermosphere (Vadas and Fritts 2009; Vadas et al. 2009).
We ignored dissipation, the effects of a compressible atmosphere, and the reflection of downgoing waves from the ocean surface (except for the case of Fig. 7). We also omitted Fourier components that are evanescent at the ocean surface, even though they may become propagating a short distance above the ocean surface. The present results suggest that there would be significant wave transmission through a near-surface evanescent layer caused, for example, by a locally small value of N as might occur in a well-mixed boundary layer.
This research was partially supported by the Office of Naval Research through the NRL 6.1 work unit “The Boundary Paradox.” We thank Prof. M. Hickey for an early discussion.