修士論文発表会

非静力学対流モデルの定式化に関する再検討

地形を考慮した雲解像モデルの開発へ向けて

Reconsideration of the formulation of a nonhydrostatic convection model

Toward a cloud resolving model with terrain-following coordinates

高橋こう子

北海道大学大学院理学研究科地球惑星科学専攻地球惑星流体科学講座地球流体力学研究室

2004年2月4日

1. 背景

最近活発な研究が行われているメソ気象

メソ気象研究の数値実験によく利用される非静力学数値モデル

既存の数値モデルを利用する場合の問題点

• 完成度が高い →不必要な項目が入っている可能性がある

• 数理モデルで使われているパラメタリゼーションが明記されていない場合 →解析する際, 何の効果が効いているのか判別しにくい

• 多くの既存の数値モデルは複雑なため, 解読が困難 →新しい概念を導入することが難しい

→ [figure: 自分で数値モデルを開発する]

2. 研究目的

本研究では

1. 準圧縮モデルの定式化を示す.

2. 準圧縮モデルの離散化を示す.

• 音波減衰項の係数の具体的な値を決める.

• 数値粘性項の係数の具体的な値を決める.

3. 簡単な対流を表現できる2次元の数値モデルを作成する.

を行い,

→ [figure: 地形を考慮した雲解像モデル開発の基礎とする]

3. 準圧縮モデルとは

準圧縮モデル

• 多くのメソ気象研究の数値実験で用いられている非静力学モデル

多くの準圧縮モデルで基にされている数理モデル

• Klemp and Wilhelmson (1978) (以下 KW1978)

準圧縮モデルの特徴

• 鉛直方向の運動方程式を静水圧平衡の式で近似しない.

• 連続の式を線形化する(密度の時間変化項を含む).

\frac{∂ρ}{∂t} + \frac{∂(ρ_{0}u)}{∂x} + \frac{∂(ρ_{0}v)}{∂y} + \frac{∂(ρ_{0}w)}{∂z} = 0  (1)

様かつ静止状態.

[figure]

{以下では, 準圧縮の標準的なモデルである KW1978 に \atop 従って定式化を行う.}

[figure]

u, v, w は x, y, z 方向の速度, ρ は密度, ρ_{0} = ρ_{0}(x, y, z) は基本場の密度. 基本場は水平一

4. 乾燥大気における準圧縮モデルの定式化

予報変数: u, v, w (x, y, z 方向の速度), π(≡ (p/p_{0})^{R_{d}/c_{p}} , エクスナー関数),

θ(≡ T (p_{0}/p)^{R_{d}/c_{p}} , 温位).

運動方程式:

\frac{du_{i}}{dt} + c_{p}Θ\frac{∂π}{∂x_{i}} = δ_{i3}g [\frac{θ}{Θ} − 1] + D_{u_{i}},  (2)

圧力方程式:

\frac{∂π}{∂t} + \frac{c_{0}^{2}}{c_{p}ρ_{0}Θ^{2}} \frac{∂}{∂x_{j}} (ρ_{0}Θu_{j} ) = 0,  (3)

熱力学の式:

\frac{dθ}{dt} = D_{θ}.  (4)

D_{u_{j}}, D_{θ} はサブグリッドスケールの乱流に伴う拡散項, c_{0} = c_{p}R_{d}ΠΘ/c_{v} は基本場の音速, c_{p} は定圧比熱, c_{v} は定積比熱, R_{d} 乾燥空気の気体定数, g は重力加速度, Θ は基本場の温位, Π は基本場のエクスナー関数. 基本場は静水圧平衡が成り立つ.

5. サブグリッドスケールの拡散

渦粘性モデルを用いて, 基礎方程式系を以下のように近似する:

D_{u_{i}} = −\frac{∂}{∂x_{j}} u′_{i}u′_{j} ∼ −\frac{∂}{∂x_{j}} \{−K_{m} (\frac{∂u_{j}}{∂x_{i}} + \frac{∂u_{i}}{∂x_{j}} ) + \frac{2}{3}δ_{ij}E\} ,  (5)
D_{θ} = −\frac{∂}{∂x_{j}} u′_{j}θ′ ∼ −\frac{∂}{∂x_{j}} (K_{h}\frac{∂θ}{∂x_{j}} ) .  (6)

ϕ はグリッドスケールの平均成分, ϕ′ はサブグリッドスケールの変動成分, K_{m} は運動量渦拡散係数, K_{h} は熱渦拡散係数, E は乱流運動エネルギー.

\frac{dK_{m}}{dt} = −\frac{3g}{θ̄} \frac{∂θ}{∂z} + 2C_{m}^{2}l^{2} (\frac{∂u_{j}}{∂x_{j}} )^{2} + \frac{1}{2} |ε_{ijk}|C_{m}^{2}l^{2} (\frac{∂u_{j}}{∂x_{i}} + \frac{∂u_{i}}{∂x_{j}} )^{2}

−\frac{1}{3}δ_{ij}K_{m} (\frac{∂u_{j}}{∂x_{j}} ) + \frac{1}{2}δ_{ij} (\frac{∂^{2}K_{m}^{2}}{∂x_{i}∂x_{j}} ) + δ_{ij} (\frac{∂K_{m}}{∂x_{i}}) (\frac{∂K_{m}}{∂x_{j}}) − (\frac{C_{ε}}{2C_{m}l^{2}}) K_{m}^{2}, (7)
E = (\frac{K_{m}}{C_{m}l})^{2} ,  (8)
K_{h} = 3K_{m}.  (9)

l = (ΔxΔyΔz)^{1/3} は混合距離, C_{ε} = C_{m} = 0.2 (Deardorff, 1975).

6. 空間微分の離散化空間微分

• 中心差分で離散化

• 運動方程式中に音波減衰項 (Skamarock and Klemp, 1992)を追加

αDiv = α (\frac{∂u}{∂x} + \frac{∂v}{∂y} + \frac{∂w}{∂z}) ,  (10)
α ≤ 5 × 10^{−6} × \frac{Δx^{2}}{Δτ}  (11)

\begin{matrix}目的 & 1.音波の計算 & {2.条件(11)を満たす音波減衰項 \atop の係数αを決める} \\ 解く方程式 & 運動方程式(2), 圧力方程式(3) & 運動方程式(2), 圧力方程式(3) \\ 音波減衰項 & なし & あり \\ {考慮しない項 \atop (上記以外)} & 移流項,浮力項,拡散項 & 移流項,浮力項,拡散項 \\ 基本場 & 等温(300 K),等圧大気 & 等温(300 K),静水圧平衡 \\ 初期擾乱 & (数値実験1の結果の節で示す) & 2 次元ガウス型の圧力 \\ その他の設定 & Δτ = 1.0 s & \begin{matrix}Δτ = 1.0 s, \\ α ≤ 2.5 & (12) \end{matrix} \end{matrix}

6-1. 音波の計算

実験: case 1-1: Δτ = 1 s, x 方向一様なガウス型圧力の初期擾乱. case 1-2: Δτ = 1 s, z 方向一様なガウス型圧力の初期擾乱.

case 1-1 case 1-2

結果: case 1-1 より, 音速が 333 m/s (≃ 5.0 × 10^{3} m/15 s). case 1-2 より, 音速が 343 m/s (≃ 1.2 × 10^{4} m/35 s).

→[figure: 音波が計算できる.]

6-2. 音波減衰項の係数を決める

実験: case 2-1 : α = 0 (音波減衰項なし).

case 2-2 : α = 2.5 × 10^{−2}.

case 2-1 case 2-2

結果: case 2-1 より, 時間の経過ともにエクスナー関数の偏差が急増. case 2-2 より, 時間が十分経過してもエクスナー関数の偏差が安定. → [figure: 音波減衰項の係数は α = 2.5 × 10^{−2}.]

7. 時間積分の離散化時間積分

• 音波に関する特別な工夫として,時間分割法と陰解法を用いる

– 音波に関係する項を短い時間間隔Δτ で解く ∗ 水平陽解-鉛直陰解法: 水平方向を前進差分, 鉛直方向を後退差分

– その他の項を長い時間間隔Δt で解く: leap-frog 法

• 運動方程式と熱力学の式に数値粘性項を追加

Diff.ϕ = ν_{h} [\frac{∂^{2}ϕ}{∂x^{2}} + \frac{∂^{2}ϕ}{∂y^{2}}] + ν_{v} \frac{∂^{2}ϕ}{∂z^{2}},  (13)
ν_{h} ≤ \frac{1}{8}\frac{ΔxΔy}{Δt}, ν_{v} ≤ \frac{1}{8}\frac{Δz^{2}}{Δt}  (14)
目的条件(14)を満たす数値粘性の係数ν_{h}, ν_{v} を決める
解く方程式\frac{∂ϕ}{∂t} + u_{j} \frac{∂ϕ}{∂x_{j}} + Diff.ϕ = 0
数値粘性項あり
考慮しない項(上記以外) 拡散項
初期擾乱  2次元ガウス型の温度
その他の設定\begin{matrix}速度一定(u = 10 m/s, w = 0 m/s). Δt = 1.0 s, \\ ν_{h} ≤ 1.25 × 10^{5}, ν_{v} ≤ 3.125 × 10^{4} & (15) \end{matrix}

7-1. 数値粘性項の係数を決める

実験: case 3-1 : ν_{h} = 0, ν_{v} = 0(数値粘性項なし).

case 3-2 : ν_{h} = 5.0, ν_{v} = 1.25. case 3-3 : ν_{h} = 5.0 × 10^{2}, ν_{v} = 1.25 × 10^{2} . case 3-4 : ν_{h} = 5.0 × 10^{3}, ν_{v} = 1.25 × 10^{3} .

case 3-1  case 3-2  case 3-3  case 3-4

各図中において, 山型の曲線は左から t = 0, 250, 500, 750, 1000 s での温位偏差の伝播を表す.

結果: 温位偏差の後ろに生じた擾乱が実験結果の中でもっとも小さいのは case 3-2. → [figure: 数値粘性項の係数は ν_{h} = 5.0, ν_{v} = 1.25.]

8. 数値モデルへの移流非線形項の組み込み 1

• x, z 方向の 2 次元

• 水平計算領域 2.4 ×10^{4} m

• 鉛直計算領域 1.0 ×10^{4} m

• 水平格子間隔 Δx = 1.0 × 10^{3} m

• 鉛直格子間隔 Δz = 5.0 × 10^{2} m

\begin{matrix}目的 & 簡単な対流であるサーマルの計算 \\ 解く方程式 & 運動方程式(2), 圧力方程式(3), 熱力学の式(4) \\ 音波減衰項 & あり \\ 数値粘性項 & あり \\ 基本場 & {温度は高度z = 5.0 × 10^{3} mまでは乾燥断熱減率で \atop 減少, 高度z = 5.0 × 10^{3} mから等温.} \\ 初期擾乱 & 2次元ガウス型の温度 \\ その他の設定 & {Δτ = 1.0 s, Δt = 10 s. α, ν_{h}, ν_{v} は先に決めた値, \atop K_{m} は定数.} \end{matrix}

8-1. 数値モデルへの移流非線形項の組み込み 2

実験: case 4-1 : K_{m} = 10 m^{2}/s. case 4-2 : K_{m} = 50 m^{2}/s.

case 4-3 : K_{m} = 100 m^{2}/s.

case 4-1 : t = 0 s  t = 500 s  t = 800 s  t = 1000 s

case 4-2 : t = 0 s  t = 500 s  t = 800 s  t = 1000 s

8-2. 数値モデルへの移流非線形項の組み込み 3

case 4-3 : t = 0 s  t = 500 s  t = 800 s  t = 1000 s

結果: case 4-1〜4-3より, 下部境界の中心に与えたガウス型の温度擾乱は圏界面までサーマルが上昇し, 圏界面で横に広がっている. 運動量渦拡散係数を変えてもほぼ同様.

→ [figure: サーマルを表現できる.]

9. 数値モデルへの渦粘性モデルの組み込み (未完成)

K_{m} 一定にした数値モデルに, 運動量渦拡散係数の時間発展方程式(7)を導入

→ [figure: 乾燥大気における準圧縮の数値モデルの完成](未完成) KW1978 中の離散化において, 運動量渦拡散係数の式(7)中の2項が抜けている:

• 式(7)の第4項

\frac{1}{3}δ_{ij}K_{m} (\frac{∂u_{j}}{∂x_{j}} )  (16)

• 式(7)の第6項

δ_{ij} (\frac{∂K_{m}}{∂x_{i}}) (\frac{∂K_{m}}{∂x_{j}})  (17)

これらの項の効果を調べるためにも,

→ 式(7)を導入した数値モデルが必要 10. まとめ

1. 準圧縮モデルの定式化を示した.

2. 準圧縮モデルの離散化を示した.

• 音波減衰項の係数の具体的な値を決めた.

• 数値粘性項の係数の具体的な値を決めた.

3. K_{m} 固定で,簡単な対流のサーマルを表現できる2次元非静力学対流数値モデルを作成した.

運動量渦拡散係数の時間発展方程式(7)を導入すると,

乾燥大気における準圧縮モデルが完成する.

[figure]

地形を考慮した準圧縮の雲解像モデルができる.^{さらに, 湿潤大気と地形を考慮すると,}

[figure]

11. 参考文献

• Klemp J. B. and Robert B. Wilhelmson, 1978: The simulation of three-dimensional convective storm dynamics. J. Atmos. Sci., 35, 1070-1096.

• Deardorff, J. W., 1975: The development of boundary-layer turbulence models for use in studying the severe storm environment. Proc. SESAME Meeting, Boulder, NOAA-ERL, 251-261.

• Skamarock, Wiliam C. and Joseph B. Klemp, 1992: The stability of time-split methods for the hydrostatic and the nonhydrostatic elastic equation. Mon. Wea. Rev., 120, 2109-2127.

• 坪木和久, 榊原篤志, 2001: CReSS ユーザーズガイド第 2 版. http://www.tokyo.rist.or.jp/CReSS Fujin/CReSS.top.html