9  低緯度の循環

この章では、赤道付近での熱源に対する応答を Gill (1980) に基づいて述べる。

9.1 変数分離

\(\log p\) 座標系での摂動方程式系は次のように書ける。

\[ \frac{\partial u}{\partial t} - fv = -\frac{\partial \phi}{\partial x} \tag{9.1}\]

\[ \frac{\partial v}{\partial t} + fu = -\frac{\partial \phi}{\partial y} \tag{9.2}\]

\[ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{1}{\rho_0}\frac{\partial \rho_0 w}{\partial z} = 0 \tag{9.3}\]

\[ \frac{\partial}{\partial t}\left(\frac{\partial \phi}{\partial z}\right) + N^2 w = 0 \tag{9.4}\]

Equation 9.3 と Equation 9.4 から \(w\) を消去すると

\[ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = \frac{1}{\rho_0}\frac{\partial }{\partial z} \left[\frac{\rho_0}{N^2}\frac{\partial}{\partial z} \left(\frac{\partial \phi}{\partial t}\right)\right] \tag{9.5}\]

となる。 Equation 9.5 の右辺が \(\partial\phi/\partial t\) に比例すると仮定し、 次のように水平と鉛直とを分離する。

\[ \begin{pmatrix} u\\ v\\ \phi \end{pmatrix} = \begin{pmatrix} u_n\\ v_n\\ \phi_n \end{pmatrix} (x,\,y,\,t)\hat{\phi}(z) \]

これを Equation 9.1, Equation 9.2, Equation 9.5 に代入すると水平構造方程式

\[ \frac{\partial u_n}{\partial t} - fv_n = -\frac{\partial \phi_n}{\partial x} \tag{9.6}\]

\[ \frac{\partial v_n}{\partial t} + fu_n = -\frac{\partial \phi_n}{\partial y} \tag{9.7}\]

\[ \frac{\partial \phi_n}{\partial t} + gh\left(\frac{\partial u_n}{\partial x} + \frac{\partial v_n}{\partial y}\right) = 0 \tag{9.8}\]

及び鉛直構造方程式

\[ \frac{1}{\rho_0}\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{\rho_0}{N^2} \frac{\mathrm{d}\hat{\phi}}{\mathrm{d}z}\right) = -\frac{1}{gh}\hat{\phi} \]

を得る。 \(h\) は等価深度と呼ばれている。

9.2 熱源に対する大気の応答

ここでは、(Gill 1980) に従って赤道付近で局在化された熱源に対する大気の応答を考える。 鉛直方向には傾圧第一モード（ \(n=1\) ）のみを考え、 Section 9.1 で求めた水平構造方程式を用いる。 赤道 \(\beta\) 平面近似（ \(f=\beta y, f_0 = 0\) ）を仮定し、 水平スケールを \(l_n \equiv\sqrt{c/2\beta},\,c=\sqrt{gh}\) 、 時間を \(l_n/c\) 、 \((u, v), \phi\) をそれぞれ \(c,\,c^{2}\) でスケールし \(\phi\) を \(p\) と書くと、 無次元化された方程式系

\[ \frac{\partial u}{\partial t} - \frac{1}{2}yv = -\frac{\partial p}{\partial x} \]

\[ \frac{\partial v}{\partial t} + \frac{1}{2}yu = -\frac{\partial p}{\partial y} \tag{9.9}\]

\[ \frac{\partial p}{\partial t} + \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = -Q \] \[ w = \frac{\partial p}{\partial t} + Q \]

を得る。 ここで \(Q\) は加熱率に比例し、正の時に加熱、 \(u,\,v,\, p\) の符号は地表面での符号を表す。

レイリー摩擦とニュートン冷却を仮定し、その係数 \(\varepsilon\) は等しいとする。 定常解は \(\partial/\partial t + \varepsilon\) を \(\varepsilon\) で置き換える。 Equation 9.9 の \(\varepsilon v\) は \(\varepsilon\) が小さく、 強制の東西スケールが \(2\varepsilon\) に比べて小さいこと、すなわち \(2\varepsilon k \ll 1\) より無視できる。

\[ \varepsilon u - \frac{1}{2}yv = -\frac{\partial p}{\partial x} \tag{9.10}\]

\[ \frac{1}{2}yu = -\frac{\partial p}{\partial y} \tag{9.11}\]

\[ \varepsilon p + \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = -Q \tag{9.12}\]

\(q = p + u,\,r = p - u\) とおく。 Equation 9.12 \(+\) Equation 9.10 より

\[ \varepsilon q + \frac{\partial q}{\partial x} + \frac{\partial v}{\partial y} - \frac{1}{2}yv = -Q \tag{9.13}\]

Equation 9.12 \(-\) Equation 9.10 より

\[ \varepsilon r - \frac{\partial r}{\partial x} + \frac{\partial v}{\partial y} + \frac{1}{2}yv = -Q \tag{9.14}\]

Equation 9.11 より

\[ \frac{\partial q}{\partial y} + \frac{1}{2}yq +\frac{\partial r}{\partial y} - \frac{1}{2}yr = 0 \tag{9.15}\]

が得られる。 ところで、エルミート多項式 \(H_n(y)=(-1)^ne^{y^2/2}\frac{\mathrm{d}^n}{\mathrm{d}y^n} e^{-y^2/2}\) を用いて、放物柱函数は

\[ D_n(y) \equiv e^{-y^2/4}H_n(y) \tag{9.16}\]

と定義され、次のような性質がある。

\[ \frac{\mathrm{d}D_n}{\mathrm{d}y} + \frac{1}{2}yD_n = nD_{n-1} \tag{9.17}\] \[ \frac{\mathrm{d}D_n}{\mathrm{d}y} - \frac{1}{2}yD_n =-D_{n+1} \tag{9.18}\]

\(D_n(y)\) の定義 Equation 9.16 と Equation 9.17, Equation 9.18 を用いると

\[ (D_0, D_1, D_2, D_3) = \left[1, y, y^2 - 1, y(y^2-3)\right]e^{1/4y^2} \tag{9.19}\]

と具体的な形が得られる。 次のように Equation 9.13, Equation 9.14, Equation 9.15 の \(q, r, v, Q\) が \(y\) 方向に \(D_n(y)\) で展開できるとする。

\[ \begin{pmatrix} q\\r\\v\\Q \end{pmatrix} =\sum_{n=0}^\infty \begin{pmatrix} q_n(x)\\r_n(x)\\v_n(x)\\Q_n(x) \end{pmatrix} D_n(y) \tag{9.20}\]

Equation 9.20 を Equation 9.13, Equation 9.14, Equation 9.15 に代入し、 Equation 9.17,Equation 9.18 を用いると

\[ \varepsilon q_0 + \frac{\mathrm{d}q_0}{\mathrm{d}x} = -Q_0 \tag{9.21}\]

\[ q_1 = 0 \tag{9.22}\]

\[ (2n+1)\varepsilon q_{n+1} - \frac{\mathrm{d}q_{n+1}}{\mathrm{d}x} = -Q_{n-1} - nQ_{n+1},\; n\ge 0 \tag{9.23}\]

が得られる。 \(Q_0\) または \(Q_1\) だけを考えるので、以下 Equation 9.23 の右辺の \(nQ_{n+1}\) は考慮しない。

熱源 \(Q_0\) または \(Q_1\) の東西方向の構造は

\[ F(x) = \begin{cases} \cos kx & |x|<L \\ 0 & |x|>L \end{cases} \]

とする。 ここで \(k = \pi/2L\) である。

Equation 9.21 の解は

\[ q_0 = \begin{dcases} 0 & x \le -L\\ -\frac{1}{\varepsilon^2 + k^2}\left[\varepsilon\cos kx + k\sin kx + ke^{-\varepsilon(x+L)}\right] & x < |L|\\ -\frac{k}{\varepsilon^2 + k^2}\left\{1 + \exp(-2\varepsilon L)\right\}e^{\varepsilon(L-x)} & x \ge L \end{dcases} \]

Equation 9.23 の解は

\[ q_{n+1} = \begin{dcases} -\frac{k}{(2n+1)^2\varepsilon^2 + k^2}e^{(2n+1)\varepsilon(x+L)}\left[1 + e^{-2(2n+1)\varepsilon L}\right] & x \le -L\\ -\frac{1}{(2n+1)^2\varepsilon^2 + k^2}\left[(2n+1)\varepsilon\cos kx - k\sin kx + ke^{(2n+1)\varepsilon(x-L)}\right] & x < |L|\\ 0 & x \ge L\\ \end{dcases} \]

\(q_0,\,q_{n+1}\) の東西平均値は

\[ \int_{-\infty}^{\infty}(q_0, q_{n+1})\mathrm{d} x = -\frac{1}{\varepsilon}\left(1, \frac{1}{2n+1}\right)I \]

ここで

\[ I = \int_{-\infty}^{\infty}F(x)\mathrm{d} x = \frac{2}{k} = \frac{4L}{\pi} \]

9.2.1 赤道に対称な熱源に対する応答

赤道に対称な熱源 \(Q_0 = F(x)\)

\[ Q(x,y) = Q_0D_0(y) = F(x)e^{-y^2/4} \tag{9.24}\]

について微分方程式 Equation 9.21–Equation 9.23 を解く。

赤道対称な熱源に対する応答は二つに分かれる。 そのひとつは \(q_0\) だけで決まり、東進するとともに減衰するケルビン波を表す。 位相速度は \(1\) で減衰率は \(\varepsilon\) である。 西に情報が伝わらないので、 \(x<-L\) で振幅は0である。

\[\begin{aligned} u &= p = \frac{1}{2}q_0(x)e^{-y^2/4}\\ v &= 0\\ w &= \dfrac{1}{2}\left[\varepsilon q_0(x) + F(x)\right]e^{-y^2/4} \end{aligned}\]

この解は強制域で上昇し、上空で西風、下層で東風が赤道に沿って吹くウォーカー循環を表す。 これに対応する赤道に沿った低圧部があり、強制域の方に向かって気圧が西ほど低くなっている。

Code 
L <- 2.0
k <- pi / (2.0 * L)
eps <- 0.1
nx <- 251
ny <-  81
nz <-  30
xmin <- -10.0
ymin <-  -4.0
ymax <-   4.0
xmax <-  15.0
zmin <-   0.0
zmax <-    pi
I <- 4.0 * L / pi

lon <- seq(xmin, xmax, length.out = nx)
lat <- seq(ymin, ymax, length.out = ny)
lev <- seq(zmin, zmax, length.out = nz)

q0 <- function(x) {
  ekr <- 1.0 / (eps^2 + k^2)
  kx <- k * x
  ifelse(x <= L, 
    ifelse(x <= -L, 0.0,
      -ekr * (eps * cos(kx) + k * (sin(kx) + exp(-eps * (x + L))))
    ), -ekr * k * (1.0 + exp(-2.0 * eps * L)) * exp(eps * (L - x))
  )
}

qn <- function(n, x) {
  nn <- 2 * n[1] + 1
  ekr <- 1.0 / ((nn * eps)^2 + k^2)
  kx <- k * x
  ifelse(x <= L,
    ifelse (x <= -L,
      -ekr * k * (1.0 + exp(-2 * nn * eps * L)) * exp(nn * eps * (x + L)),
       ekr * (-nn * eps * cos(kx) + k * (sin(kx) - exp(nn * eps * (x - L))))
    ), 0.0
  )
}

ffunc <- function(x) {
  ifelse(abs(x) < L, cos(k * x), 0.0)
}

symmetric <- function(x, y) {
  q0x <- q0(x)
  q2x <- qn(1, x)
  Fx <- ffunc(x)

  y2 <- y^2
  ey2 <- exp(-0.25 * y2)

  p <- 0.5 * (outer(q0x, ey2) + outer(q2x, (1 + y2) * ey2))
  u <- 0.5 * (outer(q0x, ey2) + outer(q2x, (y2 - 3) * ey2))
  v <- outer(Fx + 4 * eps * q2x, y * ey2)
  w <- outer(0.5 * eps * q0x + Fx, ey2) + 0.5 * eps * outer(q2x, (1 + y2) * ey2)

  list(p = p, u = u, v = v, w = w)
}

antisymmetric <- function(x, y) {
  q3x <- qn(2, x)
  Fx <- ffunc(x)

  y2 <- y^2
  y3 <- y2 * y
  ey2 <- exp(-0.25 * y2)

  p <- 0.5 * outer(q3x, y3 * ey2)
  u <- 0.5 * outer(q3x, (y3 - 6 * y) * ey2)
  v <- outer(6 * eps * q3x,  (y2 - 1) * ey2) + outer(Fx, y2 * ey2)
  w <- 0.5 * eps * outer(q3x, y3 * ey2) + outer(Fx,  y * ey2)

  list(p = p, u = u, v = v, w = w)
}

# kai: clockwise + (meteorology, oceanography as opposed to fluid dynamics)
hadley_symmetric <- function(y, z) {
  y2 <- y^2
  eyy <- I * exp(-0.25 * y2)

  p <- -(4 + y2) / (6 * eps) * eyy

  u <- outer(-y2 / (6 * eps) * eyy, cos(z))
  kai <- outer(y / 3 * eyy, sin(z))

  list(p = p, u = u, kai = kai)
}

hadley_antisymmetric <- function(y, z) {
  y2 <- y^2
  y3 <- y2 * y
  eyy <- I * exp(-0.25 * y2)

  p <- -y3 / (10 * eps) * eyy
  u <- outer((6 * y - y3) / (10 * eps) * eyy, cos(z))
  kai <- outer((y2 - 6) / 5 * eyy, sin(z))

  list(p = p, u = u, kai = kai)
}

walker <- function(x, z) {
  q0x <- q0(x)
  q2x <- qn(1, x)

  p <- (q0x + 3 * q2x) * sqrt(pi)
  kai <- outer((q2x - q0x) * sqrt(pi), sin(z))

  list(p = p, kai = kai)
}
Code 
drawvec <- function(x, y, u, v, length = 0.05, uvref = NA, ...) {
  nx <- length(x)
  ny <- length(y)
  usr <- par("usr")
  d <- min((usr[2] - usr[1]) / nx, (usr[4] - usr[3]) / ny)
  if (is.na(uvref)) uvref <- max(abs(c(u, v)))
  d <- d / uvref
  u <- 0.5 * u * d
  v <- 0.5 * v * d
  X <- rep(x, times = length(y))
  Y <- rep(y, each  = length(x))
  tol <- 1e-8
  suppressWarnings(
    arrows(X - u, Y - v, X + u, Y + v, length = length, ...)
  )
}
Code 
res_symmetric <- symmetric(lon, lat)
res_antisymmetric <- antisymmetric(lon, lat)
res_both <- list(
  p = res_symmetric$p + res_antisymmetric$p,
  u = res_symmetric$u + res_antisymmetric$u,
  v = res_symmetric$v + res_antisymmetric$v,
  w = res_symmetric$w + res_antisymmetric$w)
res_hadley_symmetric <- hadley_symmetric(lat, lev) 
res_hadley_antisymmetric <- hadley_antisymmetric(lat, lev) 
res_hadley_both <- list(
  p = res_symmetric$p + res_antisymmetric$p,
  u = res_symmetric$u + res_antisymmetric$u,
  kai = res_symmetric$kai + res_antisymmetric$kai)
res <- res_symmetric
Code 
oldpar <- par(mfrow = c(3, 1), mar = c(4, 4, 2, 1), oma = c(3, 0, 0,0), mgp = c(2, 1, 0))

plot.new()
plot.window(c(-10, 15), c(-4, 4))

.filled.contour(lon, lat, res$w,
  levels = c(min(res$w), -0.1, 0.1, max(res$w)),
  col = c("lightblue", NA, "lightpink"))
clev_p <- pretty(res$p, 10)
if (any(clev_p %in% 0)) clev_p <- clev_p[-which(clev_p == 0)]
contour(lon, lat, res$p,
  levels = clev_p, add = TRUE)
ivec <- seq(1, nx, 5)
jvec <- seq(1, ny, 5)
drawvec(lon[ivec], lat[jvec], res$u[ivec, jvec], res$v[ivec, jvec],
  uvref = 1, angle = 30, length = 0.02)

axis(1)
axis(2)
title("(a) vector: (u, v), shade: w, contour: p", "", "longitude", "latitude")
box()

plot.new()
plot.window(c(-10, 15), c(0, pi))

wlk <- walker(lon, lev)
clev_kai <- pretty(wlk$kai, 10)
if (any(clev_kai %in% 0)) clev_kai <- clev_kai[-which(clev_kai == 0)]
lty_kai <- rep(1, length(clev_kai))
lty_kai[clev_kai < 0] <- 2
contour(lon, lev, wlk$kai, levels = clev_kai, lty = lty_kai,
        add = TRUE)

axis(1)
axis(2, c(0, pi / 2, pi), c("0", "D/2", "D"))
title("(b) meridionally averaged streamfunction", "", "longitude", "height")
box()

plot.new()
plot.window(c(-10, 15), c(-10, 0))

lines(lon, wlk$p)

axis(1)
axis(2)
title("(c) meridionally averaged pressure", "", "longitude", "pressure")
box()

Code 
par(oldpar)

もう一つの応答は \(n=1\) とした Equation 9.23 の解

\[\begin{aligned} u &= \frac{1}{2}(1+y^2)e^{-y^2/4}\\ p &= \frac{1}{2}(y^2-3)e^{-y^2/4}\\ v &= \left[F(x) + 4\varepsilon q_2(x)\right]ye^{-y^2/4}\\ w &= \dfrac{1}{2}\left[F(x) + q_2(x)(1+y^2)\right]e^{-y^2/4} \end{aligned}\]

である。 \(n=1\) の西進惑星波を表し、位相速度が\(1/3\), 減衰率が\(3\varepsilon\) である。

強制の西側（ \(x<-L\) ）では正味の西風となっている。

\[ \int_{-\infty}^\infty u\mathrm{d} y = -\sqrt{\pi}q_2(x) \]

強制域（ \(|x|<L\) ）では、南北風が熱源から極向きに吹いている。

\(\varepsilon\rightarrow 0\) の極限で

\[ (v,w)\rightarrow(y,1)F(x)e^{-1/4y^4} \]

となる。 Equation 9.10 及び Equation 9.11 から渦度方程式を作ると、 \(\varepsilon\rightarrow 0\) で

\[ y\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right) = 0 \]

となり収束（あるいは加熱）による渦度の生成が赤道付近からの惑星渦度の移流とバランスしている。

Code 
res <- res_hadley_symmetric
Code 
oldpar <- par(mfrow = c(3, 1), mar = c(4, 4, 2, 1), oma = c(3, 0, 0,0), mgp = c(2, 1, 0))

plot.new()
plot.window(c(-4, 4), c(0, pi))

clev_u <- pretty(res$u, 10)
if (any(clev_u %in% 0)) clev_u <- clev_u[-which(clev_u == 0)]
lty_u <- rep(1, length(clev_u))
lty_u[clev_u < 0] <- 2

contour(lat, lev, res$u, levels = clev_u, lty = lty_u, add = TRUE)

axis(1)
axis(2, c(0, pi / 2, pi), c("0", "D/2", "D"))
title("(a) Hadley circulation: zonally averaged zonal wind", "", "latitude", "height")
box()

plot.new()
plot.window(c(-4, 4), c(0, pi))

clev_kai <- pretty(res$kai, 10)
if (any(clev_kai %in% 0)) clev_kai <- clev_kai[-which(clev_kai == 0)]
lty_kai <- rep(1, length(clev_kai))
lty_kai[clev_kai < 0] <- 2

contour(lat, lev, res$kai, levels = clev_kai, lty = lty_kai, add = TRUE)

axis(1)
axis(2, c(0, pi / 2, pi), c("0", "D/2", "D"))
title("(b) Hadley circulation: zonally averaged streamfunction", "", "latitude", "height")
box()

plot.new()
plot.window(c(-4, 4), range(res$p))

lines(lat, res$p)

axis(1)
axis(2)
title("(c) Hadley circulation: zonally averaged pressure", "", "latitude", "pressure")
box()

赤道に対称な熱源対するハドレー循環。(a) \(u\)、(b) \(u,\,w\) に対する流線函数 、(c) \(p\)
Code 
par(oldpar)

9.2.2 赤道反対称な熱源に対する応答

赤道に反対称な熱源 \(Q_1 = F(x)\)

\[ Q(x,y) = Q_1D_1(y) = F(x)ye^{-y^2/4} \tag{9.25}\]

について微分方程式 Equation 9.21–Equation 9.23 を解く。

Equation 9.25 に対する応答も二つの部分から構成されている。

ひとつは \(n = 0\) の混合ロスビー重力波 Equation 9.22 で、熱源の外には伝播しない。 もうひとつは \(n = 2\) の惑星波である。 これらを合わせた解は次のように書ける。

\[\begin{aligned} p &= \frac{1}{2}q_3(x)y^3e^{-y^2/4}\\ u &= \frac{1}{2}q_3(x)(y^2-6)ye^{-y^2/4}\\ v &= \left[6\epsilon q_3(x)(y^2-1) + F(x)y^2)\right]e^{-y^2/4}\\ w &= \left[\frac{1}{2}q_3(x)y^2 + F(x)\right]ye^{-y^2/4} \end{aligned}\]

Code 
res <- res_antisymmetric
Code 
oldpar <- par(mar = c(4, 4, 2, 1), oma = c(3, 0, 0,0), mgp = c(2, 1, 0))

plot.new()
plot.window(c(-10, 15), c(-4, 4))

.filled.contour(lon, lat, res$w,
  levels = c(min(res$w), -0.1, 0.1, max(res$w)),
  col = c("lightblue", NA, "lightpink"))
clev_p <- pretty(res$p, 10)
if (any(clev_p %in% 0)) clev_p <- clev_p[-which(clev_p == 0)]
contour(lon, lat, res$p,
  levels = clev_p, add = TRUE)
ivec <- seq(1, nx, 5)
jvec <- seq(1, ny, 5)
drawvec(lon[ivec], lat[jvec], res$u[ivec, jvec], res$v[ivec, jvec],
  uvref = 1, angle = 30, length = 0.02)

axis(1)
axis(2)
title("vector: (u, v), shade: w, contour: p", "", "longitude", "latitude")
box()

赤道に反対称な熱源に対する応答。 矢印は \((u,\,v)\)、陰影は \(w\) 、等値線は \(p\)
Code 
par(oldpar)
Code 
res <- res_hadley_antisymmetric
Code 
oldpar <- par(mfrow = c(3, 1), mar = c(4, 4, 2, 1), oma = c(3, 0, 0,0), mgp = c(2, 1, 0))

plot.new()
plot.window(c(-4, 4), c(0, pi))

clev_u <- pretty(res$u, 10)
if (any(clev_u %in% 0)) clev_u <- clev_u[-which(clev_u == 0)]
lty_u <- rep(1, length(clev_u))
lty_u[clev_u < 0] <- 2

contour(lat, lev, res$u, levels = clev_u, lty = lty_u, add = TRUE)

axis(1)
axis(2, c(0, pi / 2, pi), c("0", "D/2", "D"))
title("(a) Hadley circulation: zonally averaged zonal wind", "", "latitude", "height")
box()

plot.new()
plot.window(c(-4, 4), c(0, pi))

clev_kai <- pretty(res$kai, 10)
if (any(clev_kai %in% 0)) clev_kai <- clev_kai[-which(clev_kai == 0)]
lty_kai <- rep(1, length(clev_kai))
lty_kai[clev_kai < 0] <- 2

contour(lat, lev, res$kai, levels = clev_kai, lty = lty_kai, add = TRUE)

axis(1)
axis(2, c(0, pi / 2, pi), c("0", "D/2", "D"))
title("(b) Hadley circulation: zonally averaged streamfunction", "", "latitude", "height")
box()

plot.new()
plot.window(c(-4, 4), range(res$p))

lines(lat, res$p)

axis(1)
axis(2)
title("(c) Hadley circulation: zonally averaged pressure", "", "latitude", "pressure")
box()

赤道に反対称な熱源に対するハドレー循環。(a) \(u\)、(b) \(u,\,w\) に対する流線函数 、(c) \(p\)
Code 
par(oldpar)
Note課題
  1. 赤道反対称な熱源に対する応答について解を求めよ。
  2. 赤道を挟んで発生した二つの熱帯低気圧の事例を探してみよ。
  3. 熱源に対する応答の観点からモンスーン循環について議論せよ。
  4. 東西一様な熱源に対する応答を求め、ITCZについて議論せよ。