% 表題 DCPAM5 力学過程 -- 支配方程式系 % % 履歴 %\Drireki{2010/04/15 高橋芳幸} %\Drireki{2009/10/05 高橋芳幸} %\Drireki{2009/03/03 森川靖大} %\Drireki{2008/06/15 森川靖大} %\Drireki{2008/03/25 森川靖大・石渡正樹} %\Drireki{2005/04/04 石渡正樹} %\Drireki{1993/06/16 沼口敦・保坂征宏} % %\section{支配方程式} \section{数理表現} ここでは力学過程の支配方程式系の数理表現を示す. この方程式系の詳細に関しては, Haltiner and Williams (1980) もしくは 別紙『 \htmladdnormallinkfoot{支配方程式系の導出に関する参考資料} {http://www.gfd-dennou.org/library/dcpam/dcpam5/dcpam5\_latest/doc/derivation/pub/derivation.pdf}』 の『力学過程の支配方程式系の導出』を参照せよ. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{連続の式} % \vspace*{-5ex} \begin{align} \Deqlab{連続の式} \DP{\pi}{t} + \Dvect{v}_H \cdot \nabla_{\sigma} \pi &= - D - \DP{\dot{\sigma}}{\sigma}. \end{align} またこれより, 地表気圧の時間変化式と $\sigma$系での鉛直速度$\dot{\sigma}$を求める診断式 % \begin{align} \Deqlab{気圧傾向} \DP{\pi}{t} = - \int_{0}^{1} \Dvect{v}_{H} \cdot \nabla_{\sigma} \pi d \sigma - \int_{0}^{1} D d \sigma , \end{align} % \begin{align} \Deqlab{鉛直速度} \dot{\sigma} = - \sigma \DP{\pi}{t} - \int_{0}^{\sigma} D d \sigma - \int_{0}^{\sigma} \Dvect{v}_{H} \cdot \nabla_{\sigma} \pi d \sigma , \end{align} % が導かれる. さらに, これらより % \begin{align} \omega & \equiv \DD{p}{t} = p_s \left\{ \sigma \Dvect{v}_H \cdot \nabla_\sigma \pi - \int_0^\sigma \left( \Dvect{v}_H \cdot \nabla_\sigma \pi + D \right) d\sigma \right\} \end{align} % が得られる. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{静水圧の式} \vspace*{-5ex} \begin{align} \Deqlab{静水圧の式} \DP{\Phi}{\sigma} = - \frac{RT_v}{\sigma}. \end{align} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{運動方程式} \vspace*{-5ex} \begin{align} \Deqlab{渦度方程式} \DP{\zeta}{t} \ &= \ \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{V_A}{\lambda} - \DP{U_A}{\mu} \right) + {\cal D}(\zeta), \\ \Deqlab{発散方程式} \DP{D}{t} \ &= \ \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{U_A}{\lambda} + \DP{V_A}{\mu} \right) - \nabla^{2}_{\sigma} ( \Phi + R \overline{T} \pi + \mbox{\sl KE} ) + {\cal D}(D). \end{align} % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{熱力学の式} \vspace*{-5ex} \begin{align} \begin{split} \Deqlab{熱力学の式} \DP{T}{t} \ &= \ - \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{U T^{\prime}}{\lambda} + \DP{V T^{\prime}}{\mu} \right) + T^{\prime} D \\ \ & \ \quad - \dot{\sigma} \DP{T}{\sigma} + \kappa T_v \left( \DP{\pi}{t} + \Dvect{v}_H \cdot \Dgrad_{\sigma} \pi + \frac{ \dot{\sigma} }{ \sigma } \right) + \frac{Q}{C_p} + {\cal D}(T) + {\cal D}^{\prime}(\Dvect{v}). \end{split} \end{align} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{水蒸気の式} \vspace*{-5ex} \begin{align} \begin{split} \Deqlab{水蒸気の式} \DP{q}{t} \ &= \ - \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{U q}{\lambda} + \DP{V q}{\mu} \right) + q D \\ \ & \quad - \dot{\sigma} \DP{q}{\sigma} + S_{q} + {\cal D}(q). \end{split} \end{align} ここで, 独立変数は以下の通りである. % \begin{align} \varphi & : \quad 緯度 [\mathrm{deg.}], \\ \lambda & : \quad 経度 [\mathrm{deg.}], \\ \sigma & \equiv p/p_s, \\ t & : \quad 時間 [\mathrm{s}]. \end{align} % ここで, $p$は気圧, $p_s$ は地表面気圧である. また $\mu \equiv \sin \varphi$ である. モデルで時間発展を計算することとなる予報変数は以下の通りである. % \begin{align} \pi\ (\varphi, \lambda) & \equiv \ln p_s, \\ T\ (\varphi, \lambda, \sigma) & : \quad 気温 [\mathrm{K}], \\ q\ (\varphi, \lambda, \sigma) & : \quad 比湿 [\mathrm{kg}\ \mathrm{kg}^{-1}], \\ \Deqlab{渦度定義} \zeta\ (\varphi, \lambda, \sigma) & \equiv \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{V}{\lambda} - \DP{U}{\mu} \right) : \quad 渦度 [\mathrm{s}^{-1}], \\ \Deqlab{発散定義} D\ (\varphi, \lambda, \sigma) & \equiv \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{U}{\lambda} + \DP{V}{\mu} \right) : \quad 発散 [\mathrm{s}^{-1}]. % u\ (\varphi, \lambda, \sigma) % & : \quad 東西風速 [\mathrm{m}\ \mathrm{s}^{-1}], \\ % v\ (\varphi, \lambda, \sigma) % & : \quad 南北風速 [\mathrm{m}\ \mathrm{s}^{-1}]. \end{align} % ここで, % \begin{align} U(\varphi, \lambda, \sigma) &\equiv u(\varphi, \lambda, \sigma) \cos \varphi, \\ V(\varphi, \lambda, \sigma) &\equiv v(\varphi, \lambda, \sigma) \cos \varphi, \\ u & : \mbox{東西風速}, \\ v & : \mbox{南北風速} \end{align} % である. 流線関数$\psi$と速度ポテンシャル$\chi$を導入すると, $U$, $V$, $\zeta$, $D$はそれぞれ以下のように表わされる. % \begin{align} \Deqlab{ψとχによる東西風速の表現} U &= \Dinv{a} \left( \DP{\chi}{\lambda} - (1-\mu^2) \DP{\psi}{\mu} \right), \\ \Deqlab{ψとχによる南北風速の表現} V &= \Dinv{a} \left( \DP{\psi}{\lambda} + (1-\mu^2) \DP{\chi}{\mu} \right), \\ \Deqlab{ψとχによる渦度の表現} \zeta &= \Dlapla \psi, \\ \Deqlab{ψとχによる発散の表現} D &= \Dlapla \chi. \end{align} % 各時間ステップで診断的に求められる変数は以下の通りである. % \begin{align} \Phi & \equiv gz : \quad ジオポテンシャル高度 [\mathrm{m}^{2}\ \mathrm{s}^{-2}], \\ \dot{\sigma} & \equiv \DD{\sigma}{t} \ \equiv \ \DP{\sigma}{t} + \frac{u}{a \cos \varphi} \DP{\sigma}{\lambda} + \frac{v}{a} \DP{\sigma}{\varphi} + \DP{\sigma}{\sigma}, \\ \overline{T}\ (\sigma) & : \quad 基準温度 [\mathrm{K}], \\ T^{\prime}\ (\varphi, \lambda, \sigma) & \equiv T - \overline{T}, \\ T_v\ (\varphi, \lambda, \sigma) & \equiv T \left\{ 1 + \left(\epsilon_v^{-1} - 1\right) q \right\}, \\ T_v^{\prime}\ (\varphi, \lambda, \sigma) & \equiv T_v - \overline{T}, \\ % \theta & \equiv T \left( p/p_{s} \right)^{- \kappa}, \\ \Deqlab{東西運動量移流項} U_A\ (\varphi, \lambda, \sigma) & \equiv ( \zeta + f ) V - \dot{\sigma} \DP{U}{\sigma} - \frac{R T_v^{\prime}}{a} \DP{\pi}{\lambda} + {\cal F}_{\lambda} \cos \varphi, \\ % \Deqlab{南北運動量移流項} V_A\ (\varphi, \lambda, \sigma) & \equiv - ( \zeta + f ) U - \dot{\sigma} \DP{V}{\sigma} - \frac{R T_v^{\prime}}{a} (1-\mu^2) \DP{\pi}{\mu} + {\cal F}_{\varphi} \cos \varphi, \\ % \Deqlab{πの移流} \Dvect{v}_H \cdot \nabla_{\sigma} \pi &\equiv \frac{U}{a (1 - \mu^2)} \DP{\pi}{\lambda} + \frac{V}{a} \DP{\pi}{\mu} \\ % \Deqlab{ラプラシアン} \nabla^{2}_{\sigma} & \equiv \frac{1}{a^{2} (1-\mu^2)} \DP[2]{}{\lambda} + \frac{1}{a^{2}} \DP{}{\mu} \left[ (1-\mu^2) \DP{}{\mu} \right], \\ % \Deqlab{運動エネルギー項} \mbox{\sl KE}\ (\varphi, \lambda, \sigma) & \equiv \frac{U^{2}+V^{2}}{2 (1-\mu^2) } \\ % {\cal D}(\zeta) & : \qquad 渦度の水平拡散とスポンジ層における散逸, \Deqlab{渦度水平拡散} \\ {\cal D}(D) & : \qquad 発散の水平拡散とスポンジ層における散逸, \Deqlab{発散水平拡散} \\ {\cal D}(T) & : \qquad 熱の水平拡散, \Deqlab{温度水平拡散} \\ {\cal D}(q) & : \qquad 水蒸気の水平拡散, \Deqlab{比湿水平拡散} \\ {\cal F}_\lambda \ (\varphi, \lambda, \sigma) & : \qquad 小規模運動過程 (経度方向), \Deqlab{小規模運動過程 (経度方向)} \\ {\cal F}_\varphi \ (\varphi, \lambda, \sigma) & : \qquad 小規模運動過程 (緯度方向), \Deqlab{小規模運動過程 (緯度方向)} \\ Q \ (\varphi, \lambda, \sigma) & : \qquad 放射, 凝結, 小規模運動過程等による加熱・温度変化, \Deqlab{放射, 凝結, 小規模運動過程等による加熱・温度変化} \\ S_q \ (\varphi, \lambda, \sigma) & : \qquad 凝結, 小規模運動過程等による水蒸気ソース, \Deqlab{凝結, 小規模運動過程等による水蒸気ソース} \\ {\cal D}' \ (\Dvect{v}) & : \qquad 摩擦熱 \Deqlab{摩擦熱}. \end{align} % 各水平拡散\Deqref{渦度水平拡散}〜\Deqref{比湿水平拡散} に関しては\ref{波数依存型の水平拡散項}節で説明される. % % 摩擦熱の説明. 物理過程に移される. %\begin{align} % {\cal D}^{\prime} (\Dvect{v}) % = - \frac{1}{C_p} % \Dvect{v} \cdot ( \DP{\Dvect{v}}{t} )_{diff} . %\end{align} % %$( \DP{\Dvect{v}}{t} )_{diff} $ は, %水平および鉛直の拡散による $u,v$ の時間変化項である. % 定数は以下の通りである. % \begin{align} a & : \quad 惑星半径 [\mathrm{m}], \\ R & : \quad 乾燥大気の気体定数 [\mathrm{J\ kg}^{-1}\ \mathrm{K}^{-1}], \\ C_p & : \quad 乾燥大気の大気定圧比熱 [\mathrm{J\ kg}^{-1}\ \mathrm{K}^{-1}], \\ f & : \quad コリオリパラメータ [\mathrm{s}^{-1}], \\ \kappa & \equiv R/C_p, \\ \epsilon_v & : \quad 水蒸気分子量比. % \end{align} % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{境界条件} 鉛直流に関する境界条件は % \begin{align} \dot{\sigma} = 0 \ \ \ at \ \ \sigma = 0 , \ 1 . \end{align} % である. %ただし熱的境界条件については %%別紙『地表面過程』 %\Dchapref{地表面過程} %において記述する. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %\section{水平拡散項} %\label{水平拡散節} %\subsection{波数依存型の水平拡散項} \subsection{水平拡散とスポンジ層} \label{波数依存型の水平拡散項} 水平拡散とスポンジ層における渦度と発散の散逸は次のように表現する. % \begin{align} {\cal D}(\zeta) &= {\cal D_{HD}}(\zeta) + {\cal D_{SL}}(\zeta) \\ {\cal D}(D) &= {\cal D_{HD}}(D) + {\cal D_{SL}}(D) \\ {\cal D}(T) &= {\cal D_{HD}}(T) + {\cal D_{SL}}(T) \\ {\cal D}(q) &= {\cal D_{HD}}(q) \end{align} % ここで, ${\cal D_{HD}}$, ${\cal D_{SL}}$ はそれぞれ水平拡散とスポンジ層における 散逸を表す. 水平拡散項は, 次のように $\nabla^{N_D}$ の形で計算する. % \begin{equation} \label{水平拡散} {\cal D_{HD}}(\zeta) = - K_{HD} \left[ (-1)^{N_D/2} \nabla^{N_D} - \left( \frac{2}{a^2} \right)^{N_D/2} \right] \zeta , \end{equation} % \begin{equation} {\cal D_{HD}}(D) = - K_{HD} \left[ (-1)^{N_D/2} \nabla^{N_D} - \left( \frac{2}{a^2} \right)^{N_D/2} \right] D , \end{equation} % \begin{equation} {\cal D_{HD}}(T) = - (-1)^{N_D/2} K_{HD} \nabla^{N_D} T , \end{equation} % \begin{equation} {\cal D_{HD}}(q) = - (-1)^{N_D/2} K_{HD} \nabla^{N_D} q . \end{equation} % %この水平拡散項は計算の安定化のための意味合いが強い. 小さなスケールに選択的な水平拡散を表すため, 慣例として $N_D$ には 4$\sim$16 を用いることが多い. スポンジ層における運動量の散逸項は, 東西平均成分を減衰させる場合とさせない場合の 2 通りの 計算法を導入する. % 東西平均成分も減衰させる場合には, % \begin{align} {\cal D_{SL}}(\zeta) &= - \gamma_M \zeta, \\ {\cal D_{SL}}(D) &= - \gamma_M D, \end{align} % となる. ここで, $\gamma_M$ はスポンジ層における運動量の減衰係数である. 東西平均成分を減衰させない場合には, % \begin{align} {\cal D_{SL}}(\zeta) &= - \gamma_M ( \zeta - \bar{\zeta} ), \\ {\cal D_{SL}}(D) &= - \gamma_M ( D - \bar{D} ), \end{align} % となる. ここで, $\bar{}$ は, 東西平均を表す. スポンジ層内の温度擾乱の減衰には以下の項を導入する. % \begin{align} {\cal D_{SL}}(T) &= - \gamma_H ( T - \bar{T} ), \end{align} % % ここで, $\gamma_H$ はスポンジ層における温度擾乱の減衰係数である. 減衰係数 $\gamma_M$, $\gamma_H$ の $\sigma$ 依存性に一般形はないが, \Dmodel では 下のような $\sigma$ 依存性を考慮する. % \begin{align} \gamma_M = \left\{ \begin{array}{ll} \gamma_{M,0} \left( \frac{\sigma_0}{\sigma} \right)^{N_{SL}}, & \text{($\sigma \le \sigma_{lim}$)} \\ 0 . & \text{($\sigma > \sigma_{lim}$)} \end{array} \right. \\ \gamma_H = \left\{ \begin{array}{ll} \gamma_{H,0} \left( \frac{\sigma_0}{\sigma} \right)^{N_{SL}}, & \text{($\sigma \le \sigma_{lim}$)} \\ 0 . & \text{($\sigma > \sigma_{lim}$)} \end{array} \right. \end{align} % ここで, $\gamma_{M,0}$, $\gamma_{H,0}$, $N_{SL}$, $\sigma_{lim}$ はそれぞれ, $\sigma = \sigma_0$ における減衰係数, $\sigma$ 依存性の指数, スポンジ層の 下限の $\sigma$ である. \Dmodel では, $\sigma_0$ はモデル最上層の $\sigma$ としている. \subsection{水平拡散係数の値を決めるための判断材料} \label{水平拡散係数の値} 水平拡散係数 $K_{HD}$ の値は問題に応じて試行錯誤して決めることになる. 判断規準の 1 つは, エネルギースペクトル \begin{eqnarray} \mathcal{E}_n = \frac{1}{4} \frac{a^2}{n(n+1)} \sum^{n}_{m=-n} (|\tilde{\zeta}^{m}_{n}|^2 + |\tilde{D}^{m}_{n}|^2) \end{eqnarray} において (このエネルギースペクトルの表式は Koshyk and Hamiltion, 2001 による), 高波数領域におけるエネルギーの急激な減衰やエネルギーの 蓄積が現れないようにするということである. 考え方の一つとして, 高波数領域におけるエネルギースペクトルが $n^{-5/3}$ の 指数則に従っていれば良かろう, と判断することもできるだろう. Takahashi et al. (2006) は AFES を用いて地球大気に関する高分解能計算を 行い, 200hPa における運動エネルギーのスペクトルが 低波数領域 (およそ $n<80$) では $n^{-3}$ の 指数則に, 高波数領域では $n^{-5/3}$ 則に従う 水平拡散係数の値を決定している. かれらは水平拡散のオペレータとして $K_{HD} \nabla^4$ を用いた場合 ($N_D=4$ とした場合に対応する) について, T79L24, T159L24, T319L24, T639L24, T639L48 と分解能を変更した 計算を行った. その結果, 拡散係数の値としては \begin{eqnarray} K_H = 1.2 \times 10^{21} n_t^{-3.22} \ \ \ [{\rm m}^4~{\rm sec}^{-1}] \end{eqnarray} を与えるのが良いとしている. ここで, $n_t$ は切断波数である. この式から, 各種の水平分解能に応じて $K_{HD}$ の値とモデルで表現される 最小スケールの減衰率を計算すると \Dtabref{水平拡散係数} となる. \begin{table}[h] \begin{center} \begin{tabular}{c|c|c|c} 切断波数 & $K_H$ (m$^4$ sec$^{-1}$) & 減衰率 (1/days) & 減衰時間 (day) \\ \hline T21 & $7 \times 10^{16}$ & 0.7 & 1.4 \\ T42 & $7 \times 10^{15}$ & 1.2 & 0.8 \\ T79 & $9 \times 10^{14}$ & 1.9 & 0.5 \\ T159 & $1 \times 10^{14}$ & 3.2 & 0.3 \\ T319 & $1 \times 10^{13}$ & 5.6 & 0.2 \\ T639 & $1 \times 10^{12}$ & 9.6 & 0.1 \\ \end{tabular} \end{center} \caption{Takahashi et al. (2006) の結果から得られた水平拡散係数の値と 最小スケールの減衰率. $N_D=4$ とした場合の結果を示す. ただし, T21 と T42 の計算は Takahashi et al. (2006) ではなされ ていない. 減衰率は ${\displaystyle K_{HD} \left\{ \frac{n(n+1)}{a^2} \right\}^2}$ を用いて計算した. } \Dtablab{水平拡散係数} \end{table} 更に, Takahashi et al. (2006) は, 水蒸気無し・地形無しの設定の もとで Held and Suarez 実験 (Held and Suarez, 1994) も行っている. この場合, T639L24 とした時に 高波数域におけるエネルギースペクトルが指数則に従うようにするためには 上式で与えられる $K_{HD}$ の値の 0.5 倍を使うのが良いという結果を得た. このように, エネルギースペクトルが指数則に従うようにするためには, 分解能と実験設定に応じて試行錯誤で拡散係数を決定する必要がある. その際には, 上記の数値を目安として用いるのが良いだろう. \subsection{発散減衰} 下のような発散減衰 (divergence damping) を導入する. % \begin{align} {\cal D}_{dd} &= - \gamma D \end{align} % ここで, $\gamma$ は減衰係数である. \footnote{ 以下は, 本来別のところに書くべきこと. \\ 任意の初期値から始める計算において, 積分初期の急激な変化を抑制することを 目的として, $\gamma$ は, 実際には下のように与える. % \begin{align} \gamma = \frac{1}{\Delta t} \frac{T_{dd} - t}{T_{dd}} \end{align} % ここで, $\Delta t$ は時間ステップ, $T_{dd}$ は計算開始から発散減衰が 適応される時間, $t$ は時刻である. } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %\subsection{波数非依存型の水平拡散項} % % 未導入のため, 導入されるまでコメントアウト. % %水平拡散を波数に依存しない一様な値にすることもできる. %詳細省略.