準地衡風 2 層モデルを用いた中緯度大気循環の研究

{河合佑太 \atop 神戸大学地球惑星科学科}

平成 23 年 2 月 28 日

要旨

本研究は, 世界初の大気大循環の数値実験である「Phillips(1956) の数値実験」の再現を行う. Phillips は, 大気大循環をモデル化するために, 準地衡風 2 層モデルを用いて, 南北幅 10000 km の東西方向に開いたチャンネル内の流れを求めた. その際に, 系に簡単な熱源と冷源を与え, また地表面摩擦を考慮している. Phillips の数値実験は, 対流圏大気の特徴である中緯度の西風ジェット, 中緯度の地表の偏西風帯の存在や, 3 つの子午面循環などといった当時知られていた大気大循環の基本的な特徴を数値的に再現した. Phillips の数値実験の再現する本研究の意義は, 中緯度大気循環のより深い理解と数値予報の基本技術を得ることである.

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 1}

目次

第 1 章イントロダクション 4 1.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Phillips の数値実験の概要 . . . . . . . . . . . . . . . . . . . . . . . 5

1.2.1 数値予報に至るまでの歴史的背景 . . . . . . . . . . . . . . . 5 1.2.2 Phillips の数値実験の概要 . . . . . . . . . . . . . . . . . . . 7

第 2 章数値予報モデルの構築 : 支配方程式系 9

2.1 準地衡風予報方程式系 . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.1 準地衡風予報方程式系の導出 . . . . . . . . . . . . . . . . . 9 2.1.2 準地衡風予報方程式系による数値予報 . . . . . . . . . . . . 17 2.1.3 準地衡風予報方程式系を用いる理由 . . . . . . . . . . . . . . 18

2.2 Phillips の数値実験で使われた数値予報モデルの詳細 . . . . . . . . 19 2.2.1 2 層モデルの導入 . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.2 準地衡風予報方程式系の拡張 . . . . . . . . . . . . . . . . . 20 2.2.3 系の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2.4 非断熱加熱の表現 . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.5 数値実験で用いられた予報方程式系と境界条件 . . . . . . . 29

\overline{main.tex 2011/02/28(河合佑太)}

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 2}

2.3 系のエネルギー論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.3.1 エネルギー方程式 . . . . . . . . . . . . . . . . . . . . . . . . 32 2.3.2 系のエネルギー変換 . . . . . . . . . . . . . . . . . . . . . . 35

2.4 子午面循環 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

第 3 章数値予報モデルの構築 : 数値計算手法 40 3.1 予報方程式系の離散化 . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 数値計算手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3 エネルギーパラメータと子午面循環の数値的な評価法 . . . . . . . . 44

第 4 章数値実験の手順 47 4.1 数値実験の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.2 数値実験の手順 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.3 パラメータの設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

第 5 章数値実験の結果 50 5.1 数値実験 1 (渦なし状態における帯状流の発達) の結果 . . . . . . . . 50 5.2 数値実験 2 (渦あり状態における流れの発達) の結果 . . . . . . . . . 52

第 6 章数値実験の結果の考察 61 6.1 数値実験 1 (渦なし状態における帯状流の発達) の考察 . . . . . . . . 61 6.2 数値実験 2 (渦あり状態における流れの発達) の考察 . . . . . . . . . 62

第 7 章結論 70 \underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 3} 付録 A 鉛直座標変換に伴う微分演算子の変換則 73 付録 B 圧力座標系におけるプリミティブ方程式の導出 76 付録 C 準地衡風予報方程式系の導出における各式の変形 80 付録 D Phillips の数値予報モデルに関する数式の導出 83 付録 E Phillips の数値実験におけるエネルギー方程式の導出の詳細 85

第1章イントロダクション

1.1 はじめに

大気大循環の研究の始まりは 17 世紀までさかのぼるが, 数値実験による定量的な大気大循環の研究が始まったのは約半世紀前からである. 最初の大気大循環の数値実験は, Phillips(1956) によって行われた. この数値実験は, 数値予報が実用化され始めた 1950 年代の短期間の数値予報技術の応用である.

Phillips(1956) では, 簡単な大気大循環の構造をモデル化するために, 南北幅 10000 km の東西方向のチャンネル内の流れを求めた. この数値実験は歴史的に「Phillips の数値実験」と呼ばれている. その際に数値予報モデルとしては, 準地衡予報方程式系を適用した 2 層モデルを利用した. この数値予報モデルでは南北加熱差を生じさせるために, 非常に簡単で経験的な非断熱加熱や摩擦の効果を取り入れるが, 詳細にそのような物理過程^{*1} をモデル化しているわけではない. Phillips の実験の概要については, 後ほど述べることにする.

本研究では, Phillips の数値実験を再現する. Phillips の数値実験で用いられた準地衡風予報方程式系の 2 層モデル (さらにチャネルモデル) であれば, 数値予報モデルの構築は比較的簡単であり, かつ流体の傾圧性を表現することができる. それは, 自分の興味の一つである傾圧不安定波の発達の様子を, 2 層モデルという物理的に簡単化されたモデルにおいてではあるが, 数値シミュレーションすることを可能にする. Phillips の数値実験の再現を通して, 中緯度大気循環のより深い理解と数値予報の基本技術の習得が期待される.

Phillips(1956) を基にした本研究の目的は, 以下の通りである.

• Phillips(1956) の読書ノートを作成する. したがって, 本論文で使われる表現

*1大気大循環モデルの研究者は習慣的に, 移流項, コリオリ項, 気圧傾度項, 断熱項の計算と時間変化の項の計算を「力学過程」と呼ぶ. 一方で, 太陽放射や地球放射による加熱, 水の潜熱の放出・吸収, 乱流の熱輸送の収束, 分子粘性による熱拡散などの効果の計算を「物理過程」と呼ぶ.

\overline{main.tex 2011/02/28(河合佑太)}

の一部には, 原文の表現をそのまま和訳したものも含まれる.

• Phillips の数値実験で用いられた数値予報モデルを再現する中で, 数値予報モデル作成の基本技術を学ぶ.

– モデル構築の際には, 自作の数値計算ライブラリ「GMathLib」を使用し, その有効性を確認する. また, 実用的に使用する中で改良点を浮き彫りにする.

• 作成した数値予報モデルの中で表現される大気大循環の構造を解析し, 現実の地球におけるそれと比較する. より具体的には,

– 2 層モデルの中で発生する傾圧不安定波を研究する. また, そのような波が平均場に及ぼす影響を検証する.

– 系のエネルギーサイクルについて考察を行う.

1.2 Phillips の数値実験の概要

実際に再現した数値予報モデルの説明に入る前に, Phillips の数値実験の概要を述べたいと思う. また, Phillips の数値実験が数値予報の歴史を語る上で重要な出来事であることを説明するために, この節では最初に数値予報に至るまでの歴史的背景を簡単にまとめることにする^{*2}.

1.2.1 数値予報に至るまでの歴史的背景

人類の活動の上で, 日常の大気現象を把握することはとても重要であるために, 天気予報の歴史はとても長い. Reed(1977) は, 天気予報の早期の歴史において, 1643 年の Torricelli による気圧計の発明が科学としての気象学の始まりであると考られていると述べている. しかしながら, 初めて地上天気図が描かれたのは 1820 年である. しかも, それは約 50 年前の観測値を用いたものであった. 天気予報に対する天気図の可能性はすみやかに認知されたものの, 観測データを迅速に収集する手段が未発達であることが壁となった. この問題は 1845 年の電信の発明により解決され, 実況的な総観規模の天気図を作成できるようになった. そして, 天気予報や嵐への警戒情報を提供するために, 多くの国で国家事業もまた整備された.

*2この数値予報の歴史的背景は, 書籍 [2] と書籍 [12] のイントロダクションをまとめたものである.

\overline{main.tex 2011/02/28(河合佑太)}

1850-1920 年の期間における天気予報は, 経験的なものであった. 地上気圧の天気図は少し上層の大気の情報を知るために用いられたが, その手法はまったく経験則による. その背後では基礎的な理論研究が続けられていた. その中でも 1904 年の V. Bjerknes による認識は, 科学的な天気予報のその後の発展の道筋を決める上で, 重要かつ画期的な出来事であった. その認識は次のようになる.

• 天気予報は基本的には数学的には初期値問題である. さらに, 解かねばならない基礎方程式系は既に定式化されている.

• その方程式系は高次の非線形偏微分方程式であり, 一般的には解析解を得るのは不可能である. 大幅に簡単化された形式においてはそれは例外となり得るが, その結果は天気予報には直接的には使えない.

• 初期条件を決める上で観測値をそのまま利用することは不適切である.

彼自身はその後, 数値予報を実際に行うことはなかった. L. F. Richardson は, 第一次世界大戦の間, 手回し計算機を使って方程式系を数学的に解く方法を追い求めていた. その先駆的な努力は残念ながら失敗といえる結果に終わった^{*3}. 計算された圧力の時間変化は, 観測されたものより一つも大きなオーダーとなってしまったのである. 当時, 彼の行った研究の基礎的な重要性に気づく人は少なく, また, その失敗の要因を特定するために数カ月におよぶ計算を再試しようと考える人もいなかった. 一方で, ラジオゾンデの発明は 1930 年代の上層大気の観測値の広範囲の利用につながり, 急成長する航空産業の応用の増加とともに, 天気予報における経験則のよりいっそうの改良に導いた.

理論的な研究は, Rossby, Petterssen, J. Bjerknes, Charney, Eady, Eliassen, Fj ørtoft, Obukhov 他により続けられ, 実用的な問題に対する応用へと導いた. より重要なことは, これらの研究が純粋な経験論とは根本的に異なり, 物理学を基礎においている点である. しかし, 次の躍進への突破口を開いたのは 1940 年代後半の電子計算機の発明であった. 電子計算機の発展は, 数値予報に必要な膨大な演算処理を可能にした. そのような中で, Charney, Fjørtoft, von Neumann (1950) は, Richardson のモデルよりもずっと簡単化された等価順圧渦度モデルを用いた 500hPa 面の力学的な数値予報を初めて成功させた. その精度は, 主観的かつ経験的手法を使う熟年の予報官による予報と比べ得るものであった. こうして, 数値的手法による力学的な数値予報の時代が幕を開けたのである. 一般的には, 数値天気予報(numerical weather prediction(NWP)) と呼ばれ, 現代の天気予報において数値予報は主要な基盤となった. なお, 本論文ではこの数値天気予報のことを数値予報と略すことにする.

*3Richardson は, 1922 年にこの研究内容を書籍 [2] として出版している. その書籍の中で彼はこう述べている.

「いつか遠い将来, 天気の進行よりも早く計算を行うことが可能になり, 気象情報の入手のためのコストも安くなるだろう. しかし, それは夢である」

この言葉から, 数値予報はそれ以降「リチャードソンの夢」と呼ばれるようなった. さまざまな文献にも書かれていることではあるが, 当時の計算機資源や観測網などを考えれば, 結果に関わらず彼の業績は先駆的であり, その先見性は今の数値予報の発達の歴史が裏付けるところであるだろう.

\overline{main.tex 2011/02/28(河合佑太)}

その後, 数値予報が天気予報技術として業務化されていく中で, 数値予報モデルは少しづつ高度化・精密化されていく. もちろんその背景には, 計算機や観測手法の発展が密接に関わっている. 現在の予報精度に至るまでの, 数値予報モデルの発展の歴史はここでは触れないことにする^{*4}. ところで, Phillips の数値実験が行われた 1950 年代中期は数値予報業務化への準備段階にあった. 当時業務化された数値予報は, 数日先までの天気予報を目的としており, 数日程度であれば, 初期に与えたエネルギーで運動を維持したり低気圧を生み出せるので, 外部からエネルギーの補給のないモデルでも良かった. 一方で, Phillips の数値実験でも用いられた大気大循環モデルは, 当時の短期間予報用のモデルとは異なる特徴を持っている. 大気大循環モデルは, それ自身に熱源・冷源を備えており, モデルの中で循環を持続させる機構を持たせることで, 循環の統計的振る舞いを再現することを目的としているのが特徴である. なお, 現在では数値予報モデルの物理過程の精密化に伴って, 予報モデルと大気大循環モデルの違いはほとんどなくなっている.

1.2.2 Phillips の数値実験の概要

1956 年, Phillips は, 静止した大気を初期状態として 2 層準地衡風モデルを用いた (当時としては) 長期間の数値予報を行った. 支配方程式には, 摩擦および緯度に対して線形的な非断熱な効果が含まれている. この実験において, 原理的な経験的要素は, 加熱強度・鉛直安定度の値・摩擦による散逸の種類の 3 点である. 低緯度に対応する場所では熱源, 高緯度では冷源を与えるような非断熱加熱をモデルに設定することにより, ランダムな微小運動を初期値として時間積分を進める中で, 現実の大気に見られるような発達を伴う流れのパターンが数値的に実現された. その流れのパターンには, 大規模な擾乱の成長やそのような傾圧不安定波の作用による中緯度帯の西風ジェット, 地表の西風領域が含まれていた.

この数値実験は, 簡単な熱源や冷源を与えることにより, 対流圏大気の特徴である中緯度の西風ジェット, 中緯度の地表の偏西風帯の存在や, ３つの子午面循環 (ハドレー循環, フェレル循環, 極夜循環) などといったそれまでに知られていた大気大循環の基本的な特徴を数値的に再現した最初の実験である. この実験は, 大気大循環に関する最初の数値的実験とされている.

*4その歴史に興味をもたれた方は, 書籍 [12] 等を参考にしていただきたい.

\overline{main.tex 2011/02/28(河合佑太)}

Phillips はこの実験で約 1 カ月間の時間積分を行ったが, 非線形不安定によりエネルギーが少しずつ増大し, 平衡状態に至るまで時間積分を行うことができなかった. 要因は, 格子点法を用いる際に, エリアシングによる非線形不安定への配慮を行っていなかったからだとされる (Phillips, 1959)^{*5} .

*5非線形項の存在は, より高波数の波を表す項を生み出す. 例えば ψ_{1} = sin mx と ψ_{2} = sin nx で作られる非線形項を考えるとき,

ψ_{1} \frac{∂ψ_{2}}{∂x} = \frac{n}{2} \{sin[(m + n)x] + sin[(m − n)x]\}

となり, 非線形項はもとの波より高波数な波数 (m + n) の波を生み出すことになる. しかしながら, 格子点法を用いる場合, 格子幅で決定される最大波数までしか分解できず, それよりも高波数な波は, 実際よりも小さな波数の波として見なされてしまう. このことをエリアシングと呼ぶ.

また, このエリアシングを原因として生じる計算の不安定性は, 非線形不安定と呼ばれる.

この非線形不安定を防ぐためには, 支配方程式系の差分表現においてもエネルギーやエンストロフィが保存するようにする必要がある. なお, スペクトル法を用いる場合は, エリアシングが生じないように非線形項をスペクトル切断することによって, 非線形不安定を回避できる.

第2章{数値予報モデルの構築 : 支 \atop 配方程式系}

Phillips(1956) では, 準地衡風２層モデルが使われた. この章では, 最初にプリミティブ方程式系から準地衡風方程式系を導き, 後者による数値予報の方法を述べる. その後, Phillips の数値実験で用いられた準地衡風２層モデルの詳細についての説明を行う.

2.1 準地衡風予報方程式系

この節では, 始めに圧力座標系におけるプリミティブ方程式系から準地衡風予報方程式系を導出する. その後, 準地衡風予報方程式系を用いた数値予報について説明する. また, 準地衡風予報方程式と数値予報の歴史の関わりについてもこの節の最後で触れることにする.

2.1.1 準地衡風予報方程式系の導出

圧力座標系におけるプリミティブ方程式系

大気の運動を記述する基礎方程式系において, 静水圧平衡を仮定した方程式系は, プリミティブ方程式(primitive equation) と呼ばれる. 乾燥大気の運動を記述するプリミティブ方程式は, 圧力座標系において, 以下の運動方程式, 静水圧平衡の式, 連続の式, 理想気体の状態方程式, 熱力学的エネルギー保存の式によって構成される^{*1}.

\frac{dv}{dt} = −∇Φ − fk × v + F ,  (2.1.1)
\frac{∂Φ}{∂p} = −\frac{RT}{p} ,  (2.1.2)
\frac{∂u}{∂x} + \frac{∂v}{∂y} + \frac{∂ω}{∂p} = 0,  (2.1.3)
p = ρRT,  (2.1.4)
\frac{dln θ}{dt} = \frac{1}{c_{p}T} Q̇.  (2.1.5)

ここで, v = (u, v) は水平速度ベクトル, ω ≡ dp/dt は圧力座標系における鉛直速度, ρ は密度, p は圧力, T は温度, α(= ρ^{−1}) は比容, Φ(= gz) はジオポテンシャル, f(= 2Ω sin θ) はコリオリパラメータ, R は乾燥大気の気体定数, c_{v} は乾燥大気の定圧比熱, F は摩擦力, Q̇ は非断熱加熱による加熱率, k は鉛直上向きの単位ベクトルである. また,

\frac{d}{dt} = \frac{∂}{∂t} + u\frac{∂}{∂x} + v \frac{∂}{∂y} + ω \frac{∂}{∂p}

である.

総観スケールの擾乱に対するスケール解析

高層天気図に見られるような中緯度帯の大規模な波動擾乱に対して, 運動方程式 (2.1.1) の各項の大きさを見積もることから始めたいと思う.

ここでは, そのような波動擾乱がもつ特徴的なスケールとして次のようにとるとする.

• L = 水平方向の長さスケール ∼ 10^{6} m

• P = 鉛直方向の圧力の変動スケール ∼ 10^{5} Pa

• U = 水平速度スケール ∼ 10 m/s

• c = 擾乱の伝播速度 ∼ 10 m/s

• f = コリオリパラメータ ∼ 10^{−4}s^{−1} ここで, 擾乱の伝播速度 c が水平速度スケール U と同じ程度であるところが特徴である. これにより重力波や慣性重力波^{*2}のような中緯度帯の大規模な波動擾乱とは特徴の異なる波動を, 方程式系からフィルタリングすることができる.

*1プリミティブ方程式における幾何的な座標系から圧力座標系への鉛直座標変換の導出については, 付録 B に掲載してある.

\overline{main.tex 2011/02/28(河合佑太)}

運動方程式 (2.1.1) の中の各偏微分の大きさは, 上の特徴的なスケールを用いれば

\begin{aligned}∂/∂x & ∼ ∂/∂y ∼ 1/L ∼ 10^{−6} m^{−1} \\ ∂/∂p & ∼ 1/P ∼ 10^{−5} Pa^{−1} \\ u∂/∂x & ∼ v∂/∂y ∼ U/L ∼ 10^{−5} s^{−1} \\ ∂/∂t & ∼ c∂/∂x ∼ U/L ∼ 10^{−5} s^{−1} \end{aligned}  (2.1.6)

程度であると見積もられる. また, 連続の式 (2.1.3) から ω の大きさは最大でも,

ω ∼ P U/L ∼ 1 Pa s^{−1}  (2.1.7)

程度である.

よって, 運動方程式各項の大きさは次のように見積もられることが分かる.

\frac{∂v}{∂t} , (v · ∇)v, ω \frac{∂v}{∂p} ∼ \frac{U^{2}}{L} ∼ 10^{−4} m/s^{2}

k × fv ∼ fU ∼ 10^{−3} m/s^{−2}

上の各項の大きさの見積もりから, 水平方向の運動方程式において水平加速度の項とコリオリ加速度の項の比は, U/fL(≡ R_{o}) であり, 0.1 程度である. ここで R_{o} という無次元量は, ロスビー数(Rossby number) と呼ばれる重要なパラメータである. よって, 大規模でゆったりとした流れに対するスケーリングにおいて, 水平加速度の項はコリオリ加速度と比べて相対的に小さいことが分かる. コリオリ加速度の項は, 残りの項である水平気圧傾度力によってバランスされなければならないことになり, ロスビー数が 1 よりも十分小さい状況においては, 次の地衡風平衡が近似的に成り立つことになる.

fu_{g} = −\frac{∂Φ}{∂y}, fv_{g} = \frac{∂Φ}{∂x}  (2.1.8)

地衡風近似は, ロスビー数が大きくなるにつれて成り立たなくなる. ロスビー数が 0.1 であるような状況においては, 地衡風は実際の水平風と比べて 10% 程度の誤差を含んでいる.

*2重力波や慣性重力波の位相速度は 200∼350 m/s である.

\overline{main.tex 2011/02/28(河合佑太)}

準地衡風理論

大規模でゆったりした運動ではロスビー数は小さく, 第 1 近似としては水平風は地衡風であることを上で述べた. しかし, 地衡風平衡の式 (2.1.8) は時間微分を含まないため, この式からは速度場や圧力場の時間発展を予報することはできない. よって, プリミティブ方程式系の第 1 次近似が地衡風を与え, 第 2 次近似が地衡風の時間変化を, すなわち地衡風を時間的に変化させる大気の循環を表現するような予報方程式系を導く必要があるだろう. そのような方程式系を導くために, p 系のプリミティブ方程式系 (2.1.1)∼(2.1.5) を次のように書くとする^{*3}.

\frac{∂v}{∂t} + (v · ∇) v + ω \frac{∂v}{∂p} = −∇Φ − fk × v,  (2.1.9)
\frac{∂Φ}{∂p} = −α,  (2.1.10)
∇ · v + \frac{∂ω}{∂p} = 0,  (2.1.11)
\frac{∂}{∂t} (\frac{∂Φ}{∂p}) + v · ∇\frac{∂Φ}{∂p} + \frac{ω}{p} [\frac{∂}{∂p} (p\frac{∂Φ}{∂p} − κΦ)] = 0.  (2.1.12)

ここで, κ = R/c_{p} である. 上の方程式系では, 予報方程式の導出を簡単化するために摩擦力や非断熱加熱の項が省かれているが, これらの項の影響が無視できるという意味ではないことに注意されたい.

次に, 前に挙げた特徴的なスケールを使って変数を無次元化する. ω の特徴的なスケールには, (2.1.7) を用いることとする. このとき, 無次元変数を ∗ を使って表わせば, 有次元の変数は次のように無次元化される.

\begin{aligned}x^{∗} & = x/L, y^{∗} = y/L, p^{∗} = p/P, \\ t^{∗} & = t(U/L), v^{∗} = v/U, ω^{∗} = ω/(P U/L) \end{aligned}  (2.1.13)

これにより, ∗ がつく変数はその大きさが O(1) となる. ジオポテンシャル Φ は, 上のように無次元化する代わりに, 次のように表現する.

Φ(x, y, p, t) = Φ̄(p) + (f_{0}UL)Φ^{∗}(x^{∗}, y^{∗}, p^{∗}, t^{∗})  (2.1.14)

ここで, Φ̄ は, x, y, t についての Φ の平均であり, p のみに依存する関数である. そのため, x, y, p, t についての Φ の変化は右辺第 2 項において表現される. この第 2 項の大きさを f_{0}UL でスケーリングしたのは, ジオポテンシャル場の第 1 近似が地衡風を得るようなジオポテンシャル場であるようにしたいからである. この式において, f_{0} = 2Ω sin ϕ_{0} であり, 以下ベータ平面近似を用いることにする. すなわち, ある緯度 ϕ_{0} の惑星渦度勾配をつかって, 惑星渦度場を

*3(2.1.12) の導出については, 付録 C に掲載してある.

\overline{main.tex 2011/02/28(河合佑太)}

f = f_{0} + βy, β = (\frac{df}{dy} )_{ϕ=ϕ_{0}}  (2.1.15)

と線形的に近似するとする.

ここで, 無次元化された方程式系に含まれる 3 つの無次元パラメータを定義しておくことにしよう.

R_{o} = U/f_{0}L,  (2.1.16)
β^{∗} = L/(aR_{o} tan ϕ_{0}),  (2.1.17)
R_{i} = p\frac{d}{dp} [p\frac{dΦ̄}{dp} − κΦ̄] \frac{1}{U^{2}} = −\frac{p^{2}ᾱ}{U^{2}} \frac{∂ln θ}{∂p}.  (2.1.18)

R_{o} は既に出てきたロスビー数である. (2.1.17) の中の a は, 地球の半径である. R_{i} はリチャードソン数^{*4}に対応しており, 大気の静的安定度を表している. その値は地表近くで約 80, 対流圏界面で約 65, 成層圏で約 200 くらいである. 今回は対流圏内の現象を扱うことにして, 以下 R_{i} を一定だとし, その大きさは 10^{2} ぐらいとみなすとする. また, 上で定義した β^{∗} を用いれば惑星渦度 (2.1.15) は,

\begin{aligned}f & = f_{0} (1 + \frac{β}{f_{0}} Ly^{∗}) = f_{0} (1 + \frac{L}{a tan ϕ_{0} R_{o}} R_{o}y^{∗}) \\ & = f_{0} (1 + R_{o}β^{∗}y^{∗}) \end{aligned}  (2.1.19)

と書き換えられる. 上で定めた特徴的なスケールを用いれば, 無次元パラメータの値はそれぞれ R_{o} ∼ 10^{−1}, β^{∗} ∼ 1, R_{i} ∼ 10^{2} となる.

R_{i} = g \frac{∂ln θ}{∂z} (\frac{∂u}{∂z})^{−2}

*4リチャードソン数は, のように定義される.

(2.1.18) の中辺から最右辺への式変形の詳細は, 付録 C に掲載している.

\overline{main.tex 2011/02/28(河合佑太)}

準備が整ったので, (2.1.9), (2.1.10), (2.1.11) を無次元変数で書き直せば^{*5},

R_{o}(\frac{∂u^{∗}}{∂t^{∗}} + u^{∗} \frac{∂u^{∗}}{∂x^{∗}} + v^{∗} \frac{∂u^{∗}}{∂y^{∗}} + ω^{∗} \frac{∂u^{∗}}{∂p^{∗}} )
+ \frac{∂Φ^{∗}}{∂x^{∗}} − (1 + R_{o}β^{∗}y^{∗})v^{∗} = 0,  (2.1.20)
R_{o}(\frac{∂v^{∗}}{∂t^{∗}} + u^{∗} \frac{∂v^{∗}}{∂x^{∗}} + v^{∗} \frac{∂v^{∗}}{∂y^{∗}} + ω^{∗} \frac{∂v^{∗}}{∂p^{∗}} )
+ \frac{∂Φ^{∗}}{∂y^{∗}} + (1 + R_{o}β^{∗}y^{∗})u^{∗} = 0,  (2.1.21)
\frac{∂ω^{∗}}{∂p^{∗}} = − (\frac{∂u^{∗}}{∂x^{∗}} + \frac{∂v^{∗}}{∂y^{∗}} ) ,  (2.1.22)
R_{o}(\frac{∂}{∂t^{∗}} + u^{∗} \frac{∂}{∂x^{∗}} + v^{∗} \frac{∂}{∂y^{∗}} ) (\frac{∂Φ^{∗}}{∂p^{∗}} )
+ \frac{ω^{∗}}{(p^{∗})^{2}} [R_{i}R_{o}^{2} + R_{o}p^{∗} \frac{∂}{∂p^{∗}}(p^{∗} \frac{∂Φ^{∗}}{∂p^{∗}} − κΦ^{∗})] = 0. (2.1.23)

この方程式において, u^{∗}, v^{∗}, ∂/∂x^{∗}, ∂/∂y^{∗}, ∂/∂t^{∗}, β^{∗}, R_{i}R_{o}^{2} は O(1) である. 一方で, ω^{∗}, Φ^{∗} の大きさに関しては, 最大が 1 程度であるというだけで, どれくらいの大きさを持つ量なのかはまだ分からない.

次に, 全ての従属変数を R_{o} のべき級数に展開する. 例えば, u^{∗} については,

u^{∗} = u_{0}^{∗} + R_{o}u_{1}^{∗} + R_{o}^{2}u_{2}^{∗} + ⋯  (2.1.24)

となる. ここで, u_{0}^{∗}, u_{1}^{∗}, u_{2}^{∗} などはそれぞれ O(1) の大きさを持つ. 他の従属変数に対しても同様のべき級数展開を行い, それらを (2.1.20) ∼ (2.1.23) へ代入する. そして, R_{o} についての各次の項を別々に集める.

R_{o} の 0 次の項

v_{0}^{∗} = \frac{∂Φ_{0}^{∗}}{∂x^{∗}} ,  (2.1.25)
u_{0}^{∗} = \frac{∂Φ_{0}^{∗}}{∂y^{∗}} ,  (2.1.26)
\frac{∂ω_{0}^{∗}}{∂p^{∗}} = − (\frac{∂u_{0}^{∗}}{∂x^{∗}} + \frac{∂v_{0}^{∗}}{∂y^{∗}} ) ,  (2.1.27)
ω_{0}^{∗} = 0.  (2.1.28)

(2.1.25), (2.1.26) は, 地衡風平衡の式 (2.1.8) とほとんど同じである. しかし, 地衡風平衡の式に現われる緯度に依存するコリオリパラメータが, 上の式においては定数 (ある規準の緯度 ϕ_{0} での値) であることに注意が必要である. したがって, R_{0} の 0 次の水平速度場は非発散となる. このとき (2.1.27) の右辺はゼロとなり, ∂ω_{0}^{∗}/∂p^{∗} = 0 が得られる. また, 大気上端 (p^{∗} → 0) において ω_{0}^{∗} → 0 となることも考慮すれば, ω_{0}^{∗} = 0 であることが示される. この結果は, (2.1.28) と矛盾しない. 逆を言えば, (2.1.25), (2.1.26) の地衡風の式においてコリオリパラメータを一定に取らなかった場合, (2.1.25)∼(2.1.28) の間で相互矛盾が生じてしまう.

*5(2.1.23) の導出は, 付録 C を参照.

\overline{main.tex 2011/02/28(河合佑太)}

R_{o} の 1 次の項

(\frac{∂}{∂t^{∗}} + u_{0}^{∗} \frac{∂}{∂x^{∗}} + v_{0}^{∗} \frac{∂}{∂y^{∗}} ) u_{0}^{∗} + \frac{∂Φ_{1}^{∗}}{∂x^{∗}} − β^{∗}y^{∗}v_{0}^{∗} − v_{1}^{∗} = 0,  (2.1.29)
(\frac{∂}{∂t^{∗}} + u_{0}^{∗} \frac{∂}{∂x^{∗}} + v_{0}^{∗} \frac{∂}{∂y^{∗}} ) v_{0}^{∗} + \frac{∂Φ_{1}^{∗}}{∂y^{∗}} − β^{∗}y^{∗}v_{0}^{∗} + u_{1}^{∗} = 0,  (2.1.30)
\frac{∂ω_{1}^{∗}}{∂p^{∗}} = −(\frac{∂u_{1}^{∗}}{∂x^{∗}} + \frac{∂v_{1}^{∗}}{∂y^{∗}} ),  (2.1.31)
(\frac{∂}{∂t^{∗}} + u_{0}^{∗} \frac{∂}{∂x^{∗}} + v_{0}^{∗} \frac{∂}{∂y^{∗}} ) (\frac{∂Φ_{0}^{∗}}{∂p^{∗}} ) + ω_{1}^{∗}σ(p^{∗}) = 0.  (2.1.32)

ここで,

σ(p^{∗}) = R_{i}R_{o}^{2}/(p^{∗})^{2}

である. 今, R_{i}R_{o}^{2}, p^{∗} ∼ O(1) であるので, σ は O(1) である.

続いて, R_{o} について 0 次の速度から求まる渦度

ζ_{0}^{∗} = \frac{∂v_{0}^{∗}}{∂x^{∗}} − \frac{∂u_{0}^{∗}}{∂y^{∗}} = ∇^{2}Φ_{0}^{∗}  (2.1.33)

に対する渦度方程式を導く. R_{o} の 1 次についての方程式系 (2.1.29)∼(2.1.32) において, u_{1}^{∗}, v_{1}^{∗}, Φ_{1}^{∗} を消去するために, (2.1.30) に ∂/∂x, (2.1.29) に −∂/∂y を作用させて, 両辺を加え合わせる. そして, 添え字 1 の風速場の発散に (2.1.31) を代入し, 地衡風の渦度 (2.1.33) を導入すれば,

\frac{∂ζ_{0}^{∗}}{∂t^{∗}} + u_{0}^{∗} \frac{∂ζ_{0}^{∗}}{∂x^{∗}} + v_{0}^{∗} (β^{∗} + \frac{∂ζ_{0}^{∗}}{∂y^{∗}} ) − \frac{∂ω_{1}^{∗}}{∂p^{∗}} = 0  (2.1.34)

を得る^{*6}.

(2.1.34) と (2.1.32) を有次元の変数で書き戻せば,

\frac{∂ζ_{g}}{∂t} = −v_{g} · ∇(ζ_{g} + f) + f_{0} \frac{∂ω}{∂p} ,  (2.1.35)
\frac{∂}{∂t} (−\frac{∂Φ′}{∂p} ) = −v_{g} · ∇ (−\frac{∂Φ′}{∂p} ) + σω.  (2.1.36)

*6(2.1.34) の導出は, 付録 C を参照.

\overline{main.tex 2011/02/28(河合佑太)}

ただし,

ζ_{g} = ∇^{2}Φ′/f_{0}, v_{g} = k × ∇Φ′/f_{0}
σ = −\frac{ᾱ}{θ̄} \frac{dθ̄}{dp} = \frac{1}{p} \frac{d}{dp} [p\frac{dΦ̄}{dp} − κΦ̄]  (2.1.37)

であり, Φ′ は, 等圧力面上のジオポテンシャルの平均値からのずれを表す. ふつうは, σ > 0 である.

(2.1.35) と (2.1.36) から成る方程式系について考えると, 未知変数は Φ′ と ω の二つ, 方程式は二本あるので方程式系は閉じている. この方程式系は, 準地衡風予報方程式系(quasi-geostrophic forecast equations) と呼ばれる. この 2 式において渦度, 渦度および層厚 (温度) の移流には, 実際の風でなく地衡風が用いられるが, 水平発散 (∇ · v = −∂ω/∂p) の計算に関してはそうなっていない. もし水平発散を地衡風で計算したならばその値は必ずゼロになってしまい, 連続の式からモデルの中で鉛直流は生じないことになる. その場合には絶対渦度 (ζ_{g} + f) がラグランジュ保存量となり, 流体要素がもつ絶対渦度は時間的に常に保存することになるが, 水平発散項の変化 (ゼロとは限らない) により流体要素の絶対渦度は変化する^{*7}.

次に, 準地衡風方程式系におけるポテンシャル渦度方程式を導こう. (2.1.35) と (2.1.36) から ω を消去すれば, 唯一つの従属変数 Φ′ に対する方程式

\frac{∂q}{∂t} + v_{g} · ∇q = 0  (2.1.38)

が得られる. ここで,

q = f + ζ_{g} + f_{0} \frac{∂}{∂p} (\frac{1}{σ} \frac{∂Φ′}{∂p} )  (2.1.39)

であり, 準地衡風方程式系におけるポテンシャル渦度である. (2.1.38) は, 非断熱効果および摩擦の効果がないときには, ポテンシャル渦度 q が保存することを示している.

また, (2.1.35) と (2.1.36) から時間微分の項 ∂Φ′/∂t を消去すれば,

\begin{aligned}(∇^{2} + \frac{f_{0}^{2}}{σ} \frac{∂^{2}}{∂p^{2}} ) ω = & \frac{f_{0}}{ω} \frac{∂}{∂p} [v_{g} · ∇ (\frac{1}{f_{0}} ∇^{2}Φ′ + f)] \\ & + \frac{1}{ω} ∇^{2} [v_{g} · ∇ (−\frac{∂Φ′}{∂p} )] \end{aligned}  (2.1.40)

*7プリミティブ方程式から導かれる渦度方程式には, 発散項の他に立ち上がり項も含まれ, それが流体要素の絶対渦度の変化に寄与することが分かる. 立ち上がり項は, 鉛直流の水平分布が一様でない場合に, 渦管が起き上がって水平渦度が鉛直渦度に変換される効果を表現する. しかしながら, 準地衡風世界の低次のオーダーに対する渦度方程式には, 立ち上がり項は現れない.

\overline{main.tex 2011/02/28(河合佑太)}

が得られる. この式は, オメガ方程式(omega equation) と呼ばれる. この式により, ある時刻のポテンシャル渦度の擾乱場 Φ′ が与えられたならば, その時刻の ω の場を診断的に求めることができる.

2.1.2 準地衡風予報方程式系による数値予報

準地衡風予報方程式系を用いて数値予報を行う方法には, 次の二種類がある.

方法 1

• 時刻 t = t_{0} の Φ′ の空間分布が格子点上のデータとして与えられる. このとき, この時刻の地衡風 v_{g} が計算される.

• Φ′, v_{g} からオメガ方程式 (2.1.40) の右辺が計算され, 適切な境界条件のもとでこの方程式を解けば, 時刻 t = t_{0} における ω の場が求められる.

• Φ′, ω から, 渦度方程式 (2.1.35) の右辺を計算し, 時刻 t_{0} における ∂ζ_{g}/∂t を求める.

• 時刻 t_{0} での ζ_{g} の時間変化率が分かったので, 時刻 t_{0} + δt の ζ_{g} が計算される.

• ∇^{2}Φ′ = f_{0}ζ_{g} を解けば, 時刻 t_{0} + δt における Φ′ の場が計算される.

• 求めた時刻 t_{0} + δt の値を新たな初期値として, 最初の手順へと戻る. 目標の時刻までこのルーチンを繰り返す.

方法 2

• 時刻 t = t_{0} の Φ′ の空間分布が格子点上のデータとして与えられる. このとき, この時刻の地衡風 v_{g}, その相対渦度 ζ_{g} が計算される.

• 与えた Φ′ と求めた ζ_{g} から, (2.1.39) より時刻 t = t_{0} のポテンシャル渦度 q が計算される.

• 求めた v_{g}, q を用いて (2.1.38) から, 時刻 t = t_{0} における ∂q/∂t が求まる.

• 時刻 t_{0} での q の時間変化率が分かったので, 時刻 t_{0} + δt の q が計算される.

• 時刻 t_{0} + δt の q を用いて (2.1.39) を逆に解けば, 時刻 t_{0} + δt の Φ′ が計算される.

• 求めた時刻 t_{0} + δt の値を新たな初期値として, 最初の手順へと戻る. 目標の時刻までこのルーチンを繰り返す.

なお, Phillips の数値実験では方法 2 が使われた.

2.1.3 準地衡風予報方程式系を用いる理由

プリミティブ方程式系をそのまま予報方程式系として用いるのではなく, 1950 年代から 1960 年代にかけて準地衡風予報方程式系が用いられた理由について, まず説明したいと思う.

前節までに挙げたプリミティブ方程式系は, 数日周期で天気を変化させるような低周波の波動から, 重力波や慣性重力波といった高周波の波動まで解として含んでいる. 数日周期の天気変化を数値予報したい際に, 離散化手法によっては後者の高周波の波動が問題となる. 高速な音波や重力波といった高周波の波動はたいてい振幅が小さく, 総観的な天気変化に対して本質でないことが多い. しかし, 数値計算手法の選び方によっては, この高周波の波に対して偽の増幅を予報してしまうことがある. その場合, 数日周期の天気の変化に関連するような擾乱は, 数値計算によって発生したこのノイズによってぼやけてしまう.

Richardson の失敗の要因は, そのようなところにあったとされる. すなわち, 数日周期の天気変化を予報するために, Richardson はプリミティブ方程式系を支配方程式として選んだのであるが, その際適用したスキームは, 解として含まれる高周波の波^{*8}に対して偽の増幅を予報してしまうものであった. 結果, 本来予報したかった低周波の擾乱はノイズの中に埋もれてしまい, Rhichardson の先駆的な仕事は失敗に終わってしまった.

Charney et al. (1948) の等価順圧モデルでは, この問題を解決するために力学方程式に地衡風近似・静水圧近似を系統的に導入した. これにより支配方程式系は簡単化され, 数値計算上問題となる音波や重力波は方程式の解から除去される. かくして, Charney たちの数値予報は成功したのである. なお, このフィルター近似により得た予報方程式系は, 準地衡風予報方程式系と本質的には等価である. そのような歴史的経緯の中で, 数値予報の歴史のごく初期には準地衡風予報方程式系が用いられた^{*9}.

*8Richardson は, 初期値として誤差を含む観測値をそのまま用いた. 初期値として与えた風速のうち非地衡風成分は, 高周波の慣性重力波として伝播していく.

\overline{main.tex 2011/02/28(河合佑太)}

2.2 Phillipsの数値実験で使われた数値予報モデルの詳細

この節では, まず 2 層モデルについて簡潔に説明することにする. また, 前節で導いた摩擦や非断熱効果を無視した準地衡風予報方程式系を拡張し, Phillips の数値実験で用いられた摩擦や非断熱効果を考慮した予報方程式系を導く. その後, Phillips の数値実験で設定された系についての説明を行うことにする.

2.2.1 2 層モデルの導入

{図 2.2.1: 2 層モデルにおける鉛直方向の \atop 従属変数の配置 } 2 層モデルは大気を 2 分割して, 支配方程式系の鉛直方向の離散化を行う. その結果, 連続成層における傾圧不安定波の議論に必要とされるような複雑な数学的取り扱いを簡単化することができる. 具体的には, 図 2.2.1 のように連続に成層した大気を 0, 2, 4 と番号を付けた等圧面で区切られた２つの層で分割する. これらの等圧面はそれぞれ気圧が 0, 500, 1000 hPa であるとする. そして, レベル 1 と 3 では水平速度 (流線関数の摂動成分), レベル 0, 2, 4 では鉛直速度を与え, レベル

1 と 3 において準地衡風渦度方程式を, レベル 2 において熱力学的エネルギー保存の式を適用する. レベル 2 や 4 で風速場の情報が必要になった場合は, レベル 1 と 3 における風速場の情報を使って線形補間する.

対流圏大気の運動の解析を簡単化するために, 2 層モデルの支配方程式の定式化は圧力座標系で行われる. 大気上端と下端の境界条件として, 鉛直速度は

*9しかしながら, 準地衡風近似は低緯度では成立しない上に, フィルターによって除去されたスケールの運動を予報できないために, 数値予報の歴史が進むにつれて, 再びプリミティブ方程式が予報方程式系として使われるようになった. プリミティブ方程式では静水圧近似を仮定しているが, これを仮定しない非静力学モデルもある.

\overline{main.tex 2011/02/28(河合佑太)}

ω = 0 at p = 0 and p_{4}  (2.2.41)

のように取られる. 成層圏の大きな静的安定度は鉛直運動を抑制するので, 大気上端の鉛直速度をゼロとすることは自然である. 一方で, 大気下端の鉛直速度をゼロと仮定することは数学的な取り扱いを簡単にするためである.

2.2.2 準地衡風予報方程式系の拡張

Phillips の数値実験では, 地面による摩擦, 渦拡散による散逸, 非断熱効果も数値予報モデルの中で表現している. よって, 前節では考慮しなかった摩擦および非断熱効果を準地衡風予報方程式系に組み入れることになる. その場合, 準地衡風近似を施した運動方程式, 連続の式, 熱力学的エネルギー保存の式は,

\frac{∂u}{∂t} + v · ∇u − fv = −\frac{∂Φ}{∂x} + A_{v}∇^{2}u + g \frac{∂τ_{zx}}{∂p} ,  (2.2.42)
\frac{∂u}{∂t} + v · ∇v + fu = −\frac{∂Φ}{∂y} + A_{v}∇^{2}v + g \frac{∂τ_{zy}}{∂p} ,  (2.2.43)
∇ · v + \frac{∂ω}{∂p} = 0,  (2.2.44)
\frac{dln θ}{dt} = \frac{1}{c_{p}T} \frac{dQ}{dt}  (2.2.45)

となる. ここで, v は p 系における水平速度ベクトル, A_{v} は水平方向の動渦粘性係数, τ_{zx}, τ_{zy} は水平面に沿って働く摩擦応力の x, y 成分である. また, ベータ平面近似

f = f_{0} + βy  (2.2.46)

を適用している.

この方程式系における各変数は, 前節と同じく R_{o} の 1 次までの精度を持っている. すなわち, 前節の表記を用いるならば

u = U(u_{0}^{∗} + R_{o}u_{1}^{∗}), v = U(v_{0}^{∗} + R_{o}v_{1}^{∗}), ω = \frac{P U}{L} R_{o}ω_{1}^{∗}, ω_{0}^{∗} = 0

Φ = Φ̄(p) + f_{0}UL(Φ_{0}^{∗} + R_{o}Φ_{1}^{∗}).

である. 今, R_{o} の 0 次の風速は地衡風平衡

f_{0}u_{g} = −\frac{∂Φ}{∂y}, f_{0}v_{g} = −\frac{∂Φ}{∂x}

にあることに注意されたい. また, 前節の議論から鉛直速度が R_{o} の 1 次のオーダーであるために, 運動方程式中の鉛直方向の移流項は, 水平方向の移流項に対して無視されている.

ここで, 速度 v を地衡風成分 v_{g} と非地衡風成分 v_{a}(R_{o} の 1 次の風速に対応) に分割して, 準地衡風の渦度方程式を再度導こう^{*10} .

v = v_{g} + v_{a}  (2.2.47)

を, (2.2.42)∼(2.2.45) に代入し, 地衡風平衡で釣り合っている項を消去すれば,

\frac{d_{g} u_{g}}{dt} − βyv_{g} − f_{0}v_{a} = A_{v}∇^{2}u + g \frac{∂τ_{zx}}{∂p} ,  (2.2.48)
\frac{d_{g} v_{g}}{dt} + βyu_{g} + f_{0}u_{a} = A_{v}∇^{2}v + g \frac{∂τ_{zy}}{∂p} ,  (2.2.49)
∇ · v_{a} + \frac{∂ω}{∂p} = 0,  (2.2.50)
\frac{dln θ}{dt} = \frac{1}{c_{p}T} \frac{dQ}{dt}  (2.2.51)

を得る^{*11}. ただし,

\frac{d_{g}}{dt} = \frac{∂}{∂t} + u_{g} \frac{∂}{∂x} + v_{g} \frac{∂}{∂y}  (2.2.52)

である. 前節と同様にして (2.2.49) に ∂/∂x, (2.2.48) に −∂/∂y を作用させて加え合わせる. そして, 非地衡風成分の水平発散を (2.2.50) で消去すれば,

\frac{d_{g}}{dt}(f + ζ_{g}) = f_{0} \frac{∂ω}{∂p} + A_{v}∇^{2}ζ_{g} + g \frac{∂}{∂p} (\frac{∂τ_{zy}}{∂x} − \frac{∂τ_{zx}}{∂y} )  (2.2.53)

を得る. ここまでは, 地衡風成分であることを添え字 g を使って強調してきたが以後この添え字を省略することにする. したがって, 今後予報方程式系に出てくる水平速度の変数は, 基本的には地衡風の風速を意味すると思ってもよい.

次は, この渦度方程式と熱力学的エネルギーの式を鉛直方向に離散化することになるが, その前に, 予報変数によって表現しなけばならない応力について説明する.

*10ここでなされている準地衡風に対する渦度方程式の導出には, 前節とは異なり摩擦項や非断熱加熱の項が含まれている. しかしながら, 前節の摩擦および非断熱加熱の効果を無視した準地衡風理論が踏まえられており, その理論の本質は大きく変わらないことを強調しておく.

*11運動方程式において, 非地衡風成分のコリオリ加速度の項の係数としてかかっているのは, コリオリパラメータ f でなく f_{0} となっている. これは, 上の運動方程式が非地衡風成分のオーダー (R_{o} の 1 次に対応するオーダー) に関するつり合いの式だからである. つまり, f = f_{0} + βy のうち惑星渦度勾配に関する項と非地衡風成分の積は R_{O} の 2 次に対応する項であり, 上式のつり合いの中には現れない.

\overline{main.tex 2011/02/28(河合佑太)}

応力の表現

次に接線応力の表現について考える. p = 0 の圧力面に働く接線応力 τ_{zx}, τ_{zy} は必然的にゼロであるが, p = 500 hPa における値は実際のところ分からない. よって, 500 hPa 面に働く接線応力は 1000 hPa のものに比べてずっと小さいと仮定し, 500 hPa における接線応力をゼロにとることにする.

一般的に, 地表における接線応力は, 地上風 C(x 方向成分を u_{s}, y 方向成分を v_{s} とする) と逆向きで, その大きさが観測計のある高度における地上風速の 2 乗に比例するように取られる. すなわち,

|τ_{4}| = κ_{c}ρC^{2}.  (2.2.54)

あるいは成分で書けば,

(τ_{zx})_{4} = −κ_{c}ρCu_{s}, (τ_{zy})_{4} = −κ_{c}ρCv_{s}.  (2.2.55)

ここで, κ_{c} は無次元定数であり約 0.003 程の値である. また, 添え字の 4 はレベル 4 での応力であることを示している. 以下, 各レベルの変数であることを添え字の数字で表現することにする. (観測計のある高度での) 地上風速が地衡風の約 70 % だと仮定すれば, 鉛直拡散の項は次のように書くことができるだろう^{*12}.

\frac{g}{p_{2}} (\frac{∂τ_{zy}}{∂x} − \frac{∂τ_{zx}}{∂y} )_{4} = −kζ_{4}.  (2.2.58)

ただし, 定数 k は,

k = 0.98 κ_{c} \frac{gC̄_{g}}{RT_{4}}  (2.2.59)

であり, C̄_{g} は地表における地衡風の特徴的な値である. また, ζ_{4} は地表における地衡風の渦度である. C̄_{g} が約 11 m/s 程のとき, k の値は 4 × 10^{−6} s^{−1} となる.

支配方程式系の鉛直離散化

応力の取り扱い方が定まったので, 渦度方程式 (2.2.53) と熱力学的エネルギー保存の式 (2.2.45) の鉛直離散化を行う. 「2 層モデルの導入」の小節で述べたように,

\frac{gρ_{4}}{p_{2}} = 2g \frac{ρ_{4}}{p_{4}} = \frac{2g}{RT_{4}}  (2.2.57)

*12(2.2.58) を導出する.

今, 地上風速 u_{s}, v_{s} はレベル 4 の風速の約 70 %と考えているので, (2.2.55) はレベル 4 の変数を使って,

(τ_{zx})_{4} = −0.49 κ_{c}ρCu_{4}, (τ_{zy})_{4} = −0.49 κ_{c}ρCv_{4}  (2.2.56)

と書き直せる. さらに, であることも用いる. 最後に, C を C̄_{g} と書き換え, レベル 4 の渦度を表わす変数 ζ_{4} を導入すれば, (2.2.58) を得る.

\overline{main.tex 2011/02/28(河合佑太)}

渦度方程式はレベル 1 および 3, 熱力学的エネルギー保存の式はレベル 2 にて適用される.

渦度方程式に含まれる ω の p 微分は中心差分を使って,

(\frac{∂ω}{∂p} )_{1} = \frac{ω_{2}}{p_{2}} , (\frac{∂ω}{∂p} )_{3} = −\frac{ω_{2}}{p_{2}}  (2.2.60)

と近似する. ここで, (2.2.41) と p_{2} − p_{0} = p_{4} − p_{2} = p_{2} であることを用いた. これにより, レベル 1 とレベル 3 で適用される渦度方程式はそれぞれ,

(\frac{∂}{∂t} + v_{1} · ∇) (f + ζ_{1}) − f_{0} \frac{ω_{2}}{p_{2}} = A_{v}∇^{2}ζ_{1}  (2.2.61)
(\frac{∂}{∂t} + v_{3} · ∇) (f + ζ_{3}) + f_{0} \frac{ω_{2}}{p_{2}} = A_{v}∇^{2}ζ_{3} − kζ_{4}  (2.2.62)

となる.

続いて熱力学的エネルギー保存の式の鉛直離散化を行う. 温位の定義式と静水圧平衡の式から,

ln θ = ln (\frac{∂Φ}{∂p}) + ln F (p)  (2.2.63)

を得る. ただし,

F (p) = −\frac{p}{R} (\frac{p_{4}}{p} )^{κ}

である. レベル 2 の (2.2.63) を表すために, Φ の p 微分を中心差分で表せば,

(ln θ)_{2} = ln \frac{Φ_{3}(x, y, t) − Φ_{1}(x, y, t)}{p_{2}} + ln F (p_{2})

となる. よって, 局所時間微分と水平移流項に関して,

(\frac{∂}{∂t} + v_{2} · ∇) (ln θ)_{2} = \frac{1}{Φ_{3} − Φ_{1}} (\frac{∂}{∂t} + v_{2} · ∇) (Φ_{3} − Φ_{1})  (2.2.64)

を得る. また, 鉛直移流項については,

(ω \frac{∂}{∂p} ln θ)_{2} = \frac{ω_{2}}{θ_{2}} \frac{θ_{3} − θ_{1}}{p_{2}}  (2.2.65)

である. レベル 2 の温度は静水圧平衡の式を使って次のように得られる.

T_{2} = −\frac{p_{2}}{R} (\frac{∂Φ}{∂p})_{2} = −\frac{Φ_{3} − Φ_{1}}{R} .  (2.2.66)

結果 (2.2.64), (2.2.65), (2.2.66) を用いて, レベル 2 で適用される熱力学的エネルギー保存の式は次のように鉛直離散化される.

\frac{ω_{2}}{p_{2}} [\frac{θ_{1} − θ_{3}}{θ_{2}} ] (Φ_{1} − Φ_{3}) = (\frac{∂}{∂t} + v_{2} · ∇) (Φ_{1} − Φ_{3}) − \frac{R}{c_{p}} \frac{dQ}{dt}  (2.2.67)

2 層モデルでは, 鉛直方向の静的安定性を表現するパラメータ (θ_{1} − θ_{3})/θ_{2} の変化を予報できない. ただし, (2.2.67) の左辺において Φ_{1} − Φ_{3} の変化は ω_{2} の変化に比べてずっと小さいため, ここでは左辺の Φ_{1} − Φ_{3} を定数として取り扱うことにする. また, (2.2.67) の他の項や他の方程式にあらわれる Φ_{1} と Φ_{3} に注目してみると, いずれも時間微分や水平空間微分によって作用されるので, 結局 Φ_{1} と Φ_{3} の静水圧場からのずれのみを予報することになる.

2.2.3 系の設定

前の小節において, 2 層モデルにおける渦度方程式と熱力学的エネルギー保存の式の鉛直離散化が完了した. その際には, レベル 4 で働く応力を風速場から与えることで摩擦項の具体的な数学的表現が決められたが, 非断熱加熱による効果の数学的表現は依然として決定されていない. この表現を具体的に決定する前に, Phillips の数値実験で設定された系について先に説明することにする. その後, 非断熱加熱の数学的表現を与えたいと思う.

{図 2.2.2: 設定された系と直交座標系の \atop 導入} Phillips の数値実験は, 南北端 (y = ±W) は「壁」となっていて, 東西端 (x = 0, L) に対しては周期境界条件が満たされるような矩形領域で行われた. また, 水平方向に対しては図 2.2.2 のような直交座標系が導入された. もちろん, 今 2 層モデルを考えているので, 大気は鉛直方向に対して 2 層に分割されていることに注意されたい. また, この矩形領域にはベータ平面近似 (2.2.46) が適用されるため, 南北方向の惑星渦度勾配が線形近似で考慮されている.

2.2.4 非断熱加熱の表現

まず, 地球大気のエネルギーバランスに関する一般的知識について簡潔に述べることにする. 地球大気の温度がほぼ一定に保たれるような平衡状態にあるためには, 地球が吸収する太陽の放射エネルギーと地球が放出する赤外放射のエネルギーがバランスしなければならない. 実際のところ, 獲得する太陽放射と放出する赤外放射に関するエネルギー収支は, 低緯度ではエネルギーを獲得し, 高緯度ではエネルギーを損失している. しかし, 現実の大気で永続的に温度が低緯度で上昇したり高緯度で下降したりしないのは, 上のエネルギー収支がゼロとなるように, 低緯度で余分に獲得された熱エネルギーが高緯度に運ばれるからである. Phillips の数値実験ではそのことが考慮されるように, 非断熱加熱率を決めている. ここでは, Phillips (1956) に沿って非断熱加熱の効果の表現法を説明していく.

(2.2.67) の非断熱効果の項 dQ/dt とは, 厳密に言えば 500 hPa 面における単位質量あたりの非断熱加熱率である. しかしながら, 2 層モデルが大変理想化されたモデルであることを考慮すれば, その非断熱加熱はコラム全体で平均したものとして解釈されるべきである. 故に, 小規模な乱流による鉛直方向の熱拡散による寄与は考えないことにする. したがって, 残りの dQ/dt を決める重要な物理過程は,

• 大気成分 (主に水蒸気) の状態変化に伴う潜熱の解放

• 放射過程

• 大規模な渦よりも小さなスケールにおける水平渦拡散

となる. 今回, 3 番目の過程は数値的に計算される.

水蒸気の蒸発や凝結は, 飽和領域では事実上の鉛直安定度である (θ_{1} − θ_{3})/θ_{2} を減少させる効果を持つ. さらに, ある緯度帯において凝結量が蒸発量を上回る場合, それはその緯度帯の大気が非断熱加熱を受けたことに対応する. しかし, Phillips の数値実験ではこのような取り扱いの難しい水収支を厳密にモデル化せず, 代わりに次のように対処した. すなわち潜熱解放の効果は, (2.2.67) 中の鉛直安定度 (θ_{1} − θ_{3})/θ_{2} を標準的な大気がもつ値の 8 割に取り, 同時に, 観測される潜熱の南北輸送をおおまかに許容するように加熱の南北勾配を減少させることによってのみ定性的に与えられるようにした.

図 2.2.3 は, 想定される年平均の放射エネルギーの獲得および損失にバランスするのに必要な, 南北方向のエネルギー輸送の計算結果をいくつか示している. 全ての結果において 30 度と 40 度の間にエネルギー輸送の極大が認められるが, Houghton

図 2.2.3: 放射の年平均における超過および損失とバランスするのに必要となるような, 緯度円全体を横切る極方向へのエネルギー輸送. 単位は 10^{19} cal day^{−1} である. 曲線 H, A はそれぞれ, Houghton (1954) と Albrecht (1931) によって得られた値である. Gabites (1950), Baur and Phillips (1935), Möller (1950) によって推定された曲線の極大値の位置も, 30 度から 40 度の緯度帯に存在している. (この図は Phillips(1956) から引用した. )

と Simpson による値と Albrecht と London による値とでは, 2 倍近くの差がある. 一方, 図 2.2.4 は, 観測される顕熱および潜熱の形式での南北方向のエネルギー輸送の計算結果をいくつか示している. この 2 つの図を比較すれば, 中高緯度では必要とされるエネルギー輸送の大部分が, 渦によってなされることが分かる. それは, 例え放射がより強い夏季でさえも当てはまる. 一方で, 低緯度において必要とされる南北エネルギー輸送は, 中高緯度とは別の仕組み, すなわち直接的な子午面循環によってなされることが考えられる.

次に, これらの観測的事実をもとに非断熱加熱の数学的表現を決定していくことにする. 「放射」によるエネルギーの獲得および損失を, y について線形な関数で次のように表現することにする.

(\frac{dQ}{dt})_{rad} = −2H \frac{y}{W}  (2.2.68)

図 2.2.4: 観測された渦による顕熱 (実線) および潜熱 (点線) の極方向への輸送. 単位は 10^{19} cal day^{−1} である. M(W) と M(S) は, 冬季と夏季の Mintz による地衡風の観測値である (Mintz 1955). 実線 S & W は Starr と White (1955) による年平均の値である. また, 曲線 S & W は Starr and White (ibid.) による年平均された渦の潜熱輸送の観測値である. 点線 BEW は潜熱の全輸送に対する Benzou and Estoque (1954), Wüst による推定の平均値を示している. (この図は Phillips(1956) から引用した. )

ここで, H は定数であり, y = −W から y = 0 までの単位質量あたりの平均的な加熱率 (あるいは y = 0 から y = W までの平均的な冷却率) である. よって, この加熱は運動に依存しないことになる. この条件は回転皿実験と実際の大気いずれにおいても完全には満たされないが, 第 1 近似としては問題ないと思われる.

今, 考えている矩形領域の南半分では放射によるエネルギーの獲得, 北半分では放射による損失が存在するので, 大気の温度が平衡状態にあるためには, 南半分の領域で余計に得た熱量を熱量が失われる北半分の領域へ輸送しなければならない. その時, 必要とされる y = 0 の直線を横切る熱エネルギーの輸送量は, 単位時間当たり,

∫_{0}^{L} dx ∫_{−W}^{0} dy ∫_{p_{4}}^{0} \frac{dp}{g} (−2H \frac{y}{W} ) = \frac{HW Lp_{4}}{g}  (2.2.69)

と計算される. 先程示した図によれば, 45 度の緯度円で観測される渦の顕熱輸送は, 単位長さ当たり 10^{5} kJ s^{−1} である. よって, W = 5 × 10^{6} m(地球の極と赤道間の距離の約半分) として, H の値を次のように設定すればよいだろう.

H ≃ \frac{10^{5} kJ s^{−1}}{W p_{4}} = 2 × 10^{−6} kJ ton^{−1}s^{−1}  (2.2.70)

図 2.2.4 による渦の顕熱および潜熱の輸送の合計は図 2.2.3 の 45 度における平衡状態に必要なエネルギー輸送とおおよそ等しいので, ここで計算された H の値は, 観測される正味の放射の効果から観測される潜熱輸送の収束の効果を除いたものにおおよそ対応することになる.

次に, dQ/dt における水平渦拡散の効果は, (2.2.66) を用いて次のように表現される.

(\frac{dQ}{dt})_{diff} = c_{p}A_{T} ∇^{2}T_{2} = \frac{c_{p}}{R} A_{T} ∇^{2}(Φ_{1} − Φ_{3})  (2.2.71)

これにより, 小規模な渦による水平渦拡散の効果は定量的に数値計算される.

最後に, 水平方向の動渦粘性係数 A_{v} と水平方向の熱の渦拡散係数 A_{T} について考えることにする. 今回の場合, これらは計算に使われる格子サイズよりも小さい渦の効果を表現する. Richardson の経験則から A ∼ 0.2l^{4/3} ( Richardson 1926, Batchelor 1950) であり, 今 l = Δx ∼ 30 × 10^{7}cm とすれば, A は 10^{5} m^{2}s^{−1} ぐらいの値を持つ. 一般的な混合距離理論 A ∼ v̄l において, これは実質的な (小規模の) 渦の速度が 1 m/s よりも小さいことに対応している. A_{v} や A_{T} の選び方についての経験則や基礎理論が当時はまだ不足していたために, Phillips の数値実験では両者ともある一定の値 A が使われた. すなわち,

A_{v} = A_{T} = A = 10^{5} m^{2} s^{−1}.  (2.2.72)

なお, 上のように取られた A は, 大規模な熱交換過程 (Defant, 1921) として図 2.2.3 で必要とされる熱輸送を説明するのに必要な値の約 1/40 程しかない. 故に, 南北の温度勾配が通常の 40 倍にでもならない限り, 熱拡散効果によって高緯度の放射によるエネルギーの損失を埋めることはできないだろう.

2.2.5 数値実験で用いられた予報方程式系と境界条件

予報方程式系

前述の渦度方程式と熱力学的エネルギーの式の二本で予報方程式系を成しているが, Phillips の数値実験では, それらの方程式から鉛直速度を消去したポテンシャル渦度方程式が採用された. ここでは, そのポテンシャル渦度方程式を導き最終的な予報方程式系を決定する.

最初に地衡風に対する流線関数を定義する.

ψ_{1} = Φ_{1}/f_{0}, ψ_{3} = Φ_{3}/f_{0}  (2.2.73)

これにより, 渦度方程式に現われる地衡風の渦度は,

ζ_{1} = ∇^{2}ψ_{1}, ζ_{3} = ∇^{2}ψ_{3}  (2.2.74)

と表現される. さらに, 鉛直方向の静的安定度を表すパラメータとして,

λ^{2} ≡ \frac{f_{0}^{2}θ_{2}}{(Φ_{1} − Φ_{3})(θ_{1} − θ_{3})} = constant  (2.2.75)

を導入する.

レベル 1 とレベル 3 で適用される渦度方程式 (2.2.61) を再掲すれば,

(\frac{∂}{∂t} + v_{1} · ∇) (f + ζ_{1}) − f_{0} \frac{ω_{2}}{p_{2}} = A∇^{2}ζ_{1}  (2.2.76)
(\frac{∂}{∂t} + v_{3} · ∇) (f + ζ_{3}) + f_{0} \frac{ω_{2}}{p_{2}} = A∇^{2}ζ_{3} − kζ_{4}  (2.2.77)

である. 上式において, A = A_{v} = A_{T} としている. また, レベル 2 で適用される熱力学的エネルギー保存の式 (2.2.67) は, 前小節で決められた非断熱加熱の数学的表現と鉛直方向の静的安定度のパラメータ λ^{2} および地衡風の流線関数を用いて,

\frac{f_{0}ω_{2}}{p_{2}} = λ^{2} [(\frac{∂}{∂t} + v_{2} · ∇) (ψ_{1} − ψ_{3}) + \frac{2RH}{f_{0}c_{p}} \frac{y}{W} − A∇^{2}(ψ_{1} − ψ_{3})] (2.2.78)

と書かれる.

ここで, 相対的なポテンシャル渦度を,

\begin{aligned}q_{1} & = ζ_{1} − λ^{2}(ψ_{1} − ψ_{3}) = ∇^{2}ψ_{1} − λ^{2}(ψ_{1} − ψ_{3}) \\ q_{3} & = ζ_{3} + λ^{2}(ψ_{1} − ψ_{3}) = ∇^{2}ψ_{3} + λ^{2}(ψ_{1} − ψ_{3}) \end{aligned}  (2.2.79)

と定義する. (2.2.76) と (2.2.78) から鉛直速度 ω_{2} を消去し, 定義した q_{1}, q_{3} を用いれば今回の数値予報モデルにおけるポテンシャル渦度方程式が得られる^{*13} .

\frac{∂q_{1}}{∂t} = −v_{1} · ∇(βy + q_{1}) + A∇^{2}q_{1} + \frac{2RHλ^{2}}{f_{0}c_{p}} \frac{y}{W}  (2.2.80)

\frac{∂q_{3}}{∂t} = −v_{3} · ∇(βy + q_{3}) + A∇^{2}q_{3} − \frac{2RHλ^{2}}{f_{0}c_{p}} \frac{y}{W}

−\frac{k}{2}[3q_{3} − q_{1} − 4λ^{2}(ψ_{1} − ψ_{3})]  (2.2.81)

よって, 時刻 t の q と ψ の場を知っているならば, (2.2.80) と (2.2.81) は ∂q/∂t を与えるので, 外挿によって時刻 t + Δt の q の値を得ることができる. なお, 実際に数値計算を行うためには上の方程式系を水平方向・時間方向に対して離散化しなければならない. さらに適用する数値計算スキームに関していくつかの工夫を加えることになるが, これらのことは「数値予報モデルの構築 : 数値計算手法」の章にて説明する.

境界条件

(2.2.80) と (2.2.81) を数値計算で解く際に, 矩形領域 0 ≤ x ≤ L, −W ≤ y ≤ W の境界において境界条件が与えられる必要がある. よって, この章の最後として Phillips の数値実験で適用された境界条件について述べることにする.

前述したように, 東西端 x = 0, L に対しては周期境界条件が適用される. 一方で, 南北端 y = ±W に対しての境界条件はあまり明白でない. たいていの短期間の数値予報 (領域モデル) では水平境界に対しては任意の条件が課されるが, 不正確な境界条件のために正しく予報を行える領域は時間とともに減少する. Phillips の数値実験ではより長い間時間積分を実行することになるのでそのようなことはできず, 南北端の壁において適切な境界条件が与えられなければならない. ここで言う「壁」とは主には運動を制限するための数学的な仕掛けであり, 実際の地球上で対応するものはない. ここで採用される境界条件は性質上いくらか発見的なものであるが, 大気の運動に対するエネルギーや運動量の積分に関して重要な意味を持つ. このことは, 後ほど述べることにする.

*13導出の詳細は付録 D を参照.

\overline{main.tex 2011/02/28(河合佑太)}

南北端で適用される境界条件について具体的に説明するにあたって, まず流れの場を平均場と擾乱場に分けることにする. すなわち, ψ の平均場を

ψ̄(y, t) ≡ \frac{1}{L} ∫_{0}^{L} ψ(x, y, t)dx,  (2.2.82)

ψ の擾乱場を

ψ′(x, y, t) ≡ ψ(x, y, t) − ψ̄(y, t)  (2.2.83)

と定義する. なお, この定義は q, ζ, u, v, ω など他の従属変数に対しても同様に適用される.

最初に剛体壁に与えられる一般的な境界条件として,

ψ′ ≡ 0 at y = ±W  (2.2.84)

を課す. これは y = ±W において, 壁に垂直な地衡風成分はゼロであることに対応する. 次に与える境界条件は, いくらか任意ではあるが壁において渦度がゼロとなることを課すことにする. すなわち,

ζ′ ≡ 0 at y = ±W.  (2.2.85)

次に, y = ±W においてレベル 1 および 3 の (2.2.42) の東西平均をとれば^{*14},

\frac{∂ū_{1}}{∂t} = −\frac{∂}{∂t} \frac{∂ψ̄_{1}}{∂y} = A \frac{∂^{2}ū_{1}}{∂y^{2}} ,
\frac{∂ū_{3}}{∂t} = −\frac{∂}{∂t} \frac{∂ψ̄_{3}}{∂y} = A \frac{∂^{2}ū_{3}}{∂y^{2}} − \frac{k}{2}(3ū_{3} − ū_{1})  (2.2.86)
\frac{∂u_{1}}{∂t} + u_{1} \frac{∂u_{1}}{∂x} = −\frac{∂Φ_{1}}{∂x} + A∇^{2}u_{1},
\frac{∂u_{3}}{∂t} + u_{3} \frac{∂u_{3}}{∂x} = −\frac{∂Φ_{3}}{∂x} + A∇^{2}u_{3} − k \frac{3u_{3} − u_{1}}{2}

\begin{aligned}\overline{u\frac{∂u}{∂x}} & = \frac{1}{2} (\overline{\frac{∂u^{2}}{∂x}} ) \\ & = \frac{1}{2L} ∫_{0}^{L} \frac{∂}{∂x} u^{2} dx = \frac{1}{2L} [u^{2}]_{0}^{L} = 0 \end{aligned}

*14南北境界において v = 0, x = 0, L において周期境界条件が適用されることを用いれば, (2.2.42)

は.

となる. ここで, u_{1}, u_{3} を使って u_{4} を外挿した.

上の 2 式の東西平均をとるとき, 東西平均の x 微分の項は周期境界条件によってゼロとなる. 例えば, 移流項はとなる. 他の項についても同様の考慮を行えば (2.2.86) を得る.

\overline{main.tex 2011/02/28(河合佑太)}

の関係が導かれる. ここですべての時刻 t について,

\frac{∂^{2}ū}{∂y^{2}} ≡ 0 at y = ±W  (2.2.87)

であることを仮定しよう. これにより, 最初に南北境界において ū_{1}, ū_{3} がゼロならば, (2.2.86) からその後もゼロであり続ける. したがって最後に挙げる南北境界における境界条件は次のようになる.

\frac{∂ψ̄}{∂y} ≡ 0 at y = ±W.  (2.2.88)

(2.2.67) の左辺において (Φ_{1} − Φ_{3}) を定数として取り扱うことを仮定したために, それ以後の式では ψ_{1} と ψ_{3} は時間微分・空間微分, または ψ_{1} − ψ_{3} と何らかの項との積の形式でのみ式中に現れる. よって, 例えば y = −W で ψ̄_{3} の値を固定しても問題ない. 予報方程式系 (2.2.79) ∼ (2.2.81) と境界条件 (2.2.84) ∼ (2.2.88) により, 運動は完全に決定される.

2.3 系のエネルギー論

2.3.1 エネルギー方程式

摩擦・非断熱加熱を含む準地衡風予報方程式系および前節で与えた境界条件から, このモデルにおける 4 本のエネルギー方程式を導く.

最初に, (2.2.76) と (2.2.78) の東西平均をとれば,

\frac{∂^{2}}{∂y^{2}} \frac{∂ψ̄_{1}}{∂t} − ω̄_{2} \frac{f_{0}}{p_{2}} = A\frac{∂^{2}ζ̄_{1}}{∂y^{2}} − \overline{v1^{′} · ∇ζ1^{′}},  (2.3.89)
\frac{∂^{2}}{∂y^{2}} \frac{∂ψ̄_{3}}{∂t} + ω̄_{2} \frac{f_{0}}{p_{2}} = A\frac{∂^{2}ζ̄_{3}}{∂y^{2}} − \overline{v3^{′} · ∇ζ3^{′}} − kζ̄_{4},  (2.3.90)

\frac{ω̄_{2}f_{0}}{p_{2}} = λ^{2}[\frac{∂(ψ̄_{1} − ψ̄_{3})}{∂t} + \overline{v1^{′} · ∇(ψ1^{′} − ψ3^{′})} + \frac{2RH}{f_{0}c_{p}} \frac{y}{W} − A\frac{∂^{2}(ψ̄_{1} − ψ̄_{3})}{∂y^{2}}].(2.3.91)

を得る. 続いて擾乱場の方程式を得るために, (2.2.76), (2.2.78) の 3 式を上の 3 式でそれぞれを引けば,

\overline{main.tex 2011/02/28(河合佑太)}

∇^{2} \frac{∂ψ1^{′}}{∂t} − ω2^{′} \frac{f_{0}}{p_{2}} = A∇^{2}ζ1^{′} − [v_{1} · ∇(βy + ζ_{1})]′,  (2.3.92)
∇^{2} \frac{∂ψ3^{′}}{∂t} + ω2^{′} \frac{f_{0}}{p_{2}} = A∇^{2}ζ3^{′} − [v_{3} · ∇(βy + ζ_{3})]′ − kζ4^{′},  (2.3.93)
\frac{ω2^{′}f_{0}}{p_{2}} = λ^{2} (\frac{∂(ψ1^{′} − ψ3^{′})}{∂t} + [v_{3} · ∇(ψ_{1} − ψ_{3})]′ − A∇^{2}(ψ1^{′} − ψ3^{′})) (2.3.94)

となる.

系がもつ平均帯状流の運動エネルギーの時間変化を表す方程式を導く. まず (2.3.89) と (2.3.90) にそれぞれ −ψ̄_{1} と −ψ̄_{3} を掛けて加え合わせる. そして, その結果を境界条件に注意しながら領域全体で積分すれば, 平均帯状流の運動エネルギー K̄ の時間変化率に対する方程式^{*15}

\begin{aligned}\frac{∂K̄}{∂t} ≡ & \frac{∂}{∂t} ∫ \frac{ū_{1}^{2} + ū_{3}^{2}}{2} dσ \\ = & − \frac{f_{0}}{p_{2}} ∫ ω̄_{2}(ψ̄_{1} − ψ̄_{3})dσ − A ∫ (ζ̄_{1}^{2} + ζ̄_{3}^{2})dσ \\ & − k ∫ ū_{3}ū_{4}dσ − ∫ [ū_{1} \frac{∂(\overline{u1^{′}v1^{′}})}{∂y} + ū_{3} \frac{∂(\overline{u3^{′}v3^{′}})}{∂y}] dσ \end{aligned}  (2.3.95)

を得る. ここで, dσ は面素 dx dy である.

次に系がもつ擾乱の運動エネルギーの時間変化を表す方程式を上と同様の方法で導く. (2.3.92) と (2.3.93) にそれぞれ −ψ1^{′} と −ψ3^{′} を掛けて加え合わせる. そして, その結果を境界条件に注意して領域全体で積分すれば, 擾乱の運動エネルギー K′ の時間変化率に対する方程式

\begin{aligned}\frac{∂K′}{∂t} ≡ & \frac{∂}{∂t} ∫ [(∇ψ1^{′})^{2} + (∇ψ3^{′})^{2}]dσ \\ = & − \frac{f_{0}}{p_{2}} ∫ ω2^{′}(ψ1^{′} − ψ3^{′})dσ − A ∫ [(ζ1^{′})^{2} + (ζ3^{′})^{2}]dσ \\ & − k ∫ v3^{′}v4^{′}dσ + ∫ [ū_{1} \frac{∂(\overline{u1^{′}v1^{′}})}{∂y} + ū_{3} \frac{∂(\overline{u3^{′}v3^{′}})}{∂y}] dσ \end{aligned}  (2.3.96)

を得る.

続いてポテンシャルエネルギーの時間発展の式を得るには, (2.3.91) と (2.3.94) にそれぞれ (ψ̄_{1} − ψ̄_{3}), (ψ1^{′} − ψ3^{′}) を掛けて, 予報領域全体で積分すればよい. これにより, 平均流のポテンシャルエネルギー P̄ の時間発展式

*15この節のエネルギー方程式の導出の詳細は, 付録 E を参照.

\overline{main.tex 2011/02/28(河合佑太)}

\begin{aligned}\frac{∂P̄}{∂t} ≡ & \frac{1}{2}λ^{2} \frac{∂}{∂t} ∫ (ψ̄_{1} − ψ̄_{3})^{2}dσ \\ = & \frac{f_{0}}{p_{2}} ∫ ω̄_{2}(ψ̄_{1} − ψ̄_{3})dσ + λ^{2} ∫ \overline{v1^{′}(ψ1^{′} − ψ3^{′})}\frac{∂(ψ̄_{1} − ψ̄_{3})}{∂y}dσ \\ & − \frac{2RHλ^{2}}{f_{0}c_{p}} ∫ \frac{y}{W} (ψ̄_{1} − ψ̄_{3})dσ − λ^{2}A ∫ [\frac{∂(ψ̄_{1} − ψ̄_{3})}{∂y}]^{2} dσ \end{aligned}  (2.3.97)

と擾乱運動のポテンシャルエネルギー P′ の時間発展の式

\begin{aligned}\frac{∂P′}{∂t} ≡ & \frac{1}{2}λ^{2} \frac{∂}{∂t} ∫ (ψ1^{′} − ψ3^{′})^{2}dσ \\ = & \frac{f_{0}}{p_{2}} ∫ ω2^{′}(ψ1^{′} − ψ3^{′})dσ − λ^{2} ∫ \overline{v1^{′}(ψ1^{′} − ψ3^{′})}\frac{∂(ψ̄_{1} − ψ̄_{3})}{∂y}dσ \\ & − λ^{2}A ∫ [\frac{∂(ψ̄_{1} − ψ̄_{3})}{∂y}]^{2} dσ \end{aligned}  (2.3.98)

を得る.

以上 4 本のエネルギー方程式 (2.3.95)∼(2.3.98) を結合すれば, 系がもつ全エネルギーの時間発展式を得る.

\begin{aligned}\frac{∂}{∂t} (P̄ + P′ + K̄ + K′) = & − A ∫ (ζ_{1}^{2} + ζ_{3}^{2}) dσ − A ∫ [∇(ψ_{1} − ψ_{3})]^{2}dσ \\ & − k ∫ v_{3} · v_{4}dσ − \frac{2RHλ^{2}}{f_{0}c_{p}} ∫ \frac{y}{W} (ψ̄_{1} − ψ̄_{3})dσ. \end{aligned} (2.3.99)

上式によれば, このモデルにおける全エネルギーは次の要因によって変化することになる.

• 水平方向の渦粘性の効果 A によるエネルギーの損失.

• 地表摩擦の効果 k による変化. v_{3} と v_{4} は同じ方向を向く傾向があるので, ふつう地上摩擦の効果はエネルギーを減少させる.

• 非断熱加熱の効果 H による変化. 加熱が温度 (温度は ψ_{1} − ψ_{3} に比例する) と正の相関を持つとき, この効果はエネルギーを増加させる. 一方で負の相関を持つときはエネルギーを減少させる.

(2.2.68) から系に加えられる熱エネルギーの合計はゼロである. しかしながら, 領域の高温部で加熱され, 低温部で冷却される場合には, このモデルの有効位置エネルギー (available potential energy(APE)) は増加する. このとき系のエントロピーは形式的に減少する.

2.3.2 系のエネルギー変換

前小節で導いた 4 本のエネルギー方程式 (2.3.95) ∼ (2.3.98) の各項を見比べてみれば, 同じ形で符号だけ異なる項が存在することが分かる. これらの項はあるエネルギーから他のエネルギーへの変換を表現していると考えられる. そこで, このエネルギー変換をより簡潔に表現するために, {P̄ · K̄} のような記号を導入することにする. このエネルギー変換を表わす記号は, 括弧内の左の記号で表わすエネルギー形式から括弧内の右の記号で表わすエネルギー形式にエネルギーが変換されることを表現する項であることを示している.

具体的には, 次のようにエネルギー変換の表記を定義する.

\{Q̄ · P̄\} ≡ −\frac{2RHλ^{2}}{f_{0}c_{p}} ∫ \frac{y}{W} (ψ̄_{1} − ψ̄_{3}) dσ, ,
\{P̄ · P′\} ≡ −λ^{2} ∫ [\overline{v1^{′}(ψ1^{′} − ψ3^{′})}]\frac{∂(ψ̄_{1} − ψ̄_{3})}{∂y} dσ,
\{P′ · K′\} ≡ −\frac{f_{0}}{p_{2}} ∫ ω2^{′}(ψ1^{′} − ψ3^{′}) dσ,
\{K′ · K̄\} ≡ − ∫ [ū_{1} \frac{∂(\overline{u1^{′}v1^{′}})}{∂y} + ū_{3} \frac{∂(\overline{u3^{′}v3^{′}})}{∂y}] dσ,
\{P̄ · K̄\} ≡ −\frac{f_{0}}{p_{2}} ∫ ω̄_{2}(ψ̄_{1} − ψ̄_{3}) dσ,
{K̄ · A} ≡ A ∫ (ζ̄_{1}^{2} + ζ̄_{3}^{2})dσ,
{K′ · A} ≡ A ∫ [(ζ1^{′})^{2} + (ζ3^{′})^{2}]dσ,
\{P̄ · A\} ≡ λ^{2}A ∫ [\frac{∂(ψ̄_{1} − ψ̄_{3})}{∂y}]^{2} dσ,
{P′ · A} ≡ λ^{2}A ∫ [∇(ψ1^{′} − ψ3^{′})]^{2} dσ,
{K̄ · k} ≡ k ∫ ū_{3}ū_{4} dσ,
{K′ · k} ≡ k ∫ v3^{′} · v4^{′} dσ.

これらの表記を用いれば, 4 本のエネルギー方程式 (2.3.95) ∼ (2.3.98) と全エネルギー E = K̄ + K′ + P̄ + P′ の方程式 (2.3.99) は次のように書ける.

\begin{matrix}\frac{∂K̄}{∂t} = & \{P̄ · K̄\} − \{K̄ · A\} − \{K̄ · A\} + \{K′ · K̄\}, \\ \frac{∂K′}{∂t} = & \{P′ · K′\} − \{K′ · A\} − \{K′ · k\} − \{K′ · K̄\}, \\ \frac{∂P̄}{∂t} = & − \{P̄ · K̄\} − \{P̄ · P′\} + \{Q̄ · P̄\} − \{P̄ · A\}, \\ \frac{∂P′}{∂t} = & − \{P′ · K′\} + \{P̄ · P′\} − \{P′ · A\}, \\ \frac{∂K̄}{∂t} = & \{P̄ · K̄\} − \{K̄ · A\} − \{K̄ · A\} + \{K′ · K̄\}, \\ \frac{∂E}{∂t} = & − \{K · A\} − \{K · k\} − \{P · A\} + \{Q̄ · P̄\}. \end{matrix}  (2.3.100)

これにより 4 本のエネルギー方程式および全エネルギーの方程式の物理的な解釈は, ずっと簡単になる. 最初に, このモデルのポテンシャルエネルギーの表現は 500 hPa での標準大気からの温度のずれの平方の積分 (APE) で本質的には構成される.

上で表記を決めた 11 種類のエネルギー変換についてその意味を考えることにしよう.

{P̄ · P′} : 帯状平均流のポテンシャルエネルギーから擾乱のポテンシャルエネルギーへの変換. この変換は, 暖域から寒域への擾乱の南北速度 v1^{′} または v3^{′} による水平方向の顕熱の輸送によってなされる.

{P̄ ·K̄} : 帯状平均流のポテンシャルエネルギーから平均帯状流の運動エネルギーへの変換. この変換は, 子午面における鉛直方向の循環の効果を表現している.

{P′ · K′} : 擾乱のポテンシャルエネルギーから擾乱運動の運動エネルギーへの変換. この変換は, 東西平面における鉛直方向の循環の効果を表現している.

{K′ · K̄} : 擾乱の運動エネルギーから帯状平均流の運動エネルギーへの変換. 例えば, 大きな ū の値 (強い西風の平均流) をもつ領域で, 渦による運動量フラックスの東西平均\overline{u′v′} 場が南北方向に収束である場合を考える. このとき \{K′ · K̄\} は正となり, 擾乱の運動エネルギーは帯状平均流の運動エネルギーに変換される (帯状平均流は加速される).

図 2.3.5: エネルギーダイアグラム. 定義されたエネルギー変換が正である場合, エネルギーは上の図において矢印の方に流れる. 今回の数値計算では, {P̄ · K̄} を除いて全てのエネルギー変換は正として計算される.

{K̄ · A}, {K′ · A} : 水平方向の渦粘性による運動エネルギーの損失. 非圧縮な粘性流体と同じ形式となる.

{P̄ · A}, {P′ · A} : 水平方向の渦粘性によるポテンシャルエネルギーの損失. 温度勾配を平滑化するような拡散プロセスによるポテンシャルエネルギーの損失を表現している.

{K̄·k}, {K′·k} : 地面の摩擦による運動エネルギーの損失. この項は−k ∫ ψ_{3}ζ_{4} dσ の形式で書かれているが, これは摩擦層で地表摩擦によって生み出される収束や発散に対する, 摩擦層より上部の流体の間接的な応答として解釈される (エクマン・パンピング). 流れのパターンが高さとともに劇的に位相を変化させないときには, エクマン・パンピングにより結果的に各層の運動エネルギーを減少させる.

次にこの数値予報モデルにおけるエネルギー変換を, 図 2.3.5 のエネルギーダイアグラムを通して理解してみる. 4 つの箱は 4 つのエネルギーの種類に対するエネルギーのリザーバーを表し, その間に引かれている矢印は上で定義したエエネルギー形式の変換を表現している.

定常状態では, それぞれのリザーバーに蓄積されるエネルギーは一定でなければならない. つまり, それぞれのエネルギーリザーバーにおいて流入するエネルギーと流出されるエネルギーの合計がゼロでなければならない. 古典的な大循環の理論では子午面循環を維持する上で大規模な渦の効果を考慮しなかったので, 摩擦による散逸に対して平均流の運動エネルギー K̄ を維持するのに十分なエネルギーを供給するには, エネルギー変換 {P̄ · K̄} は正でなければならなかった. このようなエネルギー変換は, 高温な低緯度では ω̄_{2} < 0(上昇流), 低温な高緯度では ω̄_{2} > 0(下降流) であるような直接的な子午面循環によってなされるだろう.

しかしながら, 実際のところ観測により {P̄ · K̄} はわずかに負であることが知られている (Starr(1954)). 本モデルにおいて K̄ に対する他のエネルギー源は K′ であるので, {K′ · K̄} の測定が, エネルギーの流れの経路を決める上で重要となる. 観測によれば, {K′ · K̄} は正であり, 摩擦による K̄ の損失を埋めるのに十分な大きさであることが示されている (Arakawa 1953, Starr 1953). したがって K′ が {K′ · K̄} によって損失されるならば, K′ は P′ からエネルギーを獲得しなければならない. すなわち, 東西平面において冷たい空気は下降し, 暖かい空気は上昇する傾向が見られなければならない. しかしながら, 鉛直運動を観測するのは困難なためにこのプロセスの客観的な観測は難しい.

最後に, もし P′ が K′ にエネルギーを供給するならば, P′ はエネルギー変換 {P̄ · P′} によって P̄ からエネルギーを受け取らなければならない. この変換が正でなければならないことは容易に想像される. もう一度 {P̄ · P′} の式を見てみよう. \overline{v1^{′}(ψ1^{′} − ψ3^{′})} は渦による顕熱の極方向への輸送に比例し, それは図 2.2.4 から正であることが知られている. また, ∂(ψ̄_{1} − ψ̄_{3})/∂y は ∂T̄_{2}/∂y に比例するので負である. したがって, {P̄ · P′} が正であることが示される.

2.4 子午面循環

数値モデルにおける子午面循環は, ω に対する境界条件 (2.2.41) を用いて, 連続の式 (2.2.44) から計算される.

最初に平均南北速度を

\begin{matrix}V̄_{1} = \frac{1}{Lp_{2}} ∫_{0}^{p_{2}} dp ∫_{0}^{L} v dx, \\ V̄_{3} = \frac{1}{Lp_{2}} ∫_{p_{2}}^{p_{4}} dp ∫_{0}^{L} v dx \end{matrix}

と定義する. 連続の式 (2.2.44) をレベル 1 とレベル 3 において適用し, 境界条件 (2.2.41) に考慮して鉛直方向に離散化すれば, 連続の式は

∇ · v_{1} = −∇ · v_{3} = −\frac{ω_{2}}{p_{2}}  (2.4.101)

となる. この式の帯状平均をとれば,

\frac{∂V̄_{1}}{∂y} = −\frac{∂V̄_{3}}{∂y} = −\frac{ω̄_{2}}{p_{2}}  (2.4.102)

を得る. このとき, ψ_{1} と ψ_{3} およびその微分が既知であるならば ω̄_{2}/p_{2} は, (2.3.89), (2.3.90), (2.3.91) のうちどれかの方程式により評価される.

今, 南北端の境界条件から, y = ±W において V̄_{1,3} = 0 である. したがって, (2.4.102) を y = −W から y = W まで y について積分すれば, ∫_{−W}^{W} ω̄_{2} dy = 0 であることが示される. また, (2.4.102) より V̄_{1} + V̄_{3}(= 2V̄_{2}) は定数であるので, 南北端の境界条件をレベル 2 に適用すれば, V̄_{1} + V̄_{3} = 0 を得る. したがって, レベル 1 とレベル 3 の平均南北速度について V̄_{1} = −V̄_{3} の関係がある. 子午面循環 V̄ は, 東西平均をとった渦度方程式中の発散項にのみ暗黙的に現われる. 一方, 残りの項は地衡風的な力学場を使って評価される.

東西風 ū に対してのみ地衡風近似を導入することによって Eliassen(1952) はこれらの効果をある程度認める軸対称の渦に対する方程式を導いた. しかしながら, 完全な準地衡風モデルにおいてこれらの特性を保持しながら, (2.2.76), (2.2.78) の内部的な一貫性を維持することは難しいと思われる. よって作成したモデルから計算される V̄ は, 大きな T̄ と ū の水平勾配を持つ領域または f が小さな赤道域の大気にのみ部分的に適用されるであろう.

東西平均をとった方程式によれば, 仮定上の定常状態における子午面循環は熱や表面摩擦といった「外部的な」要素のみでなく, \overline{v′ · ∇ζ′} や\overline{v′ · ∇T′} にも依存する. 後者の値は, 少なくとも中緯度地域では前者の値と同じぐらいの大きさのオーダーを持つことが観測されている. そのため, 平均子午面循環の議論では外部的な要素と同じぐらい渦による効果を考えることが重要である.

第3章{数値予報モデルの構築 : 数値 \atop 計算手法}

この章では, 実際に用いられた数値計算手法の詳細について説明することにする.

3.1 予報方程式系の離散化

Phillips の数値実験では格子点法が用いられた. 予報領域である 2 層の矩形領域は格子によって分割され, 各物理量は各格子点上で与えられる. 最初に, 連続的な独立変数を次のように離散変数に置き換える.

x = iΔx (i = 0, 1, 2, ..., I − 1),
y = jΔy (j = 0, 1, 2, ..., J),
t = τΔt (τ = 0, 1, 2, ...).

i = 0 は x = 0, i = I は x = L に対応している. また, j = 0 は y = −W, j = J は y = W に対応している.

ラプラシアンおよびヤコビアンの有限差分形式の演算子を次のように定義することにする.

Δx^{2}∇^{2}S ≃ ℒ_{ij}S = S_{i+1 j} − 2S_{ij} + S_{i−1 j} + ε^{2}(S_{i j+1} − 2S_{ij} + S_{i j−1}),

4ΔxΔy𝒥 (\frac{R, S}{x, y} ) ≃ 𝒥_{ij}(R, S) = (R_{i+1 j} − R_{i−1 j})(S_{i j+1} − S_{i j−1})

−(R_{i j+1} − R_{i j−1})(S_{i+1 j} − S_{i−1 j}).

ここで, ε = Δx/Δy である. さらに, 便宜上無次元の定数

γ = λ^{2}Δx^{2}

を導入し, またポテンシャル渦度 q を,

η_{ij} = q_{ij} Δx^{2}

と置き換えることにする.

x 方向の平均の定義 (2.2.82) および擾乱場の定義 (2.2.83) の有限差分形式は,

ψ̄_{jτ} = I ^{−1} ∑_{i=0}^{I−1} ψ_{ij τ},  (3.1.1)
ψijτ^{′} = ψ_{ijτ} − ψ̄_{jτ}  (3.1.2)

となる. (2.2.79) の擾乱場に対する式と東西平均場の式は次のように書かれる^{*1}.

ℒ_{ij}ψ1^{′} − γ(ψ1^{′} − ψ3^{′})_{ij} − η1 ij^{′} = 0,  (3.1.3)
ℒ_{ij}ψ3^{′} + γ(ψ1^{′} − ψ3^{′})_{ij} − η3 ij^{′} = 0,  (3.1.4)
ε^{2}(ψ̄_{1 j+1} − 2ψ̄_{1 j} + ψ̄_{1 j−1}) − γ(ψ̄_{1j} − ψ̄_{3j}) − η̄_{1j} = 0,  (3.1.5)
ε^{2}(ψ̄_{3 j+1} − 2ψ̄_{3 j} + ψ̄_{3 j−1}) + γ(ψ̄_{1j} − ψ̄_{3j}) − η̄_{3j} = 0.  (3.1.6)

ここでは, τ の添え字を省略した.

η1^{′} と η3^{′} が与えられたならば, (3.1.3) と (3.1.4) は次の境界条件を用いて領域内

(0 ≤ i ≤ I − 1, 1 ≤ j ≤ J − 1) の ψ1^{′}, ψ3^{′} を決定する.

ψ′(i, j) = 0 at j = 0 and J,  (3.1.7)
ψ′(0, j) = ψ′(I, j), ψ′(−1, j) = ψ′(I − 1, j).  (3.1.8)

前者の境界条件は (2.2.85) に対応し, 後者は周期境界条件である. このとき 2 本の方程式 (3.1.3), (3.1.4) は緩和法を用いて ψ′ に対して数値的に解かれる. 今回の再現では, 擾乱場に対する連立方程式を解くために緩和法の一つであるガウス・ザイデル法を用いた.

η̄_{j} が既知であるとき, 常微分方程式の差分形式 (3.1.5), (3.1.6) の ψ̄_{j} に対する解

q̄_{1} = \frac{∂ψ̄_{1}}{∂y} − λ^{2}(ψ̄_{1} − ψ̄_{3})

q1^{′} = ∇^{2}ψ1^{′} − λ^{2}(ψ1^{′} − ψ3^{′})

*1東西平均の式を求めるためには, (2.2.82) に基づいて方程式を x 方向に積分すればよい.

たとえば (2.2.79) の 1 式目

q_{1} = ∇^{2}ψ_{1} − λ^{2}(ψ_{1} − ψ_{3})

について東西平均を取れば, 東西平均場に対する方程式を得る.

また, もとの式を東西平均の式で引けば擾乱場に対する方程式を得る.

\overline{main.tex 2011/02/28(河合佑太)}

法には次の境界条件 (ノイマン境界条件)

ψ̄(j = J) = ψ̄(j = J − 1),
ψ̄(j = J) = ψ̄(j = J − 1),  (3.1.9)

および余分な条件, 例えば,

ψ̄_{3}(j = 0) = 0  (3.1.10)

を用いる.

伝統的な緩和法を用いた (3.1.5), (3.1.6) の解法はこの場合あまり満足でない. 方程式の 1 次元的特徴や境界条件 (3.1.9) によって, 緩和法の反復過程における急速な収束が抑制されてしまう. さらに鉛直速度 ω̄_{2j} の計算を正確に行うためには, むしろ緩和法を用いない正確な数値解法が要求される. Phillips (1956) では, これらの方程式を解くために H. Goldstine が示した方法^{*2}を用いたと述べられているが, その具体的な手法は論文中に記されていない. そのため, 今回はこの東西平均場に対する連立方程式を解くために直接解法を使用することで ψ̄_{j} を正確に求めることにした.

(2.2.80),(2.2.81) の有限差分形式を得る手法には, タイムレベル τ で方程式の右辺を評価し, 左辺を中心差分 (q_{ik τ+1} − q_{ij τ−1})/2Δt として表現する方法がある (ただし, 一番最初の時間積分には前進差分を使う必要がある). しかしながら, 今回の場合 A や k を含む摩擦項のために, この手法は時間積分を安定に行えないことが示されている. すなわち, タイムレベル τ でこれらの項を評価するならば, これらの項は指数関数的に増幅する解と減衰する解の両方を生み出す. 一方で, 減衰する解のみが物理モードである. この問題は, タイムレベル τ − 1 において陽に摩擦項を評価する, またはタイムレベル τ − 1 と τ + 1 の平均として陰に評価するかのどちらかによって達成される. 後者の方法がわずかにより正確であり, 今回の計算ではこちらの方法をスキームの基礎に置くことにする.

(2.2.80),(2.2.81) に対する実際の有限差分形式は次のようになる.

[−αℒ + 1] η_{1 ij τ+1} = χ_{1 ij} = η_{1 ij τ−1} + c𝒥_{ij}(bj + η_{1}, ψ_{1})_{τ}
+αℒ η_{1 ji τ−1} + h(2j − J)/J,  (3.1.11)
[−αℒ + 1 + \frac{3kΔt}{2}] η_{3 ij τ+1} = χ_{3 ij} = η_{3 ij τ−1} + c𝒥_{ij}(bj + η_{3}, ψ_{3})_{τ}
+αℒ η_{3 ji τ−1} + h(2j − J)/J  (3.1.12)
−kΔt [3/2 η_{3 τ−1} − η_{1τ} − 4γ(ψ_{1} − ψ_{3})_{τ}]_{ij} .

*2Phillips(1956) ではその手法は, 常微分方程式の差分形式のより広い分類に対して適用される手法として Goldstine と von Neumann により後に紹介されるだろうと述べられているが, それに該当するような彼らが著者である論文は見つけられなかった.

\overline{main.tex 2011/02/28(河合佑太)}

ここで,

α = AΔt/(Δx)^{2}, b = βΔx^{2}Δy, c = Δt/(2ΔxΔy),
h = 4RHγΔt/(f_{0}c_{p}).

である.

よって, 水平拡散項 A∇^{2}η は完全に陰的に評価される一方で, 地表の摩擦項 −kζ_{4} は, 陽的・陰的の両方で評価される.

3.2 数値計算手法

(3.1.11), (3.1.12) によって η_{ij τ+1} を求める作業は 2 つのステップに分けられる. 最初に方程式右辺 χ_{ij} の計算を行い, 続いて 2 本の方程式を連立して ζ_{s ij τ+1} を求める. つまり,

[αℒ_{ij} − μ_{s}] η_{s ij τ+1} = −χ_{s ij}  (3.2.13)

を解くことになる. ここで, s = 1 に対して μ_{s} = 1, s = 3 に対して μ_{s} = 1 + \frac{3}{2}kΔt である. なお, この方程式は擾乱場と平均場に対する方程式として, 2 本の方程式に書き直される.

[αℒ_{ij} − μ_{s}] ηs ij τ+1^{′} = −χs ij^{′},  (3.2.14)
αε^{2}(η̄_{j+1} − 2η̄_{j} + η̄_{j−1})_{s τ_{+}1} − μ_{s}η̄_{s j τ+1} = −χ̄_{s j}.  (3.2.15)

(3.2.14) は緩和法 (ガウス・ザイデル法) を用いて解かれる. このとき, 境界条件として j = 0, L において η′ = 0 を, 周期境界条件として

η′(i = 0) = η′(i = I)  (3.2.16)

を採用する. (3.2.15) は, (3.1.5) や (3.1.6) で使われた同様の方法によって正確に計算される. このとき, 境界条件も同様に

ψ̄(j = 0) = ψ̄(j = 1),
ψ̄(j = J) = ψ̄(j = J − 1)  (3.2.17)

を適用する. 最後に, η′ と η̄ を加えることによって η_{ij τ+1} を得る.

運動の「系統的な記述」は, ポテンシャル渦度の擾乱場によって行われる. すなわち, ポテンシャル渦度の擾乱場が (3.1.11) と (3.1.12) によって時間的に外挿され, その後にそれぞれの時間ステップにおいて (3.1.3) と (3.1.4) から流線関数の擾乱場が決定される.

(3.1.3), (3.1.4), (3.2.14) の緩和法に Liebman プロセスが用いられる. (3.1.3), (3.1.4) 中の ψij τ+1^{′} のための最初の推測値 < ψijτ+1^{′} > には

< ψij τ+1^{′} >= 2ψij τ^{′} − ψij τ−1^{′}  (3.2.18)

が用いられる. このため ψτ−1^{′} の値を保持する必要があるが, Liebman プロセスにおいて必要とされる繰り返し数が減少するために計算速度が向上する. ψ′ に対す

る緩和プロセスは, |ψν+1^{′} − ψν^{′}|_{1 ij} + |ψν+1^{′} − ψν^{′}|_{3 ij} が予報領域全体で\frac{3}{8} × 10^{5} m^{2}s^{−1}

以下になるまで続けられる. ここで, ν は繰り返し回数を表す添え字である.

(3.2.14) の緩和法による解法は, α/μ が小さな値であるために急速に収束する (最適化された過緩和プロセスは約 10^{−3} の固有値を持つ). ここでは, 最初の推測値として < η′ >= μ^{−1}χ′ をとった. η′ に対する緩和プロセスは, ψ′ と同様に 600 m^{2}s^{−1} になるまで続けられる.

3.3 エネルギーパラメータと子午面循環の数値的な評価法

有限差分形式の予報方程式 (3.1.3)∼(3.1.6), (3.1.11)∼ (3.1.12) は正確なエネルギー積分を与えないので, 一般にエネルギー関数を有限差分で評価する際には任意性がある. よって 2 回分のタイムステップの ψ の値を使って関数を評価するためのより良い方法は自動的には決まらない. ここでは, 今回使用したエネルギー関数の計算法を簡潔に記すだけにする.

簡単化のために, 総和について以下の縮約記号を導入する.

\begin{matrix}X S_{ij} = I^{−1} ∑_{i=0}^{I−1} S_{ij}, \\ Y S_{ij} = J^{−1} ∑_{j=0}^{J−1} S_{ij}. \end{matrix}

24 時間ごとに, レベル 1 とレベル 3 の最新のタイムレベル 2 つ分の流線関数場 ψ_{ij τ}, ψ_{ij τ−1} を記録するものとする. このとき, これらの平均と差分をそれぞれ次のように定義する.

\begin{matrix}δψ_{ij} = ψ_{ij τ} − ψ_{ij τ−1}, \\ ψ_{ij}^{∗} = \frac{1}{2}(ψ_{ij τ} + ψ_{ij}). \end{matrix}

これらの量は, エネルギーパラメータや子午面循環の計算に使われる. 以下, 平均を表す添え字の ∗ を省略することにする.

ω_{2 ij} は (2.2.78) の有限差分形式から以下のように決定される.

ω_{2 ij} = \frac{λ^{2}p_{2}}{f_{0}} [\frac{(δψ_{1} − δψ_{3})_{ij}}{δt} − \frac{𝒥_{ij}(ψ_{1}, ψ_{3})}{4ΔxΔy} + \frac{2RH(2j − J)}{f_{0}c_{p}J} − \frac{A}{Δx^{2}} ℒ (ψ_{1} − ψ_{3})] .

(3.3.19)

子午面循環は次の方程式によって決定される.

V̄_{1 j+1/2} = V̄_{1 j−1/2} − \frac{Δy}{p_{2}} ω̄_{2j}, j = 1, 2, ..., J − 1.  (3.3.20)

ここで, j = 1/2 において V̄_{1} = 0 であり, ω̄_{2j} は X ω_{2 ij} によって与えられる. 正確な有限差分の解は j = J − 1/2 において V̄_{1} と V̄_{3} をゼロにするので, j = J − 1/2 において V̄_{1} の小ささから, 計算における丸め誤差の影響を確認することができる.

運動エネルギーは次の式から計算される.

\begin{matrix}K′ = & \frac{10}{2JΔx^{2}} ∑_{j=0}^{J−1} X\{[(ψj+1^{′} − ψi^{′})_{1j}^{2} + (ψi+1^{′} − ψi^{′})_{3j}^{2}] + ε^{2}[(ψj+1^{′} − ψj^{′})_{1i}^{2} + (ψj+1^{′} − ψj^{′})_{3i}^{2}]\}, \\ K̄ = & \frac{10}{2Δy^{2}} Y [(ψ̄_{j+1} − ψ̄_{j})_{1}^{2} + (ψ̄_{j+1} − ψ̄_{j})_{3}^{2}]. \end{matrix} ; (3.3.21)

そして, ポテンシャルエネルギーは次の式から計算される.

\begin{matrix}P′ = & \frac{10}{2} Y X(ψ1^{′} − ψ3^{′})_{ij}^{2}, \\ P̄ = & \frac{10}{2} Y (ψ̄_{1} − ψ̄_{3})_{j}^{2}. \end{matrix}  (3.3.22)

また, 11 種類のエネルギー変換は次の式によって計算される.

\{Q̄ · P̄\} = −\frac{2RHλ^{2}l}{f_{0}c_{p}} Y \frac{2j − J}{J} (ψ̄_{1} − ψ̄_{3})_{j},
\{P̄ · P′\} = −\frac{λ^{2}l}{4ΔxΔy} Y [(ψ̄_{1} − ψ̄_{3})_{j} X𝒥 (ψ1^{′}, ψ3^{′})_{ij}],
\{P′ · K′\} = −\frac{f_{0}l}{p_{2}} Y Xω2 ij^{′}(ψ1^{′} − ψ3^{′})_{ij},
\{K′ · K̄\} = \frac{l}{4Δx^{3}Δy} Y ([(ψ̄_{j−1} − ψ̄_{j+1}X(ψi+1^{′} − ψi−1^{′})_{j} ℒ ψij^{′}]_{1},

+ [(ψ̄_{j−1} − ψ̄_{j+1})X(ψi+1^{′} − ψi−1^{′})_{j} ℒ ψij^{′}]_{3}),

\{P̄ · K̄\} = −\frac{f_{0}l}{p_{2}} Y ω̄_{2j}(ψ̄_{1} − ψ̄_{3})_{j},
{K̄ · A} = Al Y (ζ̄_{1}^{2} + ζ̄_{3}^{2})_{j},
{K′ · A} = Al Y X[(ζ1 ij^{′})^{2} + (ζ3 ij^{′})^{2}],
\{P̄ · A\} = \frac{λ^{2}Al}{Δy^{2}}Y [(ψ̄_{1} − ψ̄_{3})_{j+1} − (ψ̄_{1} − ψ̄_{3})_{j}]^{2},
\{P′ · A\} = \frac{λ^{2}Al}{JΔx^{2}} ∑_{j=0}^{J−1}\{[(ψ1^{′} − ψ3^{′})_{i+1} − (ψ1^{′} − ψ3^{′})_{i}]_{j}^{2} + ε^{2}[(ψ1^{′} − ψ3^{′})_{j+1} − (ψ1^{′} − ψ3^{′})_{j}]_{i}^{2}\},
{K̄ · k} = −kl Y ζ̄_{4j}ψ̄_{3j},
{K′ · k} = −kl Y Xζ4 ij^{′}ψ3 ij^{′}.

ここで, l = 10 days = 8.64 × 10^{5} s, ζ_{ij} = Δx^{−2} ℒ_{ij}ψ である. 上の 4 種類のエネルギーと 11 種類のエネルギー変換の単位は, 後の再現結果の表で用いられるエネルギー単位と同じものである. なお, {K′ · K̄} の定式化には, u, v が流線関数によって与えられるときに有効である, \overline{vζ} = −∂(\overline{uv})/∂y という関係を用いている.

第4章数値実験の手順

4.1 数値実験の概要

回転系に対して静止している大気から数値実験を始めるならば, 最初は物理的な原理から大気大循環を差し引くことができると思われる. このとき, 実際に観測される大気の運動は, 物理法則の自然な結果としてどのように現れるか説明されるだろう. 今の場合の物理法則は, 簡単化された準地衡風モデルに対する方程式によって表現される. この物理法則には経験的事実も組み込まれているが, それは物理的な理論に基づかれている. よって, 物理法則の簡単化のために応用は限られるが, 作成した数値予報モデルは大気大循環の全体的な特徴について部分的な説明を与えるだろう.

今, 理想的な数値実験の初期条件は静止大気ということになる. 小規模な擾乱運動を取り除けば, 長期間の平均場に対する数値予報を行える. 最初に南北に変化する加熱および冷却は, 東西に一様な南北温度勾配を生成し始める. 南北温度勾配が特定の臨界値に達した後, ある特徴的な波長を持った大規模な渦がおそらく発達し始める (Eady, 1949). そして, 今度は発達した大規模な渦が非線形相互作用によって平均帯状流を変化させるだろう. このとき, エネルギーの入力が摩擦によるエネルギーの散逸によって釣り合わされるようなある種の準定常状態に完全に到達することが望まれる. そのとき典型的な流れのパターンは, 現実の大気に見られる流れのパターンと類似点を持つであろう.

線形摂動理論を用いれば, 初期の小さな擾乱が増幅し始める前に存在しなければならない南北温度勾配および卓越する波の規模の両方を推定できる. W を 5000 km にとれば, モデルの南北幅としておよそ地球の極と赤道間の距離をとったことになる. この場合, 今回の 2 層モデルにおける南北の臨界温度差は約 56 度であり, もっとも不安定な波長は 6000 km である (Phillips, 1954).

これらの値は, 対流圏における極と赤道間の年平均の温度差 (500 hPa では約 34 度) や実際に温帯で発生する擾乱の波長に比べたとき, あまり現実的でない. この差異は, 大気中の傾圧不安定を表現するために 2 層モデルを用いたことによる部分的な失敗であろう. このことは Kuo(1953) によって理論的に示されており, また 2 層地衡風モデルと 3 層地衡風モデルを使って得られた数値予報の違いによっても示唆されている.

4.2 数値実験の手順

Phillips の数値実験は最初に ψ_{1} = ψ_{3} = 0 とし, 静止大気から始められる. 渦なし状態においては, (2.2.70) のように選ばれた加熱率のために南北端で 1 日につき 0.23 度の加熱と冷却が生じる. よって, 傾圧不安定波を発生させるのに十分な南北方向の温度差を生成するには約 4 か月ほど必要となる. (この必要時間は, H を決める時に潜熱の輸送を考えていないことや系の幾何的な違いのために, 現実大気の平均的な放射の効果に比べると長い.) したがって, 最初の 130 日間の予報は東西方向に非依存な運動 (渦なし) に対して行われる. これにより Phillips は計算安定性を損なうことなく, 1 タイムステップあたりの時間の刻み幅を 1 日にとって時間積分を行うことができた. この数値積分の結果, 東西平均された温度場の領域平均温度からのずれと東西平均風の分布が得られる.

続いて, (x, y) 平面上で変動するがレベル 1 とレベル 3 で同じであるような微小でランダムな擾乱を系に挿入する. このとき擾乱として与える ψ′ の格子点上の値は擬似乱数を使って生成される. 渦を挿入した 130 日後からの予報は, 2 時間かそれ以下の時間の刻み幅を用いて (計算の安定性に依存する), 完全な方程式系を使って行われる. なお, 現在の計算機の能力で同様の計算を行う場合には, 1 タイムステップあたりの時間の刻み幅を常に 30 分に取っても計算時間は問題にならないので, 本研究では常に 30 分とした. ランダムな擾乱を挿入してから 1 か月ぐらいで, 計算不安定により K̄ や K′ が急激に増加し, 数値積分が続行不可能となる.

以後, 渦なしで東西平均場に対してのみ行う最初の 130 日間の数値実験を「数値実験 1」と呼ぶ. その後, 渦を系に挿入し完全な予報方程式系を用いる 130 日目以降の数値実験を「数値実験 2」と呼ぶことにする.

4.3 パラメータの設定

当時, 格子点間隔の大きさは用いられた電子計算機の能力で決まった. また, 付随する技術的問題から最終的に水平格子点数は x 方向に 16, y 方向に 17 として選ばれた. 本研究では, 当時使われた格子点間隔をそのまま使うことにした. x 方向の差分幅 Δx として大きな値を選べば距離 L = 16Δx の中にいくつかの渦をおそらく生み出し得るが, 相対的に大きな離散化誤差を生み出すだろう. 一方で, Δx を小さく選べば波長 L の大きな渦がひとつできるだけだが, 離散化誤差はより小さくなるだろう. ここでは後者を選択することにして L = 6000 km にとれば, Δx = 375 km となる.

最終的に, パラメータおよび物理定数は以下の値を用いた.

H = 2 × 10^{−3} Jkg^{−1}s^{−1}
λ^{2} = f_{0}^{2}θ^{2}[(θ_{1} − θ_{3})(Φ_{1} − Φ_{3})]^{−1} = 1.5 × 10^{−12} m^{−2}
Δx = 375 km
Δy = 625 km
R = 287 Jkg^{−1}K^{−1}
c_{p} = 1004 Jkg^{−1}K^{−1}
A = 10^{5} m^{2}s^{−1}
k = 4 × 10^{−6} s^{−1}
f_{0} = 10^{−4} s^{−1}
β = 1.6 × 10^{−11} s^{−1}m^{−1}
p_{2} = p_{4}/2 = 500 hPa
I = J = 16

ここで, λ^{2} の計算において

\begin{aligned}Φ_{1} − Φ_{3} & = 77, 499 m^{2}s^{−2} \\ θ_{2}/(θ_{1} − θ_{3}) & = 11.625 \end{aligned}

とした.

第5章数値実験の結果

5.1 数値実験 1 (渦なし状態における帯状流の発達) の結果

予報開始から 130 日目 (数値実験 1 終了後) の 250 hPa 面, 500 hPa 面, 1000 hPa 面における平均東西速度場 ū_{1}, ū_{3}, ū_{4}, 250 hPa 面における東西平均された相対渦度場 ζ̄_{1}, 500 hPa 面における帯状平均した温度場の領域平均からの偏差 T̄_{2}, 平均南北風 V̄_{1} の予報結果をそれぞれ図 5.1.1, 図 5.1.2, 図 5.1.3, 図 5.1.4, 図 5.1.5, に示す. ここで, j は緯度方向の座標を示し, 8(y + W)/W に等しい. なお, Phillips(1956) の数値実験の結果と本研究による再現結果を対比するために, 前者の結果を破線, 後者の結果を実線で同じグラフに重ねている.

図 5.1.1: 数値実験 1 終了時 (130 日目) の 250 hPa 面における帯状平均した東西速度場. 単位は m s^{−1} である.

図 5.1.2: 数値実験 1 終了時 (130 日目) の 500 hPa 面における帯状平均した東西速度場. 単位は m s^{−1} である.

図 5.1.3: 数値実験 1 終了時 (130 日目) の 1000 hPa 面における帯状平均した東西速度場. 単位は m s^{−1} である.

図 5.1.4: 数値実験 1 終了時 (130 日目) の 250 hPa 面における帯状平均した渦度場. 単位は s^{−1} である.

図 5.1.5: 数値実験 1 終了時 (130 日目) の 500 hPa 面における帯状平均した温度場の領域平均からの偏差. 単位は K である.

図 5.1.6: 数値実験 1 終了時 (130 日目) の 250 hPa 面における帯状平均した南北速度場. 単位は m s^{−1} である.

次に, 数値実験 1 終了時 (130 日目) の帯状流の運動エネルギー K̄ と東西平均場に対するポテンシャルエネルギー P̄ およびそのときの平均流と関連する 5 種類のエネルギー変換の値を表 5.1.1 に示す. 先ほどと同様に対比のために Phillips (1956) の結果も一緒に示した. ここで, K̄ の単位は, レベル 1 と 3 における全ての格子点上で速度が 1 m/s であるときに K̄ が 10 となるようにとられている. P̄ の単位についても, それと同じものが使われている. エネルギー変換の単位は, 24 時間で K̄(または P̄) が 1 単位変化するようなエネルギー変換を 1 単位としている.

表 5.1.1: 数値実験 1 終了時のエネルギーパラメータ

エネルギーパラメータ再現結果 Phillips(1956) の結果
K̄  4198  4265
P̄  23996  24368
{Q̄ · P̄}  437  440
{P̄ · K̄}  45  41
{K̄ · A}  9  9
{K̄ · k}  -27  -27
{P̄ · A}  47  48

5.2 数値実験 2 (渦あり状態における流れの発達) の結果

数値実験開始から 3 日ごとの 1000 hPa 面における流線関数場 (ψ_{1} と ψ_{3} によって線形的に外挿される) と 500 hPa 面における全領域平均温度からの温度偏差の分布 (ψ_{1} − ψ_{3} に比例する) を図 5.2.7, 図 5.2.8 に示す.

図 5.2.7: 数値実験 2 開始後 (数値実験 1 開始から 131 日目) から 16 日目までの 3 日おきの 1000 hPa 面における流線関数場と 500 hPa 面における全領域平均温度からの温度偏差の分布. 前者はトーン, 後者はコンターで表現し, 単位はそれぞれ m, m^{2} s^{−1} である.

図 5.2.8: 数値実験 2 の 19 日目から 28 日目までの 3 日おきの 1000 hPa 面における流線関数場と 500 hPa 面における全領域平均温度からのずれの温度分布. 前者はトーン, 後者はコンターで表現し, 単位はそれぞれ m, m^{2} s^{−1} である.

数値実験 2 の 22 日目における鉛直運動の場 (鉛直 p 速度) を, 図 5.2.9 において 1000 hPa 面における流線関数場と重ねて示す.

図 5.2.9: 数値実験 2 開始から 22 日目における 1000 hPa 面の流線関数場と鉛直速度場 (dp/dt) の分布. 前者はコンター, 後者はトーンで表現し, 単位はそれぞれ m^{2} s^{−1}, Pa s^{−1} である. 青色の領域が上昇流, 赤色の領域が下降流に対応する.

図5.2.10は, 数値実験 2 における平均子午面循環V̄_{1} の様子を示している. (2.4.101) の定義から, V̄_{1} とは大気の上側半分の平均的な南北速度であり, 下側半分の平均的な南部速度 V̄_{3} はその逆で −V̄_{1} となる.

図 5.2.10: 数値実験 2 における大気上半分の平均子午面循環の緯度 (j) 方向の変化と時間変化. 単位は m s^{−1} である.

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 57} 図 5.2.11 は, 500 hPa 面における東西平均温度 T̄_{2} の時間的な変化を, 数値実験 2

開始からの温度変化の積算として示している.

図 5.2.11: 数値実験 2 の 500 hPa 面における東西平均温度の緯度 (j) 方向の変化と時間変化. 数値実験 1 終了時からの温度変化を積算している. 単位は K である. \overline{main.tex 2011/02/28(河合佑太)}

250 hPa と 1000 hPa における平均帯状流 ū_{1} と ū_{4} の時間的な変化は, 図 5.2.12 と図 5.2.13 で示される.

表 5.2.2 は, 予報サイクルにおいて 24 時間ごとに得られるさまざまなエネルギーパラメータの値を示している. 単位は数値実験 1 の結果の際に用いたものと同じである. 図の最後の行には, 数値実験 2 の 5 日目から 25 日目までの平均値も記した.

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 59}

図 5.2.12: 数値実験 2 における 250 hPa 面における平均帯状流 ū_{1} の緯度 (j) 方向の変化と時間変化. 単位は m s^{−1} である. 図 5.2.13: 数値実験 2 における 1000 hPa 面における平均帯状流 ū_{4} の緯度 (j) 方向の変化と時間変化. 単位は m s^{−1} である. \overline{main.tex 2011/02/28(河合佑太)} 表 5.2.2: 数値実験 2 におけるエネルギーとエネルギー変換の値

Days  K′  K̄  P′  P̄  E {Q̄ · P̄} {P̄ · P′} {P′ · K′} {K′ · K̄} {P̄ · K̄} {K̄ · k} {K′ · k} {K̄ · A} {K′ · A} {P̄ · A} {P′ · A}
0  755 4263  0 24353 29371  440  0  0  -166  45  -27  267  9  334  48  0
1  335 4324  40 24697 29396  443  -3  21  -16  59  -27  97  9  135  48  7
2  203 4384  29 25041 29658  446  0  -9  -28  52  -28  27  10  70  49  4
3  151 4449  26 25387 30012  449  20  9  6  57  -28  7  10  43  50  2
4  137 4511  39 25719 30407  452  29  16  25  51  -28  3  10  28  50  2
5  132 4574  71 26029 30806  455  60  41  29  44  -28  4  10  19  51  2
6  192 4641 122 26302 31258  457  112  85  40  34  -27  7  10  15  52  3
7  256 4717 200 26517 31690  459  198  159  96  21  -24  14  10  15  52  5
8  385 4797 321 26634 32137  460  328  276  91  8  -18  26  10  17  52  5
9  695 4881 521 26583 32680  460  544  463  422  -4  -6  44  10  24  52  12
10  912 4991 821 26274 32997  457  840  757  448  -31  18  71  11  35  52  14
11 1589 5190 1140 25700 33619  452  1083  1102  1079  -83  70  99  13  52  51  27
12 1637 5550 1242 25046 33474  446  1002  1303  1109  -121  155  117  20  76  50  26
13 1562 6029 884 24741 33217  443  502  1008  1641  -125  219  107  35  91  51  26
14  975 6272 354 24984 32586  444  67  415  833  -90  195  75  51  89  54  15
15  822 6129 200 25322 32472  447  200  179  362  -30  134  47  60  84  56  10
16 1031 5942 440 25362 32773  447  543  337  -50  -36  99  42  69  77  56  13
17 1459 5902 881 25037 33279  445  880  701  788  -74  92  62  86  89  56  24
18 2588 6140 1376 24345 34449  439  1326  1088  161  -131  112  96  123  108  59  31
19 2547 6501 1802 23398 34248  430  1391  1332  2073  -155  168  118  148  129  60  39
20 3800 6768 1986 22489 35043  421  1281  1402  1134  -205  262  113  130  166  58  50
21 2540 7245 1621 22044 33450  416  696  1099  1642  -238  354  70  117  168  59  36
22 2517 7534 1201 22028 33280  415  294  504  2848  -130  339  34  129  150  62  35
23 2911 7170 1099 22216 33396  416  342  -79  -1363  78  180  80  144  152  67  29
24 5008 5993 2080 21566 34646  413  1602  973  1811  -40  88  172  143  215  59  53
25 6901 5070 3649 19469 35089  394  2973  2076  1715  -77  130  251  80  294  48  63
26 10373 4962 5023 16808 37166  365  2752  2443  2189  -201  265  222  45  333  40  114
27 7675 6153 4844 15187 33859  342  1343  2425  5307  -658  520  114  48  382  35  99
28 5635 8179 2584 15494 31892  337  200  1541  4477  -612  714  46  73  404  41  85
29 4808 7948 1527 16233 30516  342  -350  -531  -1401  -101  502  80  110  461  51  55
30 8268 6925 2313 16764 34270  349  210  -92  2416  -181  282  116  118  579  51  85
31 9070 5810 3662 15645 34187  343  2098  933  -1262  99  182  283  133  705  49  89
5-25(mean) 1927 5811 1048 24385 33171  439  775  725  805  -66  120  79  67  98  55  25

\overline{main.tex 2011/02/28(河合佑太)}

第6章数値実験の結果の考察

6.1 数値実験 1 (渦なし状態における帯状流の発達) の考察

静止大気から運動の発生までについて

数値実験 1 は, 回転系に対して静止している大気 (ψ_{1} = ψ_{3} = 0) から始められる. 支配方程式 (2.2.80) と (2.2.81) によれば, 非断熱加熱は東西平均されたポテンシャル渦度場 q̄_{1} と q̄_{3} に変化をもたらす. その結果, (2.2.79) で表わされるように東西平均された流線関数場 ψ̄_{1}, ψ̄_{3} が変化し運動が生じる. ψ̄ や q̄がゼロでなくなった後は, 運動は非断熱加熱の他に摩擦項 A ∂^{2}q̄/∂y^{2} と −kζ̄_{4} の影響を受ける.

数値実験 1 終了時の帯状平均場について

数値実験 1 の開始から 130 日目の運動の様子を図 5.1.1, 図 5.1.2, 図 5.1.3, 図 5.1.5 を使って考察することにする. 以下対比のために, Phillips(1956) の数値を括弧の中に示す (平均子午面循環 V̄_{1} を除いて, 今回の数値実験 1 の結果は Phillips(1956) と非常によく一致していることが分かる. ). まず温度分布はとても単純であり, y = W の - 30.1(- 30.1) 度から y = −W の 30.1(30.1) 度までほぼ線形的に増加する. 東西風の分布もまた非常に規則的であり, レベル 1 では広範囲で西風であり, 36 m/s の極大値をもつ. 一方で, 地上のレベル 4 では全域にわたって弱い東風が吹いている. 子午面循環 V̄ は, 速度場の極大値が 3.4(3.0) cm/s であるような 1 つの直接循環セルで構成される. 大気の平均的な東西方向の運動量 (ū_{2} ≡ (ū_{1} + ū_{3})/2 によって与えられる) は正であり, 地表面摩擦は地上の弱い東風に影響を与える. 現実大気の帯状流のジェットのような特徴は, ū_{1} の分布には見られない. 東西平均された相対渦度 ζ̄_{1} = −∂ū_{1}/∂y の分布は完全に単調であり, 正しいジェットの特徴である ζ̄_{1} が ū の極大域の北側で正の極大, 南側で負の極小が存在しない.

数値実験 1 終了時のエネルギー変換について

次に表 5.1.1 から, 数値実験 1 におけるエネルギー変換を考察する. このエネルギー変換の結果は, エネルギーが子午面での弱い直接循環セルの作用によって P̄ から K̄ へと移動することを示している. {K̄ · k} は負であるが, このことは地表面摩擦を要因としてレベル 1 とレベル 3 の帯状流の運動エネルギーの増加を許可することになり, 一見おかしさを感じるだろう (熱力学第 2 法則に反するように見える). このモデルの地表摩擦の表現方法では, 摩擦層上端で弱い東風であるときに摩擦層全体において正味では赤道方向への流れを生み出す. この流れを補填するために, 摩擦層の上側において正味では極方向への流れが存在しなければならない. この極方向への流れは, コリオリ力によって帯状流の加速度の変化に寄与する. 摩擦層の上の平均帯状流は西風なので, 極方向の流れはそれぞれのレベルで帯状流の運動エネルギーを増加させる. なお, 数値実験 2 では地表の東西風の分布はより現実的になるために, {K̄ · k} は正となる.

6.2 数値実験 2 (渦あり状態における流れの発達) の考察

渦について

数値実験 1 の終了後 (130 日後) に存在している平均流中に, 微小で「ランダム」な擾乱を加える. この初期擾乱は運動エネルギー K′ として 768 単位分のエネルギーを持ち, 速度の二乗平均平方根は 8.8 m/s に対応する. 初期擾乱の流線関数場はレベル 1 と 3 で同じとしている (ψ1^{′}(x, y) = ψ3^{′}(x, y)) ので, P′ は数値実験 2 の最初においてゼロとなる. 構成される流れのパターンが時間とともに発達する様子は, 完全な方程式系 (2.2.79)∼ (2.2.81) を使って計算される.

表 5.2.2 から擾乱運動の運動エネルギー K′ の変化の様子を見てみると, 数値実験 2 開始から 5 日間は減少しているものの, 以後 25 日間はほぼ一様に増加している (ただし, 25 日間中の所々で一時的に減少しているところもある). この図には示していないが, 31 日目以後 K′ は時間的に急激に増加する.

この K′ の初期の減少は, 初期の {K′ · A} と {K′ · k} の相対的な大きな値に示されるように, 小規模な水平渦粘性と地表面摩擦を原因とする. この水平渦拡散のプロセスは, 不規則な分布を持つ初期の擾乱場 ψ′ の特性のために重要である. 一方で, 地表面摩擦の効果は, ψ1^{′} と ψ3^{′} の水平分布が初期において同じであるために (順圧的であるために) 小規模な渦をスピンダウンさせる.

数値実験 2 の 5 日目以降からは, 実際の天気図に見られる大規模な渦の特性を持つようなとても規則的な波 (傾圧波動) が発達する. 図 5.2.7, 図 5.2.8 から読み取れる大きのな渦の特徴は次のようになる.

• 波は約 1500 km/day の速度で東へと移動する. この規則性は, 図に示していない日においても当てはまる.

• 波は温暖型低気圧として始まる (例えば 10 日目) が, 舌状の冷気の南下は段々と深まりつつある地表面の低気圧の中心に追いつく. また, 舌状の暖気の北上は段々と地上の高気圧の中心へと移動する. 結果, 数値実験 2 の最後の方では閉塞低気圧のような特徴が見られる.

• おもな地表のトラフやリッジは, 北側半分の領域では北西に, 南側半分では南西に向かって傾いている.

• 主な温度勾配は 1000 hPa におけるトラフの西側で発生する構造が見られ, 実際の寒冷前線や温暖前線に見られる特徴に当てはまる.

図 5.2.9 から鉛直運動について考察する. 鉛直速度 (総観スケールの擾乱に対しては w ≃ −ω/(ρg) となる) は数 cm のオーダであり, 総観スケールの擾乱の鉛直速度として現実的である. 流線関数場に対する鉛直運動の分布も現実的であり, 低気圧の領域における上昇流や高気圧の中心における下降流に加え, トラフの前後で上昇流が下降流へと急激に変化する特徴がみられる.

平均子午面循環について

図 5.2.10 から, 数値実験 2 における平均子午面循環 V̄_{1} の 30 日間の変化について考察する. 数値実験 2 開始から最初の間は, 数値実験 1 終了時の大変弱くて広範囲に及ぶ直接循環セルの特徴のみが見られる. しかしながら, 渦が発達するにつれて, 中緯度の間接循環のセルとその南側と北側でいくらか弱い直接循環のセルを伴う, 3 つのセル構造を持つ循環が明白に見られるようになる. これは, 2 層モデルにおける不安定な傾圧波の特徴的な性質である. 数値実験 2 終盤からは, 計算不安定により V̄ の場が発散傾向になる.

温度の帯状平均場の時間的な変化と系の熱収支について

図 5.2.11 から数値実験 2 における温度の帯状平均場の 30 日間の変化について考察する. 数値実験 2 開始から 5 日間における y = −W と y = W 間の温度差は, 数値実験 1 の渦なし状態の場合と同じ割合 (1 日につき約 0.4 度づつ) で増加し続けている. 擾乱が発達するにつれて, その擾乱は v′ と T2^{′} の間の相関のために潜熱を極方向へと輸送し始める. この渦による熱輸送のために非断熱加熱による南北間の温度差の増加は止み, 数値実験 2 が終わるころには数値実験 2 開始時の南北温度差を下回るようになる. また, 数値実験 2 が後半に近づくにつれて, 中緯度帯に温度勾配が集中する一方で, 領域の北と南では温度勾配は小さくなる. このような特徴が現れる前に, 南北の熱交換のプロセスは j = 3 付近に相対的に冷たい領域や j = 12 付近に相対的に暖かい領域を生み出す. これは, 現実の天気図に現われる主な西風の流れの北側と南側で見られるような, カットオフされた温暖型高気圧や寒冷型低気圧を想像させる.

この数値予報モデルに対する熱力学的エネルギー保存の式 (2.2.78) の帯状平均の形式は,

\frac{∂T̄_{2}}{∂t} = −\frac{∂\overline{(v1^{′}T2^{′}})}{∂y} − \frac{2H}{c_{p}} \frac{y}{W} + \frac{f_{0}^{2}}{Rλ^{2}p_{2}} ω̄_{2} + A\frac{∂^{2}T̄_{2}}{∂y^{2}}  (6.2.1)

となる. 数値実験 2 の 10 日目から 20 日目の 10 日間のこの方程式の 5 項の平均の大きさと y 方向の分布は, 図 6.2.1 のようになる. なお, y = (j − 8)Δy である. 非断熱加熱の寄与 (−2Hy/(c_{p}W)) はもちろん低緯度 (j < 8) で加熱, 高緯度 (j > 8) で冷却となっている. 渦による熱輸送の収束 (−∂(\overline{v1^{′}T2^{′}})/∂y) による寄与は, 非断熱加熱の寄与とは逆に働いており, 特に領域の中心で南北温度勾配を破壊する傾向がある. 他方, 平均子午面循環 (f^{2}ω̄_{2}/(Rλ^{2}p_{2})) によって生み出される温度変化は, 南北境界の近くでは非断熱加熱・冷却の効果による温度変化を打ち消している. しかし, 中緯度における間接循環のために領域の中心付近の −∂T̄_{2}/∂y の増加に大きく寄与している. 結果的に正味の温度変化 ∂T̄_{2}/∂t としては, 領域の北側と南側では南北温度勾配を減少させる傾向があり, 中心付近では南北温度勾配を増加させる傾向がある. 図にも見られるように水平拡散による正味の温度変化の寄与は小さいので, 今回の議論では無視した.

図 6.2.1 で用いられた数値は完全な平衡状態におけるものではないが, 大気中の ∂T̄_{2}/∂y が最大となる緯度帯はある程度力学的に決定されるものであり, 放射の効果の南北勾配によってはほとんど決定されない.

平均帯状流と運動量収支について

250 hPa と 1000 hPa における平均帯状流 ū_{1} と ū_{4} の数値実験 2 の 30 日間における変化について, 図 5.2.12 と図 5.2.13 を用いて考察する. これらの図において興味深い特徴は次のようになる.

• 250 hPa 面では領域中心において強い西風ジェットが形成される.

• 1000 hPa 面では最初一様で弱い東風であるが, 中緯度で西風, 領域の北側と南側で東風のパターンへと変化する.

図 6.2.1: 数値実験 2 の 10-20 日目間の熱力学的エネルギー保存の式 (6.2.1) の各項の平均的な大きさを示したグラフ. 紺色の実線は東西平均されたレベル 2 の正味の温度変化率 (左辺), 赤色の実線は非断熱加熱・冷却の効果による温度変化率, 黄緑色の実線は断熱的な加熱・冷却の効果による温度変化, 紫色の実線は渦による顕熱のフラックスの収束による温度変化, 水色の実線は熱拡散による温度変化を示している. 単位は K / day.

これらの変化のプロセスは, 大気の上側半分と下側半分の帯状流の運動量収支を調べることによって理解される.

ū_{1} の時間変化は, (2.2.42) を帯状平均することで次のように得られる.

\frac{∂ū_{1}}{∂t} = −\frac{∂(\overline{u1^{′}v1^{′}})}{∂y} + f_{0}V̄_{1} + A\frac{∂^{2}ū_{1}}{∂y^{2}} .  (6.2.2)

この式の各項の大きさを 10-20 日間に渡って平均し, y(j) の関数としてグラフに表したのが, 図 6.2.2 である.

数値実験 2 の期間中 ū_{1} は中緯度で増加し高緯度と低緯度では減少するが, この変化は明らかにコリオリ項 f_{0}V̄_{1} と大規模な渦による応力の項 −∂(\overline{u1^{′}v1^{′}})/∂y によって主に決定される. 特に後者の方が ū_{1} の変化に大きく寄与する. 子午面循環による寄与は, 一般にはレベル 1 の渦の項の効果を減少させるように観測される ū_{1} の変化と逆である. よって, 結果的に中緯度帯に帯状流の極大を生み出す.

図 6.2.2: 数値実験 2 の 10-20 日目間の大気上半分の運動量収支 (6.2.2) の各項の平均的な大きさを示したグラフ. 紺色の実線は東西平均されたレベル 1 の正味の運動量の時間変化率 (左辺), 黄緑色の実線は子午面循環のコリオリ力による運動量の変化, 紫色の実線は大規模な渦の応力の効果による運動量の変化, 水色の実線は小規模な渦に対する水平渦拡散による運動量の変化を示している. 単位は m s^{−1} day^{−1}.

大気の下側半分に対する運動量収支も, 上と同様に運動方程式の帯状平均によって議論される.

\frac{∂ū_{3}}{∂t} = −\frac{∂(\overline{u3^{′}v3^{′}})}{∂y} + f_{0}V̄_{3} + A\frac{∂^{2}ū_{3}}{∂y^{2}} − kū_{4}.  (6.2.3)

この式の各項の大きさを 12-22 日間に渡って平均し y(j) の関数としてグラフに表したのが, 図 6.2.3 である. グラフを見ると主要な項はコリオリ力の項 f_{0}V̄_{3} と地表面摩擦の効果を表す項 −kū_{4} であることが分かる. これらの 2 つの項の大きさは大体同じぐらいであり, 符号は反対である.

したがって, 増幅する傾圧波に特有である暗黙的な子午面循環を表現する今回のモデルにおいて, 帯状流の運動量バランスの重要性がかなり明らかになった. 平均子午面循環の寄与は, 大気上半分の −∂(\overline{u′v′})/∂y の大きな値と大気下半分の地表面摩擦の効果の両方をバランスさせる傾向がある. よって, 中緯度における地表の西風とジェットの持続的な発達は, 大規模な渦による運動量輸送\overline{u′v′} と子午面循環 V̄ の複合的な結果である.

図 6.2.3: 数値実験 2 の 10-20 日目間の大気上半分の運動量収支 (6.2.3) の各項の平均的な大きさを示したグラフ. 紺色の実線は東西平均されたレベル 3 の正味の運動量の時間変化率 (左辺), 黄緑色の実線は子午面循環のコリオリ力による運動量の変化, 赤色の実線は地表面摩擦の効果による運動量の変化, 紫色の実線は大規模な渦の応力の効果による運動量の変化, 水色の実線は小規模な渦に対する水平渦拡散による運動量の変化を示している. 単位は m s^{−1} day^{−1}.

図 6.2.2 と図 6.2.3 から明らかなように, 子午面循環は帯状風の鉛直シアの大きな変化 (−∂(ū_{1}−ū_{3})/∂t) を減少させる. さもなければ, この帯状風の鉛直シアは大規模な渦の応力の項と地表面摩擦のみの結果であるだろう. この結果は, Fjørtoft(1951) によって導かれた.

系のエネルギー変換について

表 5.2.2 は, 予報サイクルにおいて 24 時間ごとに得られるさまざまなエネルギーパラメータの値を示している. 系へのエネルギーの入力値である {Q̄ · P̄} を最初に見てみると, 非断熱加熱と東西平均温度場 T̄_{2} との正の相関から想像されるように, 一様に正となっている. 数値実験 2 の 5-25 日間の {Q̄ · P̄} の平均は 439 単位である. 一方でその期間中のエネルギー損失は, 6 つのエネルギー損失を表すエネルギー変換の値を総和することで 444 単位と計算され, エネルギーの入力値とおおよそ同じ大きさである.

{Q̄ · P̄} の平均値である 439 単位は, 単位時間単位面積あたり約 2.6J のエネルギーを与えたことに対応する. これは, 実際の大気において想定されるエネルギー散逸の割合の約半分ほどである. その理由として H として選んだ値がやや小さいということに帰することもできるが, モデル自体が簡単化され過ぎているために, この平均値を実際の大気により近付けることは難しい.

2 層モデルのエネルギー論において, 大気で生じるエネルギー変換に対する機構の要点を述べた. つまり, 帯状流の運動エネルギー K̄ は擾乱からのエネルギーの供給によって維持される. また, その擾乱は平均場がもつポテンシャルエネルギー P̄ からエネルギーを受け取るといったものだった (Lorentz(1955) も参照). 表から読み取れるように, {P̄ · P′}, {P′ · K′}, {K′ · K̄} の平均値は, 全て正となっている.

{P̄ · K̄} によって与えられるような, P̄ と K̄ 間の直接的なエネルギー交換の符号は負が優勢であり, 特にこのとき中緯度における地表の西風の様子に変化をもたらす. これは, 北側と南側における弱い直接循環のセルや中緯度の相対的に大きな東西平均された南北温度勾配に比べて, その真ん中にできるの間接循環セルの強度がより強いという結果である. しかしながら, {P̄ · K̄} の大きさは小さい. そのため, 低緯度における真のハドレー循環は, {K′ · K̄} に比べると依然として小さいものの, {P̄ · K̄} を正にするのに十分なほど寄与するということが考えられる.

以上より, 図 2.3.5 において予想したエネルギー変換の方向は, {P̄ · K̄} が負であることを除いて, 計算によって裏付けられた. もちろん, 今回使用している方程式系はとても簡単化されているので, このダイアグラムは大気で生じる原理的なエネルギー変化の完全な表現であるとは言えない. しかし, 予報された流れのパターンのもっともらしさは, 今回のエネルギーダイアグラムがまずまず現実的であることを強く示唆している.

第7章結論

数値実験について

本研究で得られた数値計算結果は, Phillips(1956) で記されている実験結果と大変よく一致し, Phillips の数値実験の再現に成功したと言えるだろう. 今回作成した数値モデルの中で表現される大気大循環の構造は, Phillips(1956) が結論付けるように現実の大気におけるそれといくつか類似した側面を持っている. その側面とは,

• 地表における東西風の分布 (高緯度で東よりの風, 中緯度で西よりの風, 低緯度で東よりの風)

• 中緯度におけるジェットの存在

• 地球が獲得する正味の放射エネルギーの過不足を埋めるのに必要とされる, 大規模な渦による極方向への熱輸送の存在

である. また, これらの現象を生み出すモデルの中のエネルギー変換プロセスは, 大気の原理的なエネルギー変換と定性的には一致する.

方程式系を簡単化したことによる失敗もいくつか見られる. 例えば, 平均的な南北温度勾配が大きすぎる点や, 擾乱の特定の特徴があまり現実的でない点である. これらの問題は使用する予報方程式系をより精密なもの (例えばプリミティブ方程式など) にすることで, 解決されると思われる.

数値予報モデルの実装と数値計算ライブラリ「GMathLib」について

今回数値予報モデルを実装する中で, 自作の数値計算ライブラリである「GMath-Lib」を利用した. 本研究と独立した数値計算ライブラリ (行列の取り扱い, 連立方程式のソルバ等を提供している) を利用することで, (主観的ではあるが) 実装に必要なソースコードの量を減らすことができた. それ以上の利点として, プログラムをデバッグする際に, ライブラリを利用しているソースコードの部分をデバッグの対象から外すことができたという点が挙げられる. このために, 比較的短い期間で目的の数値予報モデルを完成させることができた. したがって, 本研究の目標の一つである GMathLib の有効性が確認されたと言えるだろう. 一方で, 実際に数値予報モデルを実装する中で, たくさんの GMathLib の問題点も浮かび上がってきた. その詳細はここでは触れないが, これらの問題点を集計して GMathLib を改良する予定である.

今回作成した数値予報モデルは, 準地衡風予報方程式系を適用した 2 層モデル (矩形領域) である. 本研究では Phillips の数値実験の再現にそれを用いたが, 本質的には南北に境界をもつチャネルの中に生じる流れを数値シミュレーションするならばどのような問題にも, 今回作成したモデルは利用できるだろう (ただし, その数値シミュレーションの結果は準地衡風理論の枠組みの中でしか議論できない欠点がある). たとえば, 作成したモデルを使ってロスビー波の初期値問題を数値シミュレーションするのも興味深いと思われる.

本研究の意義について

最後に, Phillips の数値実験の再現を通して感じた本研究の意義について述べたい. この数値実験で採用された地衡風近似や幾何的な簡単化は, 簡単に取り外すことができる. しかしながら, 数値予報モデルを改良しより精密にすることで, より難しい問題にすぐに突き当たるだろう. それは, 小規模な乱流や対流, 潜熱の解放, 放射の温度依存性, 湿気や雲についての取り扱いである. Phillips の数値実験が行われた時代に比べて現在は, 計算機や数値計算の技術の向上によって, 気象現象の数値シミュレーションの技術もはるかに精巧になった. 現代的な数値予報モデルでは, そのような物理的・数値計算的に取り扱いの難しい問題も考慮するように努力されている. しかし, 大気大循環の基本構造を数値シミュレーションを通して理解したいと思ったときに, そのような精巧な数値予報モデルを用いることは適切でないかもしれない. なぜならば, 本質的でない他のメカニズムに注意を奪われることによって, 大気大循環の基本的原理の理解が妨げられることが考えられるからである. したがって, 現代においても, (準地衡風近似を適用した) 2 層モデルのような粗い数値予報モデルを用いて, 大気大循環の基本的原理を理解した本研究は, 初等的研究の意味で意義があると思われる.

謝辞

本研究を進めるにあたって, 協力頂いた多くの方々に感謝を申し上げます. 指導教官である林祥介教授と岩山隆寛准教授には, 研究発表の場で研究内容についての貴重なコメントをたくさん頂きました. 本研究の出発点は「数値予報モデル」を自作したいという欲求からであり, 岩山隆寛准教授との研究についての楽しい会話の中で, 「Phillips(1956) の再現実験」を紹介して頂きました. その際に, 当研究室に以前所属されていた田原智志さんが, Phillips(1956) の再現実験をされていたことを知りました. 本研究では数値計算の結果の正当性を検証する際に, 田原さんの研究結果を参考にさせてもらいました. 岩山准教授には, 本研究はもちろん他の研究から日頃の研究生活の悩みの相談まで, 自分を激励して下さったことを深く感謝申し上げます. 林祥介教授は授業の場で, 高橋芳幸特任助教や西澤誠也助教は基礎理論読書会などの場で, 自分が教科書内に抱いた疑問点を解説頂き, 地球流体力学の深い理解に毎度導いて下さいました. そこで, 身に付けた知識は本研究に大変生かされていると感じています. 研究室の先輩である島津通さんは, 本論文の体裁をチェックして下さいました. 最後になりましたが, 自分の研究生活をいつも暖かく見守って下さる研究室の先輩方や同輩にも, この場を借りて感謝申し上げます.

なお, 本研究の図の作成には, 電脳 Ruby プロジェクトのソフトウエア (GPhys) を利用させて頂きました.

付録A {鉛直座標変換に伴う微分演 \atop 算子の変換則}

この付録では, 鉛直座標変換において必要となる微分演算子の変換則についてまとめることにする^{*1}.

鉛直座標を表現する独立変数を高度から他の変数に取り換えたいときに, その変数が高度に対して一価関数であれば, z の代わりとなる鉛直座標として用いることができる. そのため, 鉛直方向の独立変数として圧力やエントロピーなどが選ばれる. ここでは, z について一価関数であるような任意の変数 ζ = ζ(x, y, z, t) によって, z 座標系を ζ 座標系へと鉛直座標変換する場合を考えてみることにする.

微分演算子の変換則

水平微分

図 1.0.1: 関数 ψ の空間的な変化を考える. 関数 ψ(x, y, z, t) が,

\begin{aligned}ψ(x, y, z, t) & = ϕ(x, y, z(x, y, ζ, t), t) \\ & = ψ̃(x, y, ζ, t) \end{aligned} ; (1.0.1)

となるように ψ̃(x, y, ζ, t) へと座標変換されたとする. ここで, チルダ (˜ ) は, x, y, ζ, t を独立変数とする関数であることを表している. 図 1.0.1 のように xz 平面上に 3 点 A, B, C をとり, 等 z 面上における ψ の水平微分が, 等 ζ 面上の水平微分とどのような関係を持つかを求める. 点 A, B を含む等 z 面上における点 A での x 方向の水平微分を, 図を参考にして有限差分近似で表し, またそれを変形すると,

*1この付録の作成に当たっては, 地球流体電脳倶楽部理論ノートの「鉛直座標変換」 (林祥介教授) を参考にさせていただいた.

\overline{main.tex 2011/02/28(河合佑太)}

\frac{ψ(B) − ψ(A)}{Δx} = \frac{ψ(C) − ψ(A)}{Δx} − \frac{ψ(C) − ψ(B)}{Δζ} \frac{Δζ}{Δz} \frac{Δz}{Δx}  (1.0.2)

となる. ここで, ψ(A) のように書かれている量は, 点 A における ψ の値を表している.

(1.0.2) において, Δx → 0 とすれば,

(\frac{∂ψ}{∂x} )_{z} = (\frac{∂ψ̃}{∂x})_{ζ} − \frac{∂ψ̃}{∂ζ} \frac{∂ζ}{∂z} (\frac{∂z̃}{∂x} )_{ζ}  (1.0.3)

を得る. y 方向の水平微分についても同様にすれば,

(\frac{∂ψ}{∂y} )_{z} = (\frac{∂ψ̃}{∂y})_{ζ} − \frac{∂ψ̃}{∂ζ} \frac{∂ζ}{∂z} (\frac{∂z̃}{∂y} )_{ζ}  (1.0.4)

を得る.

鉛直微分

(1.0.1) より明らかに,

\frac{∂ψ}{∂z} = \frac{∂ψ̃}{∂ζ} \frac{∂ζ}{∂z}  (1.0.5)

である.

オイラー微分

(1.0.3) を導く際に, xz 平面上における ψ の変化を考えたが, 時間微分についても同様に導ける. すなわち, zt 平面上において局所的な時間変化率を考えればよい. したがって,

(\frac{∂ψ}{∂t} )_{z} = (\frac{∂ψ̃}{∂t})_{ζ} − \frac{∂ψ̃}{∂ζ} \frac{∂ζ}{∂z} (\frac{∂z̃}{∂t})_{ζ}  (1.0.6)

を得る.

\overline{main.tex 2011/02/28(河合佑太)}

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 75}

ラグランジュ微分

(1.0.3) ∼ (1.0.6) の結果を用いて, ラグランジュ微分

\frac{dψ}{dt} = (\frac{∂ψ}{∂t} )_{z} + u (\frac{∂ψ}{∂x} )_{z} + v (\frac{∂ψ}{∂y} )_{z} + w\frac{∂ψ}{∂z}  (1.0.7)

を ζ 座標系へ変換することを考える. (1.0.3) ∼ (1.0.6) を, (1.0.7) に代入すると,

\begin{aligned}\frac{dψ}{dt} = & (\frac{∂ψ̃}{∂t})_{ζ} + ũ (\frac{∂ψ̃}{∂x})_{ζ} + ṽ (\frac{∂ψ̃}{∂y})_{ζ} \\ & − \frac{∂ψ̃}{∂ζ} \frac{∂ζ}{∂z} ((\frac{∂z}{∂t} )_{ζ} + ũ (\frac{∂z}{∂x} )_{ζ} + ṽ (\frac{∂z}{∂y} )_{ζ} − w̃) \end{aligned}  (1.0.8)

となる. ただし,

u(x, y, z, t) = ũ(x, y, ζ, t),  (1.0.9)
v(x, y, z, t) = ṽ(x, y, ζ, t),  (1.0.10)
w(x, y, z, t) = w̃(x, y, ζ, t).  (1.0.11)

である.

(1.0.8) において, ψ = ζ とすると,

\frac{dζ}{dt} = −\frac{∂ζ}{∂z} \{(\frac{∂z}{∂t} )_{ζ} + ũ (\frac{∂z}{∂x} )_{ζ} + ṽ (\frac{∂z}{∂y} )_{ζ} − w̃\}.  (1.0.12)

となるので, (1.0.8) は次のように書き換えられる.

\frac{dψ}{dt} = (\frac{∂ψ̃}{∂t})_{ζ} + ũ (\frac{∂ψ̃}{∂x})_{ζ} + ṽ (\frac{∂ψ̃}{∂y})_{ζ} + ζ̇ \frac{∂ψ̃}{∂ζ}  (1.0.13)

ここで, ζ̇ = dζ/dt である.

したがって, ζ 座標系における全微分は,

\frac{d}{dt} = (\frac{∂}{∂t} )_{ζ} + ũ (\frac{∂}{∂x} )_{ζ} + ṽ (\frac{∂}{∂y} )_{ζ} + ζ̇ \frac{∂}{∂ζ}  (1.0.14)

と書けることが示された.

付録B {圧力座標系におけるプリミ \atop ティブ方程式の導出}

この付録では, z 系で記述されたプリミティブ方程式系を p 系で記述するために鉛直座標変換を行う.

乾燥大気におけるプリミティブ方程式は, 次にあげる運動方程式, 静水圧平衡の式, 連続の式, 理想気体の状態方程式, 熱力学的エネルギー保存の式によって構成される.

\frac{dv}{dt} = −\frac{1}{ρ}∇p − fk × v + F  (2.0.1)
\frac{∂p}{∂z} = −ρg  (2.0.2)
\frac{dρ}{dt} = ρ (∇ · v + \frac{∂w}{∂z})  (2.0.3)
p = ρRT  (2.0.4)
c_{v} \frac{dT}{dt} + p\frac{dα}{dt} = Q̇  (2.0.5)

ここで, v = (u, v) は水平速度ベクトル, w は鉛直速度, ∇ = (∂/∂x, ∂/∂y), ρ は密度, p は圧力, T は温度, α(= ρ^{−1}) は比容, f(= 2Ω sin θ) はコリオリパラメータ, R は乾燥大気の気体定数, c_{v} は乾燥大気の定圧比熱, F は摩擦力, Q̇ は非断熱加熱による加熱率, k は鉛直上向きの単位ベクトルである. また, 全微分と局所微分には,

\frac{d}{dt} = \frac{∂}{∂t} + u\frac{∂}{∂x} + v \frac{∂}{∂y} + w\frac{∂}{∂z}  (2.0.6)

の関係がある.

また, 熱力学的エネルギー保存の式 (2.0.5) は, 温位

θ = T (\frac{p_{s}}{p} )^{κ}  (2.0.7)

を使って,

\frac{dln θ}{dt} = \frac{1}{c_{p}T} Q̇  (2.0.8)

と書き直せる. ここで, p_{s} は基準となる定数の圧力を任意に選べばよいが, 地上気圧 (1000 hPa) を取ることが多い. また, κ = R/c_{p} である.

次に, z 系のプリミティブ方程式 (2.0.1) ∼ (2.0.5) において, 鉛直座標を幾何的な高度から圧力に取り直すことを考える. z 座標系から任意の鉛直座標系への変換に伴う微分演算子の変換則ついては, 付録 A にまとめた. ここでは, 付録 A における任意の鉛直方向の独立変数 ζ が圧力 p であるとして, その結果を用いることにする. 付録 A の ζ̇ はこの場合 dp/dt となり, これを ω と定義する. この ω は, p 座標系での鉛直速度を表すことになる. p 座標系におけるラグランジュ微分は, 付録 A の (1.0.14) から次のように書ける.

\frac{d}{dt} = (\frac{∂}{∂t} )_{p} + ũ (\frac{∂}{∂x} )_{p} + ṽ (\frac{∂}{∂y} )_{p} + ω̃ \frac{∂}{∂p}  (2.0.9)

ここで, 添え字の p は等圧力面上での局所微分を表わしている. また, チルダのついている変数は, 独立変数 (x, y, p, t) の従属変数であることを表すとする.

運動方程式 (2.0.1) の鉛直座標変換加速度項のラグランジュ微分は, (2.0.9) によって変換される. 圧力傾度力の項は, 付録 A の (1.0.3), 静水圧平衡の式 (2.0.2) を用いれば次のように変換される.

\begin{aligned}−\frac{1}{ρ} (\frac{∂p}{∂x} )_{z} & = −\frac{1}{ρ̃} [(\frac{∂p̃}{∂x})_{p} − \frac{∂p}{∂p} \frac{∂p}{∂z} (\frac{∂z̃}{∂x})_{p}] \\ & = −g (\frac{∂z̃}{∂x})_{p} = − (\frac{∂Φ̃}{∂x})_{p} \end{aligned}  (2.0.10)

ここで, Φ̃(= gz) はジオポテンシャルである.

したがって, p 系における運動方程式は,

\frac{dṽ}{dt} = −∇_{p}Φ̃ − fk × ṽ + F̃  (2.0.11)

となる. ここで,

∇_{p} = ((\frac{∂}{∂x} )_{p} , (\frac{∂}{∂y} )_{p}) .

静水圧平衡の式 (2.0.2) の鉛直座標変換付録 A の (1.0.5) において, ψ = z, ψ̃ = z̃ とすると,

\begin{aligned}\frac{∂z}{∂z} & = \frac{∂z̃}{∂p} \frac{∂p}{∂z} \\ & = −ρ̃g \frac{∂Φ̃}{∂p} = −ρ̃\frac{∂Φ̃}{∂p} = 1 \end{aligned}

となる. よって, 理想気体の状態方程式 (2.0.4) を用いれば, p 系における静水圧平衡の式

\frac{∂Φ̃}{∂p} = −\frac{RT̃}{p}  (2.0.12)

を得る.

連続の式 (2.0.3) の鉛直座標変換連続の式 (2.0.3)

−\frac{1}{ρ} (\frac{dρ}{dt})_{z} = (\frac{∂u}{∂x})_{z} + (\frac{∂v}{∂y} )_{z} + \frac{∂w}{∂z}

の z 系での各微分に対して, (2.0.9), 付録 A の (1.0.3), (1.0.4), (1.0.5) を適用する. また, 付録 A の (1.0.12) より,

ω̃ = −\frac{∂p}{∂z}\{(\frac{∂z̃}{∂t})_{p} + ũ (\frac{∂z̃}{∂x})_{p} + ṽ (\frac{∂z̃}{∂y})_{p} − w̃\}

であるので, 静水圧平衡の式を用いれば w と ω の関係式

w̃ = −\frac{ω̃}{ρg} + (\frac{∂z̃}{∂t})_{p} + ũ (\frac{∂z̃}{∂x})_{p} + ṽ (\frac{∂z̃}{∂y})_{p}

を得る. したがって, 連続の式は,

\begin{aligned}−\frac{1}{ρ̃} & [(\frac{∂ρ̃}{∂t})_{p} + ũ (\frac{∂ρ̃}{∂x})_{p} + ṽ (\frac{∂ρ̃}{∂y})_{p} + ω̃ \frac{∂ρ̃}{∂p}] \\ = & (\frac{∂ũ}{∂x} )_{p} + (\frac{∂ṽ}{∂y} )_{p} + ρ̃g [(\frac{∂z̃}{∂x})_{p} \frac{∂ũ}{∂p} + (\frac{∂z̃}{∂y})_{p} \frac{∂ṽ}{∂p}] \\ & + \frac{∂ω̃}{∂p} − \frac{ω̃}{ρ̃} \frac{∂ρ̃}{∂p} − ρ̃g \frac{∂}{∂p} [(\frac{∂z̃}{∂t})_{p} + ũ (\frac{∂z̃}{∂x})_{p} + ṽ (\frac{∂z̃}{∂y})_{p}] \end{aligned}  (2.0.13)

となる. さらに, 上式の右辺第 6 項をさらに展開すれば,

−\frac{1}{ρ̃} [(\frac{∂ρ̃}{∂t})_{p} + ũ (\frac{∂ρ̃}{∂x})_{p} + ṽ (\frac{∂ρ̃}{∂y})_{p}] − ρ̃g [\frac{∂ũ}{∂p} (\frac{∂z̃}{∂x})_{p} + \frac{∂ṽ}{∂p} (\frac{∂z̃}{∂y})_{p}]

となるので, 多くの項は打ち消しあうことが分かる.

最終的に, p 系における連続の式は次のような簡潔な形で与えられる.

(\frac{∂ũ}{∂x})_{p} + (\frac{∂ṽ}{∂y})_{p} + \frac{∂ω̃}{∂p} = 0.  (2.0.14)

あたかも z 系での非圧縮流体に対する連続の式のようであるが, ここでは, 非圧縮流体の仮定を置いていないことに注意されたい.

熱力学的エネルギー保存の式 (2.0.5) の鉛直座標変換ラグランジュ微分は (2.0.9) によって書かれるが, p 系における熱力学的エネルギー保存の式は, 見かけ上 z 系のものと変わらない.

c_{v} \frac{dT̃}{dt} + p\frac{dα̃}{dt} = ˙ Q̃  (2.0.15)

付録C {準地衡風予報方程式系の導 \atop 出における各式の変形}

この付録では, 2.1 節の準地衡風予報方程式系の部分で省略した各式の変形の詳細について説明する.

(2.1.12) の導出

熱力学的エネルギー保存の式 (2.1.5)

\frac{dln θ}{dt} = \frac{1}{c_{p}T} Q̇

において, 非断熱効果を考慮しなければ右辺はゼロとなる. 上式の左辺を変形する際に, 静水圧平衡の式 (2.1.10), 理想気体の状態方程式

\frac{∂Φ}{∂p} = −α = −\frac{RT}{p}

を用いれば,

\frac{dln θ}{dt} = \frac{d}{dt} (ln T − κln p)

= \frac{d}{dt} [ln (\frac{pα}{R}) − κln p] = \frac{d}{dt} [ln α + (1 − κ)ln p]

= −\frac{1}{α} \frac{d}{dt} (\frac{∂Φ}{∂p}) + (1 − κ)\frac{1}{p} \frac{dp}{dt}

= −\frac{1}{α} [\frac{∂}{∂t} (\frac{∂Φ}{∂p}) + v · ∇\frac{∂Φ}{∂p} + ω \frac{∂}{∂p} (\frac{∂Φ}{∂p})] + \frac{(1 − κ)ω}{p} = 0.

(3.0.1)

と変形される. さらに, 両辺に −α = ∂Φ/∂p を乗ずれば,

\frac{∂}{∂t} (\frac{∂Φ}{∂p}) + v · ∇\frac{∂Φ}{∂p} + ω \frac{∂}{∂p} (\frac{∂Φ}{∂p}) + \frac{(1 − κ)ω}{p} \frac{∂Φ}{∂p} = 0.  (3.0.2)

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 81} となるので, 結果 (2.1.12)

\frac{∂}{∂t} (\frac{∂Φ}{∂p}) + v · ∇\frac{∂Φ}{∂p} + \frac{ω}{p} [\frac{∂}{∂p} (p\frac{∂Φ}{∂p} − κΦ)] = 0.  (3.0.3)

を得る.

(2.1.18) の最右辺の導出

(2.1.18) の最右辺を得るために, 静水圧平衡の式 (2.1.10) を用いて中辺を次のように変形する.

R_{i} = p\frac{d}{dp} [p\frac{dΦ̄}{dp} − κΦ̄] \frac{1}{U^{2}}

= \frac{p^{2}}{U^{2}} (\frac{d^{2}Φ̄}{dp^{2}} + \frac{1 − κ}{p} \frac{dΦ̄}{dp}) = −\frac{p^{2}}{U^{2}} (\frac{∂ᾱ}{∂p} + \frac{1 − κ}{p}ᾱ)

= −\frac{p^{2}ᾱ}{U^{2}} (\frac{∂ln ᾱ}{∂p} + (1 − κ)\frac{∂ln p}{∂p})

ここで, 温位の p による対数微分は理想気体の状態方程式を用いて,

\frac{∂ln θ}{∂p} = \frac{∂}{∂p} ln (T (p_{s}/p)^{κ}) = \frac{∂}{∂p} (ln ᾱ + (1 − κ)ln p)

となることに注意すれば,

R_{i} = −\frac{p^{2}ᾱ}{U^{2}} \frac{∂ln θ̄}{∂p}  (3.0.4)

となり, (2.1.18) の最右辺を得ることができる.

(2.1.23) の導出

今, ジオポテンシャルを (2.1.14)

Φ(x, y, p, t) = Φ̄(p) + f_{0}ULΦ^{∗}(x, y, p, t)

のように分割している. (2.1.12) に上の Φ を代入すれば,

\begin{aligned}f_{0}UL [\frac{∂}{∂t} (\frac{∂Φ^{∗}}{∂p} ) + v · ∇\frac{∂Φ^{∗}}{∂p} ] & + \frac{ω}{p^{2}} [p\frac{d}{dp} (p\frac{dΦ̄}{dp} − κΦ̄)] \\ & + \frac{ω}{p} f_{0}UL\frac{∂}{∂p} (p\frac{∂Φ^{∗}}{∂p} − κΦ^{∗}) = 0. \end{aligned} (3.0.5)

\overline{main.tex 2011/02/28(河合佑太)}

上式の各変数を無次元化し, さらにリチャードソン数 (2.1.18) を導入すれば,

\begin{aligned}\frac{f_{0}U^{2}}{P} & (\frac{∂}{∂t^{∗}} + u^{∗} \frac{∂}{∂x^{∗}} + v^{∗} \frac{∂}{∂y^{∗}} ) (\frac{∂Φ^{∗}}{∂p^{∗}} ) \\ & + \frac{U}{P L} \frac{ω^{∗}}{(p^{∗})^{2}} [U^{2}R_{i} − f_{0}ULp^{∗} \frac{∂}{∂p^{∗}} (p^{∗} \frac{∂Φ^{∗}}{∂p^{∗}} + κΦ^{∗})] = 0 \end{aligned}

となる. 最後に両辺に P/f_{0}^{2}UL をかけて, U/f_{0}L をロスビー数 R_{o} で置き換えれば,

\begin{aligned}R_{o} & (\frac{∂}{∂t^{∗}} + u^{∗} \frac{∂}{∂x^{∗}} + v^{∗} \frac{∂}{∂y^{∗}} ) (\frac{∂Φ^{∗}}{∂p^{∗}} ) \\ & + \frac{ω^{∗}}{(p^{∗})^{2}} [R_{i}R_{o}^{2} − R_{o}p^{∗} \frac{∂}{∂p^{∗}} (p^{∗} \frac{∂Φ^{∗}}{∂p^{∗}} + κΦ^{∗})] = 0 \end{aligned}  (3.0.6)

を得る.

(2.1.34) の導出

R_{o} の 1 次についての方程式系 (2.1.29)∼(2.1.32) において, u_{1}^{∗}, v_{1}^{∗}, Φ_{1}^{∗} を消去するために, (2.1.30) に ∂/∂x, (2.1.29) に −∂/∂y を作用させて, 両辺を加え合わせれば,

\begin{matrix}[(\frac{∂}{∂t^{∗}} + u_{0}^{∗} \frac{∂}{∂x^{∗}} + v_{0}^{∗} \frac{∂}{∂y^{∗}} ) + (\frac{∂u_{0}^{∗}}{∂x^{∗}} + \frac{∂v_{0}^{∗}}{∂y^{∗}} )] & (\frac{∂v_{0}^{∗}}{∂x^{∗}} − \frac{∂u_{0}^{∗}}{∂y^{∗}} ) \\ + β^{∗}v_{0}^{∗} + β^{∗}y^{∗} (\frac{∂u_{0}^{∗}}{∂x^{∗}} + \frac{∂v_{0}^{∗}}{∂y^{∗}} ) + (\frac{∂u_{1}^{∗}}{∂x^{∗}} + \frac{∂v_{1}^{∗}}{∂y^{∗}} ) = 0 \end{matrix}  (3.0.7)

を得る.

ところで, (2.1.26) および (2.1.25) より, 添え字 0 の風速場は非発散である. すなわち,

\frac{∂u_{0}^{∗}}{∂x^{∗}} + \frac{∂v_{0}^{∗}}{∂y^{∗}} = 0.

さらに, 添え字 1 の風速場の発散に (2.1.31) を代入し, 地衡風の渦度 (2.1.33) を導入すれば, 結局 (3.0.7) は,

\frac{∂ζ_{0}^{∗}}{∂t^{∗}} + u_{0}^{∗} \frac{∂ζ_{0}^{∗}}{∂x^{∗}} + v_{0}^{∗} (β^{∗} + \frac{∂ζ_{0}^{∗}}{∂y} ) − \frac{∂ω_{1}^{∗}}{∂p^{∗}} = 0  (3.0.8)

となり, (2.1.34) が得られる.

付録D Phillips の数値予報モデル に関する数式の導出

２層モデルにおけるポテンシャル渦度方程式 (2.2.80), (2.2.81) の導出

レベル 1 と 3 で適用される渦度方程式 (2.2.76) に含まれる鉛直速度 ω_{2} を, 熱力学的エネルギー保存の式 (2.2.78) 使って消去する. その際に, 相対的なポテンシャル渦度 (2.2.79) の定義

\begin{aligned}q_{1} & = ζ_{1} − λ^{2}(ψ_{1} − ψ_{3}) = ∇^{2}ψ_{1} − λ^{2}(ψ_{1} − ψ_{3}) \\ q_{3} & = ζ_{3} + λ^{2}(ψ_{1} − ψ_{3}) = ∇^{2}ψ_{3} + λ^{2}(ψ_{1} − ψ_{3}) \end{aligned}

を利用すれば,

\frac{∂q_{1}}{∂t} = − v_{1} · ∇(βy + ζ_{1}) + v_{2} · ∇[λ^{2}(ψ_{1} − ψ_{3})]
+ A ∇^{2}q_{1} + \frac{2RHλ^{2}}{f_{0}c_{p}} \frac{y}{W}  (4.0.1)
\frac{∂q_{3}}{∂t} = − v_{3} · ∇(βy + ζ_{3}) − v_{2} · ∇[λ^{2}(ψ_{1} − ψ_{3})]
+ A ∇^{2}q_{3} − \frac{2RHλ^{2}}{f_{0}c_{p}} \frac{y}{W} − kζ_{4}  (4.0.2)

を得る.

(4.0.1) の v_{2} と (4.0.2) の ζ_{4} は, ２層モデルにおいて予報変数ではないので直接的には分からない. よって, それらの変数は, 予報されるレベル 1 および 3 の風速場を用いて線形補間することによって与える. すなわち,

v_{2} = \frac{1}{2} (v_{1} + v_{3})  (4.0.3)
ζ_{4} = \frac{1}{2} (3ζ_{3} − ζ_{1}) = \frac{1}{2} [3q_{3} − q_{1} − 4λ^{2}(ψ_{1} − ψ_{3})]  (4.0.4)

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 84} とすることで計算される.

よって (4.0.1) の右辺第 1, 2 項は,

\begin{aligned}−v_{1} · ∇(βy + ζ_{1}) + v_{2} · ∇(ψ_{1} − ψ_{3}) = & − v_{1} · ∇(βy + q_{1}) \\ & − λ^{2} \frac{v_{1} − v_{3}}{2} · ∇(ψ_{1} − ψ_{3}) \end{aligned}

と変形される. さらに, 地衡風の流線関数を用いて上式の右辺第 2 項の変形を続ければ,

(v_{1} − v_{3}) · ∇(ψ_{1} − ψ_{3}) = J(ψ_{1} − ψ_{3}, ψ_{1} − ψ_{3}) = 0

となる. ここで, J はヤコビオペレータ

J(α, β) = \frac{∂α}{∂x} \frac{∂β}{∂y} − \frac{∂α}{∂y} \frac{∂β}{∂x}

である. (4.0.2) の右辺第 1, 2 項についても同様に変形すれば,

−v_{3} · ∇(βy + ζ_{3}) + v_{2} · ∇(ψ_{1} − ψ_{3}) = −v_{3} · ∇(βy + q_{3})

を得る.

したがって, (4.0.1) と (4.0.2) は最終的に,

\frac{∂q_{1}}{∂t} = −v_{1} · ∇(βy + q_{1}) + A∇^{2}q_{1} + \frac{2RHλ^{2}}{f_{0}c_{p}} \frac{y}{W}  (4.0.5)

\frac{∂q_{3}}{∂t} = −v_{3} · ∇(βy + q_{3}) + A∇^{2}q_{3} − \frac{2RHλ^{2}}{f_{0}c_{p}} \frac{y}{W}

−\frac{k}{2}[3q_{3} − q_{1} − 4λ^{2}(ψ_{1} − ψ_{3})]  (4.0.6)

となり, 今回の 2 層モデルにおけるポテンシャル渦度方程式を得る.

付録E Phillips の数値実験におけ るエネルギー方程式の導出 の詳細

第 4 章において系のエネルギーサイクルを説明したが, そこで示した各種方程式の導出についてこの付録にて説明する.

エネルギー方程式の導出

平均帯状流の運動エネルギーの時間発展の式 (2.3.95) の導出

(2.3.89) と (2.3.90) にそれぞれ −ψ̄_{1} と −ψ̄_{3} を掛けて加え合わし, その結果を領域全体で積分することを考えるが, 2 式は非常に似た形式を持つので, (2.3.89) についてひとまず考える. (2.3.89) に −ψ_{1} を掛けて, 領域全体で積分すれば

∫ (−ψ̄_{1} ∂_{yyt}ψ̄_{1} − ψ̄_{1}ω̄_{2}f_{0}/p_{2})dσ = ∫ (−Aψ̄_{1} ∂_{yy}ζ̄_{1} + ψ̄_{1}\overline{v1^{′} · ∇ζ1^{′}})dσ.  (5.0.1)

左辺第 1 項目の積分は, 南北壁における境界条件 (ū_{1} = 0, ∂^{2}ū_{1}/∂y^{2} = 0) に注意すれば

\begin{aligned}− ∫ ψ̄_{1} ∂_{yyt}ψ̄_{1} dσ & = ∫ ψ̄_{1} ∂_{yt}ū_{1} dσ = ∫ [∂_{y}(ψ̄_{1} ∂_{t}u_{1}) − ∂_{y}ψ̄_{1}∂_{t}ū_{1}] dσ \\ & = ∫ [ψ̄_{1} ∂_{t}ū_{1}]_{−W}^{W} dx + ∫ ∂_{t} \frac{ū_{1}^{2}}{2} dσ = ∫ ∂_{t} \frac{ū_{1}^{2}}{2} dσ. \end{aligned}

となる.

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 86} また, 右辺第 1 項の積分は

−A ∫ ψ̄_{1}∂_{yy}ζ̄_{1}dσ = −A ∫ [∂_{y}(ψ̄_{1}∂_{y}ζ̄_{1}) − ∂_{y}ψ̄_{1}∂_{y}ζ̄_{1}]dσ
= A ∫ [ψ̄_{1}∂_{y}ζ̄_{1}]_{−W}^{W} dx − A ∫ [∂_{y}(ζ̄_{1} ∂_{y}ψ̄_{1}) − ζ̄_{1} ∂_{yy}ψ̄_{1}] dσ
= −A ∫ [ζ̄_{1} ∂_{y}ψ̄_{1}]_{−W}^{W} dx + A ∫ ζ̄_{1} ∂_{yy}ψ̄_{1}dσ = A ∫ ζ̄_{1}^{2}dσ

となる.

次に右辺第 3 項の積分を変形する. まず, 上線の付いている部分は

\begin{matrix}\overline{v1^{′} · ∇ζ1^{′}} & = −\overline{∂_{y}ψ1^{′} ∂_{x}ψ1^{′}} + \overline{∂_{x}ψ1^{′} ∂_{y}ψ1^{′}} \\ = & − \overline{∂_{x}(ζ1^{′} ∂_{y}ψ1^{′})} + \overline{ζ1^{′}∂_{xy}ψ1^{′}} + \overline{∂_{y}(ζ1^{′} ∂_{x}ψ1^{′})} − \overline{ζ1^{′}∂_{xy}ψ1^{′}} \\ = & \overline{∂_{y}(ζ1^{′} ∂_{x}ψ1^{′})} \end{matrix}

となる. また, 渦による南北方向の渦度フラックスが, 水平方向の渦による運動量フラックスの収束に等しいという関係式

\begin{aligned}\overline{ζ1^{′} ∂_{x}ψ1^{′}} & = \overline{∇^{2}ψ1^{′} ∂_{x}ψ1^{′}} \\ & = \overline{∂_{xx}ψ1^{′} ∂_{x}ψ1^{′}} + \overline{∂_{yy}ψ1^{′} ∂_{x}ψ1^{′}} \\ & = \overline{∂_{x} \frac{(ψ1^{′})^{2}}{2}} + \overline{∂_{y}(∂_{x}ψ1^{′}∂_{y}ψ1^{′})} − \overline{∂_{x} \frac{(∂_{y}ψ1^{′})^{2}}{2}} \\ & = −∂_{y}(\overline{u1^{′}v1^{′}}) \end{aligned}  (5.0.2)

を用いれば, 結局のところ右辺第 3 項の積分の変形は

\begin{aligned}− ∫ (−ψ̄_{1})\overline{v1^{′} · ∇ζ1^{′}}dσ & = − ∫ ψ̄_{1}∂_{y}(\overline{u1^{′}v1^{′}})dσ \\ & = ∫ [ψ̄_{1}\overline{u1^{′}v1^{′}}]_{−W}^{W} dx + ∫ ∂_{y}ψ̄_{1}\overline{u1^{′}v1^{′}}dσ \\ & = − ∫ ū_{1}\overline{u1^{′}v1^{′}}dσ \end{aligned}  (5.0.3)

となる.

一方, (2.3.90) の右辺第 3 項と −ψ̄_{3} の積についての積分は,

\begin{aligned}∫ (−ψ̄_{3}) · (−kζ̄_{4})dσ & = −k ∫ ψ̄_{3} ∂_{y}ū_{4}dσ \\ & = −k ∫ [ψ̄_{3}ū_{4}]_{−W}^{W} dx − k ∫ ∂_{y}ψ̄_{3}ū_{4}dσ \\ & = −k ∫ ū_{3}ū_{4}dσ \end{aligned}  (5.0.4)

と計算される. (2.3.90) のそれ以外の項についての積分は, (2.3.89) についての積分と同様に行えばよい.

以上より, 平均帯状流の運動エネルギーの時間発展の式 (2.3.95) が導出される.

擾乱運動の運動エネルギーの時間発展の式 (2.3.96) の導出

(2.3.92) と (2.3.93) にそれぞれ −ψ1^{′} と −ψ3^{′} を掛けて加え合わせる. そして, その結果を境界条件に注意して領域全体で積分する. 平均帯状流の運動エネルギーを求めたときのように, ひとまず (2.3.92)× (−ψ1^{′}) についての積分について考えることにする. このとき, その式の左辺第 1 項目の積分は

\begin{aligned}− ∫ ψ1^{′}∇^{2}∂_{t}ψ1^{′}dσ & = − ∫ [ψ1^{′} ∂_{tx}ψ1^{′}]_{0}^{L} dy + ∫ ∂_{t}∂_{x} \frac{(ψ1^{′})^{2}}{2} dσ \\ & − ∫ [ψ1^{′} ∂_{t}(∂_{y}ψ1^{′})]_{−W}^{W} dx + ∫ ∂_{t}∂_{y} \frac{(ψ1^{′})^{2}}{2} dσ \\ & = ∂_{t} ∫ \frac{(∇ψ1^{′})^{2}}{2} dσ \end{aligned}  (5.0.5)

となる.

また右辺第 1 項目の積分は,

\begin{matrix}−A ∫ ψ1^{′}∇^{2}ζ1^{′}dσ = & − A ∫ [ψ1^{′}∂_{x}ζ1^{′}]_{0}^{L}dy + A ∫ ∂_{x}ψ1^{′}∂_{x}ζ1^{′}dσ \\ − A ∫ [ψ1^{′}∂_{y}ζ1^{′}]_{−W}^{W} dx + A ∫ ∂_{y}ψ1^{′}∂_{y}ζ1^{′}dσ \\ = & A ∫ (∂_{x}ψ1^{′}∂_{x}ζ1^{′} + ∂_{y}ψ1^{′}∂_{y}ζ1^{′})dσ \\ = & A ∫ [∂_{x}ψ1^{′}∇^{2}ψ1^{′}]_{0}^{L}dy + A ∫ [∂_{y}ψ1^{′}∇^{2}ψ1^{′}]_{−W}^{W} dx − A ∫ (∇^{2}ψ1^{′})^{2}dσ \\ = & − A ∫ (ζ1^{′})^{2}dσ \end{matrix} ; (5.0.6)

と計算される.

続いて右辺第 2 項目の積分計算を行う. 変数を平均場と擾乱場に分解して, 被積 \underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 88} 分関数を変形すれば

\begin{matrix}ψ1^{′}[v_{1} · ∇(βy & + ζ_{1})]′ = ψ1^{′}v_{1} · ∇(βy + ζ_{1}) − ψ1^{′}\overline{v_{1} · ∇(βy + ζ_{1})} \\ = & ψ1^{′}[v̄_{1} · ∇(\overline{βy + ζ_{1}}) + v̄_{1} · ∇ζ1^{′} + v1^{′} · ∇ζ̄_{1} + βv1^{′} + v1^{′} · ∇ζ1^{′}] \\ − ψ1^{′}\overline{v1^{′} · ∇ζ1^{′}} \\ = & ψ1^{′}v̄_{1} · ∇ζ1^{′} + ψ1^{′}v1^{′} · ∇ζ̄_{1} + βψ1^{′}v1^{′} + ψ1^{′}v1^{′} · ∇ζ1^{′} − ψ1^{′}\overline{v1^{′} · ∇ζ1^{′}} \end{matrix} (5.0.7)

となる. ここで, V̄_{1} = (−∂_{y}ψ̄_{1}, 0) と ∇\overline{(βy + ζ_{1})} = (0, β + ∂_{y}ζ_{1}) より

\overline{V̄_{1} · ∇(βy + ζ̄_{1})} = V̄_{1} · \overline{(βy + ζ_{1})} = 0

であることを用いた. さらに, (5.0.7) の各項を変形していく. まず, (5.0.7) の最右辺第 1, 2 項の積分は次のように式変形できる.

\begin{matrix}∫ (ψ1^{′} & v̄_{1} · ∇ζ1^{′} + ψ1^{′}v1^{′} · ∇ζ̄_{1})dσ = ∫ (ψ1^{′}ū_{1}ζ1^{′} + ψ1^{′}v1^{′}∂_{y} ζ̄_{1}) dσ \\ = & ∫ [ψ1^{′}ū_{1}ζ1^{′}]_{0}^{L}dy − ∫ ζ1^{′}ū_{1} ∂_{x}ψ1^{′}dσ + ∫ \frac{1}{2}∂_{x}(ψ1^{′})^{2}∂_{y}ζ̄_{1} dσ \\ = & − ∫ ζ1^{′}ū_{1} ∂_{x}ψ1^{′} dσ + \frac{1}{2} ∫ [(ψ1^{′})^{2}]_{0}^{L}∂_{y}ζ̄_{1} dy \\ = & − ∫_{−W}^{W} (∫_{0}^{L} ζ1^{′}v1^{′} dx) ū_{1} dy \\ = & − L ∫_{−W}^{W} ū_{1}\overline{ζ1^{′}v1^{′}} dy = ∫ ū_{1} ∂_{y}(\overline{u1^{′}v1^{′}}) dσ. \end{matrix}

なお, 最後の式変形に (5.0.2) を用いた. (5.0.7) の最右辺第 3 項の積分は,

∫ ψ1^{′}v_{1}′βdσ = β ∫ [\frac{(ψ1^{′})^{2}}{2} ]_{0}^{L} dy = 0

となる. (5.0.7) の最右辺第 4 項は,

∫ ψ1^{′}v1^{′} · ∇ζ_{1}dσ = ∫ ψ1^{′}(−∂_{y}ψ1^{′}∂_{x}ζ1^{′} + ∂_{x}ψ1^{′}∂_{y}ζ1^{′}) dσ
= − \frac{1}{2} ∫ ∂_{y}(ψ1^{′})^{2} ∂_{x}ζ1^{′} dσ + \frac{1}{2} ∫ ∂_{x}(ψ1^{′})^{2} ∂_{y}ζ1^{′} dσ
= − \frac{1}{2} ∫ [ζ1^{′} ∂_{y}(ψ1^{′})^{2}]_{0}^{L} dσ + \frac{1}{2} ∫ ζ1^{′} ∂_{xy}(ψ1^{′})^{2} dσ

+ \frac{1}{2} ∫ [ζ1^{′} ∂_{y}(ψ1^{′})^{2}]_{−W}^{W} dx − \frac{1}{2} ∫ ζ1^{′} ∂_{xy}(ψ1^{′})^{2} dσ

= 0.

(5.0.7) の最右辺第 5 項は,

∫ ψ1^{′} \overline{v1^{′} · ∇ζ1^{′}} dσ = L ∫ \overline{v1^{′} · ∇ζ1^{′}}ψ̄1^{′} dy = 0

\overline{main.tex 2011/02/28(河合佑太)}

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 89} となる.

以上より, (2.3.92)× (−ψ1^{′}) を領域全体で積分した結果は,

\frac{∂}{∂t} ∫ \frac{(∇ψ1^{′})^{2}}{2} dσ = −\frac{f_{0}}{c_{p}} ∫ ω2^{′}ψ1^{′} dσ − A ∫ (ζ1^{′})^{2} dσ + ∫ ū_{1} \frac{∂(\overline{u1^{′}v1^{′}})}{∂y} dσ

となる.

一方, (2.3.93) × (−ψ3^{′}) を領域全体で積分するとき, その右辺第 3 項の積分以外は上の結果の添え字を変えるだけでよい. なお, 右辺第 3 項の積分は

\begin{matrix}∫ ψ3^{′}kζ4^{′} dσ = & k ∫ ψ3^{′} ∂_{x}v4^{′} dσ − k ∫ ψ3^{′}∂_{y}u4^{′} dσ \\ = & k ∫ [v4^{′} ψ3^{′}]_{0}^{L} dy − k ∫ v4^{′} ∂_{x}ψ3^{′} dσ \\ − & k ∫ [u4^{′}ψ3^{′}]_{−W}^{W} dx + k ∫ u4^{′} ∂_{y}ψ3^{′}dσ \\ = & − k ∫ (u3^{′}u4^{′} + v3^{′}v4^{′})dσ = −k ∫ v3^{′} · v4^{′}dσ \end{matrix}

となる.

以上より, 擾乱運動の運動エネルギーの時間発展の式 (2.3.96) が導出される.

平均帯状流のポテンシャルエネルギーの時間発展の式(2.3.97)の導出

(2.3.91) に (ψ̄_{1} − ψ̄_{3}) を掛けると,

\begin{matrix}(ψ̄_{1} − ψ̄_{3})\frac{ω̄_{2}f_{0}}{p_{2}} = & λ^{2} [(ψ̄_{1} − ψ̄_{3})\frac{∂}{∂t} (ψ̄_{1} − ψ̄_{3}) + (ψ̄_{1} − ψ̄_{3})\overline{v2^{′} · ∇(ψ1^{′} − ψ3^{′})}] \\ + & (ψ̄_{1} − ψ̄_{3})\frac{2RH}{f_{0}c_{p}} \frac{y}{W} − A(ψ̄_{1} − ψ̄_{3})∂_{yy}(ψ̄_{1} − ψ̄_{3}) \end{matrix} (5.0.8)

\overline{main.tex 2011/02/28(河合佑太)}

となり, この式の領域全体での積分を計算することになる. 右辺第 2 項についての積分を計算すれば,

\begin{aligned}λ^{2} & ∫ (ψ̄_{1} − ψ̄_{3}) \overline{v2^{′} · ∇(ψ1^{′} − ψ3^{′})} dσ \\ = & λ^{2} ∫ (ψ̄_{1} − ψ̄_{3}) \overline{u2^{′} ∂_{x}(ψ1^{′} − ψ3^{′})} dσ + λ^{2} ∫ (ψ̄_{1} − ψ̄_{3}) \overline{v2^{′} ∂_{y}(ψ1^{′} − ψ3^{′})} dσ \\ = & − λ^{2} ∫ (ψ̄_{1} − ψ̄_{3}) \overline{(ψ1^{′} − ψ3^{′}) ∂_{x}u2^{′}} dσ + λ^{2} ∫ [(ψ̄_{1} − ψ̄_{3}) \overline{v2^{′}(ψ1^{′} − ψ3^{′})}]_{−W}^{W} dx \\ & − λ^{2} ∫ \overline{(ψ1^{′} − ψ3^{′}) ∂_{y}\{v2^{′}(ψ̄_{1} − ψ̄_{3})}\} dσ \\ = & − λ^{2} ∫ (ψ̄_{1} − ψ̄_{3})\overline{(ψ1^{′} − ψ3^{′})(∂_{x}u2^{′} + ∂_{y}v2^{′})} dσ − λ^{2} ∫ \overline{v2^{′}(ψ1^{′} − ψ3^{′})} ∂_{y}(ψ̄_{1} − ψ̄_{3}) dσ \end{aligned} ; (5.0.9)

となる. ところで, レベル 2 において連続の式を評価すれば

\frac{∂u_{2}}{∂x} + \frac{∂v_{2}}{∂y} = −\frac{∂ω}{∂p} = −\frac{ω_{0} − ω_{4}}{2p_{2}} = 0

となり, 速度変数を平均場と擾乱場に分けることで

\frac{∂u2^{′}}{∂x} + \frac{∂v2^{′}}{∂y} = 0

を得る. これにより (5.0.9) は結局

λ^{2} ∫ (ψ̄_{1} − ψ̄_{3}) \overline{v2^{′} · ∇(ψ1^{′} − ψ3^{′})} dσ = −λ^{2} ∫ \overline{v2^{′}(ψ1^{′} − ψ3^{′})} ∂_{y}(ψ̄_{1} − ψ̄_{3}) dσ (5.0.10)

となる.

次に右辺第 4 項の積分を計算すれば,

\begin{matrix}−λ^{2} & A ∫ (ψ̄_{1} − ψ̄_{3}) ∂_{yy}(ψ̄_{1} − ψ̄_{3}) dσ \\ = & − λ^{2}A ∫ [(ψ̄_{1} − ψ̄_{3}) · (ū_{3} − ū_{1})]_{−W}^{W} dx + λ^{2}A ∫ [∂_{y}(ψ̄_{1} − ψ̄_{3})]^{2} dσ \\ = & λ^{2}A ∫ [∂_{y}(ψ̄_{1} − ψ̄_{3})]^{2} dσ \end{matrix} (5.0.11)

となる.

以上より, 平均帯状流のポテンシャルエネルギーの時間発展の式 (2.3.97) が導出される.

\overline{main.tex 2011/02/28(河合佑太)}

\underline{準地衡風 2 層モデルを用いた中緯度大気循環の研究 91}

擾乱運動のポテンシャルエネルギーの時間発展の式 (2.3.98) の導出

(2.3.94) に (ψ1^{′} − ψ3^{′}) を掛けると,

\begin{matrix}(ψ1^{′} − ψ3^{′})\frac{ω2^{′}f_{0}}{p_{2}} = & λ^{2}\{(ψ1^{′} − ψ3^{′})\frac{∂}{∂t} (ψ1^{′} − ψ3^{′}) + (ψ1^{′} − ψ3^{′}) [v2^{′} · ∇(ψ1^{′} − ψ3^{′})]′ \\ − & A(ψ1^{′} − ψ3^{′})∇^{2}(ψ1^{′} − ψ3^{′})\} \end{matrix} ; (5.0.12)

となり, この式の領域全体での積分を計算することになる.

右辺第 2 項の積分を考える. 最初に被積分関数を変形すれば

\begin{matrix}(ψ1^{′} & − ψ3^{′}) [v_{2} · ∇(ψ_{1} − ψ_{3})]′ \\ = & (ψ1^{′} − ψ3^{′})v_{2} · ∇(ψ_{1} − ψ_{3}) − (ψ1^{′} − ψ3^{′})\overline{v_{2} · ∇(ψ1^{′} − ψ3^{′})} \\ = & (ψ1^{′} − ψ3^{′})v̄_{2} · ∇(ψ̄_{1} − ψ̄_{3}) + (ψ1^{′} − ψ3^{′})v2^{′} · (ψ1^{′} − ψ3^{′}) \\ + (ψ1^{′} − ψ3^{′})v̄_{2} · ∇(ψ1^{′} − ψ3^{′}) + (ψ1^{′} − ψ3^{′})v2^{′} · ∇(ψ̄_{1} − ψ̄_{3}) \\ − (ψ1^{′} − ψ3^{′}) [v̄_{2} · ∇(ψ̄_{1} − ψ̄_{3}) + v2^{′} · ∇(ψ1^{′} − ψ3^{′})] \\ = & (ψ1^{′} − ψ3^{′}) [v̄_{2} · (ψ1^{′} − ψ3^{′}) + v2^{′} · ∇(ψ̄_{1} − ψ̄_{3})] \end{matrix}  (5.0.13)

となる. したがって, その積分は

\begin{matrix}∫ (ψ1^{′} & − ψ3^{′}) [v_{2} · ∇(ψ_{1} − ψ_{3})]′ dσ \\ = & ∫ (ψ1^{′} − ψ3^{′}) [ū_{2} ∂_{x}(ψ1^{′} − ψ3^{′}) + v2^{′} ∂_{y}(ψ̄_{1} − ψ̄_{3})] dσ \\ = & ∫ ū_{2} [\frac{(ψ1^{′} − ψ3^{′})^{2}}{2} ]_{0}^{L} dy + ∫_{−W}^{W} [∫_{0}^{L}(ψ1^{′} − ψ3^{′})v2^{′}dx] ∂_{y}(ψ̄_{1} − ψ̄_{3}) dy \\ = & ∫ \overline{(ψ1^{′} − ψ3^{′})v2^{′}} ∂_{y}(ψ̄_{1} − ψ̄_{3}) dσ \end{matrix} ; (5.0.14)

と計算される.

なお, 右辺第 3 項の積分については, 平均帯状流のポテンシャルエネルギーの時間発展式を導出する際の右辺第 4 項の積分と同様の方法で計算される.

以上より, 擾乱運動のポテンシャルエネルギーの時間発展の式 (2.3.98) が導出される.

参考文献

[1] Charney, J.G., FJØRTOF, R., von Neumann, J., 1950: Numerical Integration of the Barotropic Vorticity Equation. Tellus, 2, 237-254.

[2] Haltiner, G.J., Williams, R.T., 1980 : Numerical prediction and dynamics meteorology. John Wiley & Sons, 477 pp

[3] Holton, J.R., 2004: An introduction to dynamic meteorology, 4th edn. The Geophysics Series, 88, 535 pp

[4] Lewis, M. J., 1998: Clarifying the Dynamics of the General Circulation: Phillips’s 1956 Experiment. Bulletin of the American Meteorological Society, 79, 39-60.

[5] Pedlosky, J., 1987: Geophysical Fluid Dynamics, 2nd edn. Springer-Verlag, New York, 710 pp

[6] Phlillips, N.A.,1956: The general circulation of the atmosphere: a numerical experiment. Quart. J. Roy. Meteor. Soc., 82, 123-164.

[7] Vallis, K.G., 2006: Atmospheric and Ocean Fluid Dynamics: Fundamentals and Large-scale Circulation. CAMBRIDGE, 745 pp

[8] 小倉義光, 1978: 気象力学通論. 東京大学出版会, 249 pp

[9] 岸保勘三郎, 田中正之, 時岡達志, 1982: 大気大循環. 大気科学講座. 東京大学出版会, 256 pp

[10] 田原智志, 2001: 簡単化したモデルを用いた, 大気大循環に関する数値的研究 [11] 時岡達志, 山岬正紀, 佐藤信夫 : 気象の数値シミュレーション. 気象の教室.

東京大学出版会, 247 pp

[12] 新田尚, 二宮洸三, 山岸米次郎, 2009: 数値予報と現代気象学. 東京堂出版, 224 pp

[13] 増田善信, 1981: 数値予報 -その理論と実際-. 東京堂出版, 278 pp