2  方程式系と大気波動

ここでは、 (Arakawa and Konor 2009) に従って、大気中に存在しうる様々な波動について調べる。

2.1 音波と重力波

2.1.1 基礎方程式系

\(f\)平面上の等温大気（\(T=T_0\)）に対する次のような基礎方程式系を考える。

\[\begin{aligned} \frac{\mathrm{d}u}{\mathrm{d}t} &= \;\;\;fv - \frac{1}{\;\rho\;}\frac{\partial p}{\partial x} \\ \frac{\mathrm{d}v}{\mathrm{d}t} &= -fu - \frac{1}{\;\rho\;}\frac{\partial p}{\partial y} \\ \frac{\mathrm{d}w}{\mathrm{d}t} &=-\frac{1}{\;\rho\;}\frac{\partial p}{\partial z} - g \\ \frac{\mathrm{d}\ln\theta}{\mathrm{d}t} &= 0 \\ \frac{\mathrm{d}p}{\mathrm{d}t} &= -c_\mathrm{s}^2\rho\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right) \end{aligned}\]

線型化された方程式系（擾乱を表す\('\)は省略）は、次のようになる。

\[\begin{aligned} \frac{\partial u}{\partial t} &= \;\;\;fv-\frac{\partial}{\partial x}\left(\frac{p}{\;\overline{\rho}\;}\right) \\ \frac{\partial v}{\partial t} &= -fu-\frac{\partial}{\partial y}\left(\frac{p}{\;\overline{\rho}\;}\right) \\ \end{aligned}\] \[ \delta\frac{\partial w}{\partial t} = -\left(\frac{\partial}{\partial z}-\underline{\underline{\frac{\kappa}{H}}}\right)\frac{p}{\overline{\rho}}+g\frac{\theta}{\;\overline{\theta}\;} \tag{2.1}\] \[ \frac{\partial}{\partial t}\left(\frac{\theta}{\overline{\theta}}\right) = -\frac{\kappa}{H}w \] \[ \frac{1}{c_\mathrm{s}^2}\frac{\partial}{\partial t}\left(\underline{\frac{p}{\;\overline{\rho}\;}}\right) = -\frac{1}{\;\overline{\rho}\;}\left\{\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\left[\frac{\partial}{\partial z}+\frac{1}{H}(\underline{\underline{\kappa}}-1)w\right]\right\} \tag{2.2}\]

準静力学（quasi-hydrostatic、プリミティブ）系

静力学平衡を仮定し、(Equation 2.1) で\(\delta=0\)とする。

非弾性（anelastic）系

熱力学変数の水平一様成分からのずれが小さいと仮定し、連続の式で密度の局所時間微分を無視する(Ogura and Phillips 1962)。 上の方程式系では、(Equation 2.1) 及び (Equation 2.2) の下線及び二重下線を無視する。

擬非圧縮（pseudo-incompressible）系

連続の式における気圧の局所時間微分を無視するので、(Equation 2.2) の下線を無視する(Durran 1989)。

統一（unified）系

連続の式の密度を静力学の密度に置き換える。 気圧を静力学成分とそれからのずれに分ける(Arakawa and Konor 2009)。 \[ p = p_\mathrm{qs}+ p' \] 以下\(p_\mathrm{qs}\)は単に\(p\)と書く。

簡単のため\(y\)方向に一様であるとし、

\[\begin{aligned} b &\equiv g\frac{\theta}{\;\overline{\theta}\;} \\ N^2 &\equiv g \frac{\mathrm{d}\ln\overline{\theta}}{\mathrm{d}z}\\ \frac{\mathrm{d}\ln\overline{\theta}}{\mathrm{d}z}&= \frac{\kappa}{\;H\;} \end{aligned}\]

を用いると、

\[ \frac{\partial^2u}{\partial t^2} = -f^2u-\frac{\partial^2}{\partial t\partial x}\left(\frac{p}{\;\overline{\rho}\;}+\frac{p'}{\;\overline{\rho}\;}\right) \tag{2.3}\] \[ \delta \frac{\partial w}{\partial t} = -\left(\frac{\partial}{\partial z} - \underline{\underline{\frac{\kappa}{\;H\;}}}\right)\left(\frac{p'}{\;\overline{\rho}\;}\right) \tag{2.4}\] \[ \frac{\partial b}{\partial t} = -N^2w \tag{2.5}\] \[ \left(\frac{\partial}{\partial z} - \underline{\underline{\frac{\kappa}{\;H\;}}}\right)\left(\frac{p}{\;\overline{\rho}\;}\right) = b \] \[ \frac{1}{c_\mathrm{s}^2}\frac{\partial}{\partial t}\left(\underline{\frac{p}{\;\overline{\rho}\;}+\varepsilon\frac{p'}{\;\overline{\rho}\;}}\right) = -\left\{\frac{\partial u}{\partial x}+\left[\frac{\partial}{\partial z}+\frac{1}{\;H\;}(\underline{\underline{\kappa}}-1)\right]w\right\} \tag{2.6}\]

が得られる。

統一系では、(Equation 2.2) において\(\epsilon=0\)とする。

表 Table 2.1 に方程式系を切り替える「スイッチ」をまとめる。

Table 2.1: 方程式系と\(\varepsilon\), \(\delta\), 下線、二重下線の関係

	\(\varepsilon\)	\(\delta\)	下線
完全圧縮系	\(1\)	\(1\)	残す
統一系	\(0\)	\(1\)	残す
擬非圧縮系	N/A	\(1\)	下線を無視
非弾性系	N/A	\(1\)	下線及び二重下線を無視
準静力学系	\(1\)	\(0\)	残す

2.1.2 波動解

波動解を仮定する。

\[\begin{aligned} \pmb{v} &= (u,w,b,\frac{p}{\;\overline{\rho}\;},\frac{p'}{\;\overline{\rho}\;})^\mathrm{T} \\ \pmb{u} &= (U,W,B,P,Q)^\mathrm{T} \\ \sqrt{\overline{\rho}}\pmb{v} &= \Re[\pmb{u}\exp(kx+mz-\omega t)] \end{aligned}\]

式 (Equation 2.3)–(Equation 2.6) は行列 \(M\) を用いて表すことができる。 \[ M\pmb{u} = 0 \]

\[\begin{aligned} (f^2-\omega^2)U-k\omega(P+Q) &= 0 \\ \delta\omega W + (-m+i\mu)Q &= 0 \\ \omega B + iN^2W &= 0 \\ -iB+(-m+i\mu)P = 0 \\ \frac{1}{c_\mathrm{s}^2}\omega(\underline{P+\varepsilon Q}) &=0 \end{aligned}\]

ここで

\[ \mu \equiv \left(\frac{1}{\;2\;}-\underline{\underline{\kappa}}\right)\frac{1}{\;H\;} \]

\[\begin{aligned} \frac{1}{\;\overline{\rho}\;}\frac{\mathrm{d}\overline{\rho}}{\mathrm{d}z} &= \frac{1}{\overline{\rho}RT_0}\frac{\partial\overline{p}}{\partial z} = -\frac{1}{H} \\ \frac{\partial}{\partial z} &= im - \frac{1}{\;2\;}\frac{1}{\;\overline{\rho}\;}\frac{\mathrm{d}\overline{\rho}}{\mathrm{d}z} = \left(im+\frac{1}{2H}\right) \\ \frac{\partial}{\partial z} - \frac{\kappa}{\;H\;}&= im + \left(\frac{1}{\;2\;}-\kappa\right)\frac{1}{\;H\;} = im + \mu \\ %\pd{}{z} + (\kappa-1)\fracl{1}{H} &= im - \mu \end{aligned}\]

\[ M=\begin{pmatrix} f^2-\omega^2 & 0 & 0 & k\omega & k\omega \\ 0 & \delta\omega & 0 & 0 & -m+i\mu \\ 0 & iN^2 & \omega & 0 & 0 \\ -k & -m-i\mu & 0 & \underline{\dfrac{\omega}{c_\mathrm{s}^2}} & \underline{\varepsilon\dfrac{\omega}{c_\mathrm{s}^2}} \\ 0 & 0 & -i & -m+i\mu & 0 \end{pmatrix} \tag{2.7}\]

2.1.3 順圧モード

\(b=0\) とすると、(Equation 2.7) は縮退し、その行列式は次のようになる。

\[\begin{aligned} \det M &= \begin{vmatrix} f^2-\omega^2 & k\omega \\ -k & \underline{\dfrac{\omega}{c_\mathrm{s}^2}} \end{vmatrix} \\ &= \underline{\dfrac{\omega}{c_\mathrm{s}^2}} (f^2-\omega^2) + k^2\omega = 0 \end{aligned}\]

各方程式系の分散関係式は次のようになる。

非弾性系、擬非圧縮系

自明な解のみ。 \[\omega=0\]

完全圧縮系、統一系、準静力学系

コリオリ力で変形されたLamb波。 \[\omega^2 = f^2 + k^2c_\mathrm{s}^2\]

2.1.4 傾圧モード

\(\det M = 0\) より次の式を得る。

\[\begin{aligned} \underline{\varepsilon\delta\omega^4}&-c_\mathrm{s}^2\left(\delta k^2 + m^2 + \mu^2 + \underline{\frac{N^2+\varepsilon\delta f^2}{c_\mathrm{s}^2}}\right)\omega^2 \nonumber \\ &+c_\mathrm{s}^2\left[f^2\left(m^2+\mu^2+\underline{\frac{N^2}{c_\mathrm{s}^2}}\right)+k^2N^2\right] = 0 \end{aligned} \tag{2.8}\]

\(\omega^2\) の2次方程式 (Equation 2.8) を解くと分散関係式を得る。

\[\begin{aligned} \omega^2 = &\frac{c_\mathrm{s}^2}{2}\left(\delta k^2 + m^2 + \mu^2 + \underline{\frac{N^2+\varepsilon\delta f^2}{c_\mathrm{s}^2}}\right) \nonumber \\ &\pm \frac{c_\mathrm{s}^2}{2}\left\{\left(\delta k^2 + m^2 + \mu^2 + \underline{\frac{N^2+\varepsilon\delta f^2}{c_\mathrm{s}^2}}\right)^2\right. \nonumber \\ &\left.-4\left[f^2\left(m^2+\mu^2+\underline{\frac{N^2}{c_\mathrm{s}^2}}\right)+k^2N^2\right]\right\}^{1/2} \end{aligned}\]

正の根は音波を表す。

\[ \omega^2 = c_\mathrm{s}^2\left(k^2 + m^2 + \mu^2 + \underline{\frac{N^2+\varepsilon\delta f^2}{c_\mathrm{s}^2}}\right),\; \delta=\varepsilon=1 \]

負の根は重力波を表す。

\[ \omega^2 = \frac{f^2\left(m^2+\mu^2+\underline{\dfrac{N^2}{c_\mathrm{s}^2}}\right)+k^2N^2}{\delta k^2 + m^2 + \mu^2 + \underline{\dfrac{N^2+\varepsilon\delta f^2}{c_\mathrm{s}^2}}} \]

各方程式系に対する分散関係式は次のようになる。

完全圧縮系

\(\varepsilon=\delta=1\)とし、下線の項は残す。 \[\begin{equation} \omega^4 -c_\mathrm{s}^2\left(k^2 + m^2 + \mu^2 + \frac{N^2+f^2}{c_\mathrm{s}^2}\right)\omega^2 + c_\mathrm{s}^2\left[f^2\left(m^2+\mu^2+\frac{N^2}{c_\mathrm{s}^2}\right)+k^2N^2\right] = 0 \end{equation}\]

統一系

\(\varepsilon = 0, \delta=1\)とし、下線の項は残す。 \[\omega^2 = \frac{f^2\left(m^2+\mu^2+\dfrac{N^2}{c_\mathrm{s}^2}\right)+k^2N^2}{k^2 + m^2 + \mu^2 + \dfrac{N^2+f^2}{c_\mathrm{s}^2}}\]

擬非圧縮系

\(\delta=1\)とし、下線の項を省略する。 \[\omega^2 = \frac{f^2\left(m^2+\mu^2\right)+k^2N^2}{k^2 + m^2 + \mu^2}\]

非弾性系

\(\delta=1\) とし、下線及び二重下線の項を省略する。 \(\mu = 1/2H\) となる。 \[\omega^2 = \frac{f^2\left(m^2+\dfrac{1}{4H^2}\right)+k^2N^2}{k^2 + m^2 + \dfrac{1}{4H^2}}\]

準静力学系

\(\delta=0\) とし、下線の項は保持する。 \[\omega^2 = \frac{f^2\left(m^2+\mu^2+N^2/c_\mathrm{s}^2\right)+k^2N^2}{m^2 + \mu^2 + \dfrac{N^2+f^2}{c_\mathrm{s}^2}}\]

2.1.5 分散関係の図示

分散関係式において、パラメタ \(L_x, n, H, D, \gamma, g\) に具体的な値を与える。

\[\begin{aligned} H = 10\,\mathrm{km},\; D = 100\,\mathrm{km},\;g=9.8\,\mathrm{m/s} \\ L_x = \frac{2\pi}{k}\;\mathrm{m} \\ m = \frac{\pi n}{D},\; n= 1, 5, 10, 20, 40, 80, 160 \\ c_\mathrm{s}^2= \gamma R T_0 = \gamma g H = (370\,\mathrm{m/s})^2\\ \mu = \left(\frac{1}{\;2\;}-\underline{\underline{\kappa}}\right)\frac{1}{\;H\;} \\ \kappa = \frac{R}{\;c_p\;} = 1-\frac{1}{\gamma},\; \gamma = \frac{c_p}{\;c_v\;} = \frac{7}{\;5\;} \\ N^2 = g\frac{\mathrm{d}\ln\overline{\theta}}{\mathrm{d}z}= \frac{g\kappa}{H} = 2.8\times10^{-4}\,\mathrm{s}^{-1} \end{aligned}\]

これらを用いて求めた分散関係を Figure 2.1 に示す。

# constants
gamma <- 7 / 5
kappa <- 1 - 1 / gamma
g <- 9.8
a <- 6.371e6
day_in_sec <- 86400
om <- 2 * pi / day_in_sec

# list for systems with switches
# d, e, u1, u2
# d: delta
# e: epsilon
# u1: underlined
# u2: double underlined
full <-      c(1, 1, 1, 1)
unified <-   c(1, 0, 1, 1)
pseudo <-    c(1, 0, 0, 1)
anelastic <- c(1, 0, 0, 0)
qhydro <-    c(0, 1, 1, 1)
systems <- data.frame(full, unified, pseudo, anelastic, qhydro)

# parameters
H <- 10.0e3           # scale height m
D <- 100.0e3          # top of the atmosphere m
cs2 <- gamma * g * H  # speed of sound squared
N2 <- g * kappa / H   # Brunt-Väisälä frequency squared
lat <- 40
f <- 2 * om * sin(lat / 180 * pi)
f2 <- f ** 2
beta <- 2 * om * cos(lat / 180 * pi) / a

nx <- 61
nmode <- 7

l <- 10 ** seq(2, 8, length.out = nx)
k <- 2 * pi / l
l <- l * 1.0e-3
n <- c(1, 5, 10, 20, 40, 80, 160)
m <- pi * n / D
omega <- matrix(NA, nmode, nx)

mu <- function(dk) {
  (0.5 - dk * kappa) / H
}

acoustic <- function(k, m) {
  sqrt(cs2 * (k^2 + m^2 + mu(1.0)^2) + (N2 + f2))
}

gravity <- function(k, m, d, e, u1, u2) {
  mu2 <- mu(u2)^2
  k2 <- k^2
  mmn <- m^2 + mu2 + u1 * N2 / cs2
  sqrt((f2 * mmn + k2 * N2) / (d * k2 + mmn + e * d * u2 * f2 / cs2))
}


barotropic <- function(k, u1) {
  -beta * k / (k^2 + u1 * f2 / cs2)
}

baroclinic <- function(k, m, u1, u2) {
  -beta * k / (k^2 + f2 * (u1 / cs2 + H / (kappa * g) * (m^2 + mu(u2)^2)))
}

plot_dispersion <- function(system, d, e, u1, u2) {
  uu <- u1 * u2
  fc <- d * e * uu

  # empty plot
  plot(NULL, NULL, log = "xy",
       xlim = c(tail(l, 1), l[1]),
       ylim = c(5e-5, 1),
       xlab = expression(L[x] ~ "km"),
       ylab = expression(abs(omega) ~ "s"^-1),
       main = system)

  # Lamb wave
  if (uu == 1) {
    omega[1, ] = sqrt(f^2 + k^2 * cs2)
    lines(l, omega[1, ], col = "black", lwd = 2, lty = 2)
  }

  # Sound waves
  if (fc == 1) {
    for (j in 1:nmode) {
      omega[j, ] <- acoustic(k, m[j])
      lines(l, omega[j, ], col = "black", lwd = 2, lty = 1)
    }
  }

  # Gravity waves
  for (j in 1:nmode) {
    omega[j, ] <- gravity(k, m[j], d, e, u1, u2)
    lines(l, omega[j, ], col = "black", lwd = 2, lty = 1)
  }
}

par(mfrow = c(3, 2), mar = c(4, 4, 3, 2))
panels <- c("(a)", "(b)", "(c)", "(d)", "(e)")
snames <- c("Fully-Compressible", "Unified",
  "Pseudo-Compressible", "Anelastic", "Quasi-Hydrostatic")
names(snames) <- names(systems)

i <- 1
for (sys in names(systems)) {
  deu1u2 <- systems[[sys]]
  d <- deu1u2[1]
  e <- deu1u2[2]
  u1 <- deu1u2[3]
  u2 <- deu1u2[4]
  
  full_title <- paste(panels[i], snames[sys])
  plot_dispersion(full_title, d, e, u1, u2)

  i <- i + 1
}

Figure 2.1: 平面上の分散関係式。(a) 完全圧縮系、(b) 統一系、(c) 擬非圧縮系、(d) 非弾性系, (e) 準静力学系。破線はLamb波、左上は音波、右下は重力波の分散関係

2.2 準地衡風方程式系

この節では、 \(\beta\) 平面上の準地衡風方程式系で静力学平衡を仮定して、ロスビー波の分散関係式を導出する。

2.2.1 基礎方程式系

\[ f_0\frac{\partial v}{\partial x} = \frac{\partial^2}{\partial x^2}\frac{p}{\overline{\rho}} \tag{2.9}\]

\[\begin{aligned} \frac{\partial}{\partial t}\frac{v}{x} &= \beta v - f_0\frac{\partial u}{\partial x} \\ \left(\frac{\partial}{\partial z}-\frac{\kappa}{H}\right)\frac{p}{\overline{\rho}} &= \frac{\theta}{\overline{\theta}}g \\ \frac{\partial}{\partial t}\frac{\theta}{\overline{\theta}} &= -\frac{\kappa}{H}w \\ \end{aligned}\]

\[ \frac{1}{c_\mathrm{s}^2}\frac{\partial}{\partial t}\frac{p}{\overline{\rho}} = -\left\{\frac{\partial u}{\partial x} + \left[\frac{\partial}{\partial z} + \frac{1}{H}(\kappa-1)\right]w\right\} \tag{2.10}\]

2.2.2 波動解

波動解を仮定する。

\[\begin{aligned} \pmb{v} &= (u,v,w,\frac{p}{\;\overline{\rho}\;})^\mathrm{T} \\ \pmb{u} &= (U,V,W,P)^\mathrm{T} \\ \sqrt{\overline{\rho}}\pmb{v} &= \Re[\pmb{u}\exp(kx+mz-\omega t)] \end{aligned}\]

(Equation 2.10)–(Equation 2.10) に代入すると、次を得る。

\[\begin{aligned} ikf_0V + k^2 P &=0 \\ ikf_0U+(\beta + k\omega)V &= 0 \\ kU + (m+i\mu)W - \frac{1}{\;c_\mathrm{s}^2\;}\omega P &= 0\\ \frac{\kappa g}{H}W + \omega(m - i\mu) P &= 0 \end{aligned}\]

2.2.3 順圧モード

完全圧縮系、統一系、準静力学系では、順圧ロスビー波の分散関係式は次の式で与えられる。

\[ \omega = -\frac{k\beta}{k^2+\underline{\left(\dfrac{f_0}{c_\mathrm{s}}\right)^2}} \]

波数 \(k\rightarrow 0\) のとき \(\omega\rightarrow 0\) である。

擬非圧縮系及び非弾性系では、下線の項が無いので、

\[ \omega = -\frac{\beta}{\;k\;} \]

となる。 波数 \(k\rightarrow 0\) のとき \(\omega\rightarrow -\infty\) となる。

2.2.4 傾圧モード

傾圧ロスビー波の分散関係式は、次のようになる。

\[ \omega = -\frac{k\beta}{k^2+f_0^2\left[\underline{\dfrac{1}{c_\mathrm{s}^2}}+\dfrac{H}{\kappa g}(m^2+\mu^2)\right]} \tag{2.11}\]

モデルの上端の高さ \(D\) がスケールハイト \(H\) 程度であれば、分母の括弧内の第2項が支配的で、 \(D\) が \(H\) の10倍でも無視できないので、順圧ロスビー波に比べて擬非圧縮系や非弾性系による変形は小さい。

plot_dispersion_qg <- function(system, u1, u2) {
  # empty plot
  plot(NULL, NULL, 
       log = "xy",                  # 両対数グラフ (loglog)
       xlim = c(tail(l, 1), l[1]),  # x軸の範囲（逆向き）
       ylim = c(1.0e-2, 1.0e-8),    # y軸の範囲（逆向き仕様）
       xlab = expression(L[x] ~ "km"), 
       ylab = expression(abs(omega) ~ "s"^-1), 
       main = system)

  # Barotropic
  omega[1, ] <- barotropic(k, u1)
  lines(l, -omega[1, ], col = "black", lty = "dashed", lwd = 2)

  # Baroclinic
  for (j in 1:nmode) {
    omega[j, ] <- baroclinic(k, m[j], u1, u2)
    lines(l, -omega[j, ], col = "black", lty = "solid", lwd = 2)
  }
}

par(mfrow = c(2, 2), mar = c(4, 4, 3, 2))
panels <- c("(a)", "(b)", "(c)")
snames <- c("Fully-Compressible", "Pseudo-Compressible", "Anelastic")
names(snames) <- names(systems)[-which(names(systems) %in%
                                       c("unified", "qhydro"))]
i <- 1
for (sys in names(systems)) {
  if (sys == "unified" || sys == "qhydro") next

  deu1u2 <- systems[[sys]]
  u1 <- deu1u2[3]
  u2 <- deu1u2[4]
  
  full_title <- paste(panels[i], snames[sys])
  plot_dispersion_qg(full_title, u1, u2)

  i <- i + 1
}

Figure 2.2: 完全圧縮系、統一系、準静力学系、(b) 擬非圧縮系、(c)非弾性系。破線は順圧、実線は傾圧ロスビー波の分散関係

2.2.5 分散関係の図示

Section 2.1.5 と同じパラメタを用いて求めた分散関係を Figure 2.2 に示す。

2.3 まとめ

• 統一系は音波を取り除くが、Lamb波、重力波、ロスビー波を変形しない。
• 擬非圧縮系は深い重力波を少し変形し、準静力学系はこれを大きく変形する。
• 擬非圧縮系と非弾性系は、波長の長いロスビー波を変形する。

Note課題

• 高解像度モデルでは、準静力学系はどのくらいの水平解像度まで用いることができるか。
• 擬非圧縮系や非弾性系を全球モデルに用いるとどのような誤差が生じるか。
• 音波を取り扱うためには、音波を取り除く方法以外にどのようなものがあるか。
• 式 (Equation 2.11) を導出せよ。

Arakawa, A., and C. S. Konor, 2009: Unification of the anelastic and quasi-hydrostatic systems of equations. Mon. Wea. Rev., 137, 710–726, https://doi.org/10.1175/2008MWR2520.1.

Durran, D. R., 1989: Improving the anelastic approximation. J. Atmos. Sci., 46, 1453–1461, https://doi.org/10.1175/1520-0469(1989)046<1453:ITAA>2.0.CO;2.

Ogura, Y., and N. A. Phillips, 1962: Scale analysis of deep and shallow convection in the atmosphere. J. Atmos. Sci., 19, 173–179, https://doi.org/10.1175/1520-0469(1962)019<0173:SAODAS>2.0.CO;2.