Testing the hypothesis that it's a lamb wave

The mathematics of lamb waves

In the atmosphere there is a special kind of acoustic wave called a Lamb wave, which is a horizontally propagating sound wave with no vertical velocity propagation (w′w' = 0). They are supported by a hydrostatic equation set. It has been observed that the amplitudes of the oscillations of a Lamb wave vary with height.

Use the equations: (∂∂t+u‾∂∂x)u′+1ρ‾∂p′∂x=0 \left(\frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x} \right)u' + \frac{1}{\overline{\rho}} \frac{\partial p'}{\partial x} = 0 (∂∂t+u‾∂∂x)p′+γp‾∂u′∂x=0 \left(\frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x} \right)p' + \gamma \overline{p} \frac{\partial u'}{\partial x} = 0

∂p′∂z=−ρ′g\frac{\partial p'}{\partial z} = -\rho' g (∂∂t+u‾∂∂x)ρ′+ρ‾∂u′∂x=0\left(\frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x} \right) \rho' + \overline{\rho} \frac{\partial u'}{\partial x} = 0

Where p‾(z)\overline{p}(z), ρ‾(z)\overline{\rho}(z) depend on height, u‾(z,ϕ)\overline{u}(z, \phi) , T‾(z,ϕ)\overline{T}(z, \phi) are taken as parameters and γ≡cp/cv\gamma \equiv c_p/c_v and eliminate ρ′\rho' in the last equation:

−1g(∂∂t+u‾∂∂x)∂p′∂z+ρ‾∂u′∂x=0-\frac{1}{g}\left(\frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x} \right) \frac{\partial p'}{\partial z} + \overline{\rho} \frac{\partial u'}{\partial x} = 0

We multiply through by −g-g to get

(∂∂t+u‾∂∂x)∂p′∂z−gρ‾∂u′∂x=0 \left(\frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x} \right) \frac{\partial p'}{\partial z} -g \overline{\rho} \frac{\partial u'}{\partial x} = 0

Before we go further that if ρ=ρ‾+ρ′\rho = \overline{\rho} + \rho' and p=p‾+p′, p = \overline{p} + p', then we find that ∂p∂z=−ρg\frac{\partial p}{\partial z} = - \rho g which gives us ∂p′∂z+∂p‾∂z=−g(ρ‾+ρ′)\frac{\partial p'}{\partial z} + \frac{\partial \overline{p}}{\partial z} = - g \left(\overline{\rho} + \rho' \right) and using what we know about ρ′\rho' and p′p' we get that ∂p‾∂z=−gρ‾\frac{\partial\overline{p}}{\partial z} = -g\overline{\rho} which should be true just from physical intuition. We note that γ≠0\gamma \neq 0 and p‾≠0\overline{p} \neq 0 and thus

p‾γ(∂∂t+u‾∂∂x)∂p′∂z−gγp‾ρ‾∂u′∂x=0 \overline{p} \gamma \left(\frac{\partial}\partial {t} + \overline{u} \frac{\partial}{\partial x} \right) \frac{\partial p'}{\partial z} - g \gamma \overline{p} \overline{\rho} \frac{\partial u'}{\partial x} = 0 gρ‾(∂∂t+u‾∂∂x)p′+gγρ‾p‾∂u′∂x=0 g \overline{\rho} \left(\frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x} \right)p' + g \gamma \overline{\rho} \overline{p} \frac{\partial u'}{\partial x} = 0

and we add these equations to get

p‾γ(∂∂t+u‾∂∂x)∂p′∂z+gρ‾(∂∂t+u‾∂∂x)p′=0 \overline{p} \gamma \left(\frac{\partial }{\partial t} + \overline{u} \frac{\partial}{\partial x} \right) \frac{\partial p'}{\partial z} + g \overline{\rho} \left(\frac{\partial }{\partial t} + \overline{u} \frac{\partial }{\partial x} \right)p' = 0 (∂∂t+u‾∂∂x)u′+1ρ‾∂p′∂x=0 \left(\frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x} \right)u' + \frac{1}{\overline{\rho}} \frac{\partial p'}{\partial x} = 0

Which allows us to solve for p′p' (note we started with 3 equations and 4 unknowns) and then back-substitute to get other values: assume that p′=p^exp⁡(ik(x−ct)) p' = \hat{p}\exp(ik(x - ct)), p^(z=0,ϕ)=p00\hat{p}(z=0, \phi) = p_{00} and get

p‾γ(−kc+u‾k)(∂p^∂z)Ψ+igρ‾(ku‾−kc)p^Ψ=0 \overline{p}\gamma\left( -kc + \overline{u}k \right)\left(\frac{\partial \hat{p}}{\partial z} \right)\Psi + ig\overline{\rho} \left( k\overline{u} - k c\right)\hat{p}\Psi = 0

  ⟹  ikp‾γ(u‾−c)(∂p^∂z)Ψ+ikgρ‾(u‾−c)p^Ψ=0 \implies ik\overline{p}\gamma\left( \overline{u} -c \right)\left(\frac{\partial \hat{p}}{\partial z} \right)\Psi + ik g\overline{\rho} \left( \overline{u} - c\right)\hat{p}\Psi = 0

  ⟹  p‾γ(∂p^∂z)+gρ‾p^=0 \implies \overline{p}\gamma\left(\frac{\partial \hat{p}}{\partial z} \right) + g\overline{\rho}\hat{p} = 0   ⟹  (∂p^∂z)=−gρ‾p‾γp^ \implies \left(\frac{\partial \hat{p}}{\partial z} \right) = -\frac{g\overline{\rho}}{\overline{p}\gamma}\hat{p}

and we use the ideal gas law to find p‾=ρ‾RdT‾  ⟹  p‾ρ‾=RdT‾ \overline{p} = \overline{\rho} R_d \overline{T} \implies \frac{\overline{p}}{\overline{\rho}} = R_d \overline{T} and therefore

(∂p^∂z)=−gRdγT‾(z)p^\left(\frac{\partial \hat{p}}{\partial z} \right) = - \frac{g}{R_d\gamma \overline{T}(z)}\hat{p}

The righthand size has units Pa/m which agrees with the lefthand side. Clearly this has solution ln⁡(p^(z))−ln⁡(p^(z=0))=−∫z=0zgRdγT‾(z) dz\ln(\hat{p}(z)) - \ln(\hat{p}(z=0)) = - \int_{z=0}^z\frac{g}{R_d \gamma \overline{T}(z)}\, \mathrm{d}z and therefore p^(z)=p00exp⁡(−∫z=0zgRdγT‾(z) dz)\hat{p}(z) = p_{00}\exp\left(- \int_{z=0}^z\frac{g}{R_d \gamma \overline{T}(z)}\, \mathrm{d}z\right)

In order to find the physical solution we take

p′ p' =Re⁡(p00exp⁡(−∫z=0zgRdγT‾(z) dz)exp⁡(ik(x−ct)))= \Re \left(p_{00}\exp\left(- \int_{z=0}^z\frac{g}{R_d \gamma \overline{T}(z)}\, \mathrm{d}z\right) \exp \left(ik(x - ct) \right) \right)
=p00exp⁡(−∫z=0zgRdγT‾(z) dz)Re⁡(exp⁡(ik(x−ct)))= p_{00}\exp\left(- \int_{z=0}^z\frac{g}{R_d \gamma \overline{T}(z)}\, \mathrm{d}z\right) \Re \left( \exp \left(ik(x - ct) \right) \right)
=p00exp⁡(−∫z=0zgRdγT‾(z) dz)Re⁡(cos⁡(k(x−ct))+isin⁡(k(x−ct)))= p_{00}\exp\left(- \int_{z=0}^z\frac{g}{R_d \gamma \overline{T}(z)}\, \mathrm{d}z\right) \Re \left( \cos(k(x - ct)) + i \sin(k(x - ct)) \right)
=p00exp⁡(−∫z=0zgRdγT‾(z) dz)cos⁡(k(x−ct)) = p_{00}\exp\left(- \int_{z=0}^z\frac{g}{R_d \gamma \overline{T}(z)}\, \mathrm{d}z\right) \cos(k(x - ct))
And we use the linearized hydrostatic equation ∂p′∂z=−ρ′g\frac{\partial p'}{\partial z} = -\rho' g to find
ρ′\rho' =−1g∂p′∂z= -\frac{1}{g} \frac{\partial p'}{\partial z}
=−p00RdγT‾(z)exp⁡(−∫z=0zgRdγT‾(z) dz)cos⁡(k(x−ct))= - \frac{p_{00}}{R_d\gamma \overline{T}(z)}\exp\left(- \int_{z=0}^z\frac{g}{R_d \gamma \overline{T}(z)}\, \mathrm{d}z\right) \cos(k(x-ct))

We return to the remaining equation (∂∂t+u‾∂∂x)u′+1ρ‾∂p′∂x=0\left( \frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x}\right)u' + \frac{1}{\overline{\rho}} \frac{\partial p'}{\partial x} = 0 and assume that u′=Re⁡[u^(z)exp⁡(ik(x−ct))]u' = \Re \left[\hat{u}(z) \exp(ik(x-ct)) \right] to find

(−ikc+iu‾k)u^(z)Ψ+1ρ‾(ik)p^Ψ (-ikc + i\overline{u}k)\hat{u}(z)\Psi + \frac{1}{\overline{\rho}} (ik) \hat{p} \Psi   ⟹  (−c+u‾)u^(z)+p^ρ‾ \implies (-c + \overline{u})\hat{u}(z) + \frac{\hat{p}}{\overline{\rho}}   ⟹  u^(z)=−p^ρ‾(u‾−c) \implies \hat{u}(z) = - \frac{\hat{p}}{\overline{\rho}( \overline{u}-c)}   ⟹  u^(z)=−p00ρ‾(u‾(z)−c)exp⁡(−gRdγT‾z)cos⁡(k(x−ct)) \implies \hat{u}(z) = - \frac{p_{00}}{\overline{\rho}( \overline{u}(z)-c)} \exp\left(-\frac{g}{R_d\gamma \overline{T}}z\right) \cos(k(x-ct))

We also return to the equation (∂∂t+u‾∂∂x)p′+γp‾∂u′∂x=0 \left(\frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x} \right)p' + \gamma \overline{p} \frac{\partial u'}{\partial x} = 0 in order to get a dispersion relation

00 =(∂∂t+u‾∂∂x)p′+γp‾∂u′∂x= \left(\frac{\partial}{\partial t} + \overline{u} \frac{\partial}{\partial x} \right)p' + \gamma \overline{p} \frac{\partial u'}{\partial x}
=(−ikc+u‾ik)p^Ψ−γp‾ikρ‾(u‾−c)p^Ψ = \left(-ikc + \overline{u}ik \right)\hat{p}\Psi - \gamma \overline{p} \frac{ik}{\overline{\rho}(\overline{u} -c)} \hat{p}\Psi
=[(−c+u‾)−γp‾1ρ‾(u‾−c)]p^Ψ= \left[\left(-c + \overline{u} \right) - \gamma \overline{p} \frac{1}{\overline{\rho}(\overline{u} -c)} \right]\hat{p}\Psi

And by rearranging this we get c=u‾+γRdT‾ c = \overline{u} + \sqrt{\gamma R_d \overline{T}} which agrees with Lamb waves being acoustic

The conclusions for our data:

p′p' decays vertically like exp⁡(−∫z=0zgRdγT‾(z) dz) \exp\left(- \int_{z=0}^z\frac{g}{R_d \gamma \overline{T}(z)}\, \mathrm{d}z\right) . We know from the dcmip document that

Tv(z,ϕ)=1τ1(z)−τ2(z)IT(ϕ)T_v(z, \phi) = \frac{1}{\tau_1(z) - \tau_2(z)I_T(\phi)} IT(ϕ)=(cos⁡ϕ)K−KK+2(cos⁡ϕ)K+2I_T(\phi) = (\cos \phi)^K - \frac{K}{K+2}(\cos \phi)^{K+2} τint,1(z)=1Γ[exp⁡(ΓzT0)−1]+z(T0−TPT0TP)exp⁡[−(gzbRdT0)2]\tau_{\mathrm{int}, 1}(z) = \frac{1}{\Gamma}\left[\exp\left(\frac{\Gamma z}{T_0} \right) - 1 \right] + z\left(\frac{T_0-T_P}{T_0T_P}\right)\exp\left[-\left( \frac{gz}{bR_dT_0}\right)^2 \right] τint,2(z)=K+22(TE−TPTETP)zexp⁡[−(gzbRdT0)2]\tau_{\mathrm{int}, 2}(z) = \frac{K+2}{2}\left( \frac{T_E-T_P}{T_E T_P}\right)z\exp\left[ -\left( \frac{gz}{bR_dT_0}\right)^2\right]

Where we note that

exp⁡(−∫z=0zgRdγT‾(z) dz)\exp\left(- \int_{z=0}^z\frac{g}{R_d \gamma \overline{T}(z)}\, \mathrm{d}z\right) =exp⁡(−gRdγ∫z=0z[τ1(z)−τ2(z)IT(ϕ)] dz)= \exp\left(- \frac{g}{R_d \gamma} \int_{z=0}^z\left[\tau_1(z) - \tau_2(z)I_T(\phi)\right] \, \mathrm{d}z\right)
=exp⁡(−gRdγ[τint,1(z)−τint,2(z)IT(ϕ)]) = \exp\left(- \frac{g}{R_d \gamma} \left[\tau_{\mathrm{int}, 1}(z) - \tau_{\mathrm{int}, 2}(z)I_T(\phi)\right] \right)

and so we want to look for this kind of decay signature in the pressure field.

  • SE increased divergence damping --> look for figures
  • Check phase speed of lamb wave
  • FV3 0.25 is probably available
  • DCMIP 2012: gravity wave test case

Tracking down the actual speed of the wave:

We'll do this in python using Omega 850 masked so that it doesn't catch the baroclinic instability.

Method:

We're looking at time 6, 12, and 18 hours after start of run.

I'm calculating the maximum value of (signed) ω850\omega_{850} at each time stamp over a small region containing just the wave.

t∈{3,4,5}t\in \{3, 4, 5\} , t∈{6,7,8,9}t\in \{6, 7, 8, 9\} , t∈{10,11,12}t\in \{10, 11, 12\} ,
ϕbottom=\phi_{\mathrm{bottom}} = -90 -90 -90
ϕtop=\phi_{\mathrm{top}} = 90 90 90
λleft=\lambda_{\mathrm{left}} = 150 190 220
λright=\lambda_{\mathrm{right}} = 220 290

And locating the gridpoint at which ω850\omega_{850} is maximized gives

t=3t=3 , t=4t=4 , t=5t=5 t=6t=6 t=7t=7 t=8t=8 t=9t=9 t=10t=10 t=11t=11 t=12t=12
ϕargmax=\phi_{\mathrm{argmax}} = 35.1562 34.453 27.421 21.796 12.65 7.0312 0.0 -8.43 -13.35 -21.09
λargmax=\lambda_{\mathrm{argmax}} = 177.89 191.95 201.09 210.23 215.85 224.29 232.03 237.65 246.79 253.82

According to the formula for greatcircle distance dgc=aarccos⁡(sin⁡ϕ1sin⁡ϕ2+cos⁡ϕ1cos⁡ϕ2cos⁡(Δλ)). d_{gc} = a \arccos {\bigl (}\sin \phi _{1}\sin \phi _{2}+\cos \phi _{1}\cos \phi _{2}\cos(\Delta \lambda ){\bigr )}.

I used the following script

View track_gw_speed.py 
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
from os.path import join


fdir = "/nfs/turbo/cjablono2/owhughes/mountain_test_case_netcdf/"
fname = "ne30_1h_output.nc"

ne30_ds = xr.open_dataset(join(fdir, fname))
t =  12
t0 = 2

bounds = [[-90, 90, 150, 220], [-90, 90, 150, 220], [-90, 90, 150, 220], [-90, 90, 150, 220],
	  [-90, 90, 190, 260], [-90, 90, 190, 260], [-90, 90, 190, 260], [-90, 90, 190, 260],
	  [-90, 90, 220, 290], [-90, 90, 220, 290], [-90, 90, 220, 290] ]
lons = ne30_ds['lon']
lon_mask = np.logical_and(lons > bounds[ (t-t0)][2], lons < bounds[ (t-t0)][3])
lats = ne30_ds['lat']
lat_mask = np.logical_and(lats > bounds[ (t-t0)][0], lats < bounds[ (t-t0)][1]) 
print(ne30_ds['time'][2 * t])
omega850 = ne30_ds['OMEGA850'][2 * t, :, :] #- ne30_ds['OMEGA850'][1, :, :]
omega850 = omega850.where(lat_mask).where(lon_mask)
res = omega850.argmax(dim=("lat", "lon"))
latmax = (omega850.lat[res['lat']].values)
lonmax = (omega850.lon[res['lon']].values)

print(f"lat: {latmax}, lon:{lonmax}")
plt.figure()


ax = plt.axes(projection=ccrs.PlateCarree())

plt.contourf(lons.where(lon_mask), lats.where(lat_mask), omega850,
             transform=ccrs.PlateCarree())
plt.text(lonmax, latmax, 'Max',c="white",
         horizontalalignment='center',
         transform=ccrs.PlateCarree())

a = 6371 #km


plt.savefig(f"{t}_omega_test.pdf")
View comp_gw_speed.py 
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
from os.path import join

fdir = "/nfs/turbo/cjablono2/owhughes/mountain_test_case_netcdf/"
fname = "ne30_1h_output.nc"

ne30_ds = xr.open_dataset(join(fdir, fname))

a = 6371e3 #km
dt = 60 * 60
gamma = 1003/(1003 - 287.3)
Rd = 287.3

lambdamax = np.array([177.89, 191.95, 201.09, 210.23, 215.85, 224.29, 232.03, 237.65, 246.79, 253.82 ])
phimax = np.array([35.1562, 34.453, 27.421, 21.796, 12.65, 7.0312, 0.0, -8.43, -13.35, -21.09])
print(ne30_ds['lon'])
subset = ne30_ds.isel(indexers={"time": [0],
			       }).sel(indexers={"lon": lambdamax, "lat": phimax}, method="nearest")
lambdamax = np.deg2rad(lambdamax)
phimax = np.deg2rad(phimax)
T = subset["T850"]
u = subset["U850"]
print("predicted speeds: ")
print(u + np.sqrt(gamma * T * Rd))

gc = a * np.arccos(np.sin(phimax[1:]) * np.sin(phimax[:-1]) + np.cos(phimax[1:]) * np.cos(phimax[:-1]) * np.cos(lambdamax[1:] - lambdamax[:-1]))

print("calculated speeds: ")
print(gc/dt)

Ok so computing this c=u‾+γRdT‾c = \overline{u} + \sqrt{\gamma R_d \overline{T}} . Taking T‾=280K \overline{T} = 280 \mathrm{K} and in the lower atmosphere u‾≈10m/s \overline{u} \approx 10\mathrm{m}/\mathrm{s} , and γ=1003/(1003−287.3)=1.40 \gamma = 1003 / (1003 - 287.3) = 1.40 , and Rd=287.3J⋅kg−1⋅K−1 R_d = 287.3 \mathrm{J}\cdot \mathrm{kg}^{-1}\cdot \mathrm{K}^{-1} we therefore find that c=10m/s+1.4⋅280K⋅287.3J⋅kg−1⋅K−1=345.5m/s. c = 10\mathrm{m}/\mathrm{s} + \sqrt{1.4 \cdot 280 \mathrm{K} \cdot 287.3 \mathrm{J}\cdot \mathrm{kg}^{-1} \cdot \mathrm{K}^{-1}} = 345.5 \mathrm{m}/\mathrm{s}.

The calculated speeds of the wave from the scripts above vary from 356.95m/s 356.95 \mathrm{m}/\mathrm{s} at hour 3, and 316.35m/s. 316.35 \mathrm{m}/\mathrm{s}. at hour 12

This seems to be compatible. This is deeply curious to me, as it indicates that the wave may be acoustic?

Idea: increase temperature by 90 degrees, i.e. TE=400K, TP=330K T_E = 400\mathrm{K} ,\ T_P = 330 \mathrm{K}

Propagation under increased temperature

In this case latmax = [33.0468, 30.937, 22.5, 14.062, 7.031, 0.7031, -11.95, -19.68, -24.60, -32.34] , lonmax = [182.812, 197.57, 206.71, 215.85, 224.29, 234.84, 237.65, 246.79, 258.75, 269.29]

Using the same calculation as above, we find that the calculated speeds are between 391m/s391 \mathrm{m} / \mathrm{s} and 372m/s 372 \mathrm{m} / \mathrm{s} at time 6.

Under a 90∘C 90^\circ\mathrm{C} change in temperature, c=10m/s+1.4⋅370K⋅287.3J⋅kg−1⋅K−1=395.5m/s. c = 10\mathrm{m}/\mathrm{s} + \sqrt{1.4 \cdot 370 \mathrm{K} \cdot 287.3 \mathrm{J}\cdot \mathrm{kg}^{-1} \cdot \mathrm{K}^{-1}} = 395.5 \mathrm{m}/\mathrm{s}.

Without being very careful to check whether changing the two temperature parameters changes anything else, this indicates that we might expect about a +50m/s + 50 \mathrm{m}/\mathrm{s} change in wave speed. This once again seems consistent.

Something important to note: changing the temperature creates an increase in the strength of the zonal jet. This complicates things somewhat.

With T+45K T + 45\mathrm{K} : latmax = [33.04, 30.23, 24.60], lonmax = [180.0, 192.65, 203.90]

Increasing divergence damping in the SE model:

If this were actually a gravity wave, we would expect that increasing divergence damping would severely curtail the ability of the wave to propagate (see lecture notes on gravity waves). In order to test this we increase se_nu_div from 0.1×10160.1 \times 10^{16} (the CAM default) by a factor of 4 to 0.4×10160.4 \times 10^{16}

This is shown below, i.e. with default divergence hyperviscosity on top and with increased hyperviscosity below. Omega 850 day 0.25 Omega 850 day 0.75 Omega 850 day 1.25

This seems to show that increasing divergence damping does little to curtail the propagation of this wave. This is evidence that we are, in fact, not observing a gravity wave.

Next steps:

  • read lin again
  • read Vallis chapter
  • Look in holton?

Parent post: