%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%表題   2 次元非静力学モデル -- 付録 C 雲微物理過程
%
%履歴   2007-06-05 小高正嗣      乱流運動エネルギー方程式の導出を追加
%       2007-06-22 小高正嗣
%       2009-09-07 山下達也      比湿で記述された Kessler パラメタリゼー
%                                ションについての記述を追加
%       2009-09-30 山下達也      雨粒の平均終端速度の導出を追加
%       2009-10-07 山下達也      雨粒の終端速度の式の導出を追加
%       2009-11-12 山下達也      終端速度の式の変更に付随して一部書き換
%                                え
%       2009-12-08 山下達也      一部加筆
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{雲微物理過程}

Kessler(1969) に基づく雲微物理パラメタリゼーションの, 終端速度 $V_D$, 雲
水の衝突併合による雨水比湿の変化率 $CL_{cr}$, 平均終端速度 $V_{term}$,
蒸発による雨水比湿の変化率 $EV_{rv}$ について解説する
\footnote{本章の内容は Ogura and Takahashi (1971), 浅井 (1983) を参考に
した. }. 


\section{雨粒の終端速度}

Newton の抵抗法則より球体の抵抗力 $F_D$ は
\begin{eqnarray}
F_D = \Dinv{8} \rho V_D^2 C_D D^2 \Deqlab{shuutan1}
\end{eqnarray}
と表される. 
ここで $C_D$ は抵力係数であり, 一般にレイノルズ数の関数である. 
雨滴の落下のようにレイノルズ数が大きい現象の場合, レイノルズ数の定義によ
り, 粘性力は流れ場にほとんど寄与しなくなる. 
このとき $C_D$ はレイノルズ数に依存しない定数となる. 
抵抗力と重力の釣合いを考えると
\begin{eqnarray}
\Dinv{6} \pi \rho_w D^3 g = \Dinv{8} \rho V_D^2 C_D D^2 \Deqlab{shuutan2}
\end{eqnarray}
となる. 
\Deqref{shuutan2} を $V_D$ について解くと, 
\begin{eqnarray}
V_D 
 = \left( \frac{ 4 \rho_w g D}{3 C_D \rho} \right)^{1/2}
 \Deqlab{shuutan3}
\end{eqnarray}
となる. 
Kessler(1969) では $\rho_w = 10^3$[kg/m$^3$], $g = 9.8$[m/s$^2$],
$\rho_0 = 1.2$ [kg/m$^3$], $C_D = 0.644$ として
\begin{eqnarray}
V_D 
 = 130 \left( \frac{\rho_0}{\rho} \right)^{0.5} D^{0.5}
 \Deqlab{shuutan4}
\end{eqnarray}
としている\footnote{
Kessler(1969) では $C_D$ をどのように決めたのかについては書かれていない. 
Gunn and Kinzer(1949) によると, レイノルズ数が 3000 程度である雨粒の
$C_D$ の値は 0.66 となるので, Kessler(1969) は系の特徴的なレイノルズ数が
3000 程度であると想定して $C_D$ の値を決めたのかも知れない. 
}. 
但し $\rho_0$ は地表面での大気密度である. 
他の惑星大気においても $C_D = 0.644$ であるとみなすと, 
\begin{equation}
V_D 
 = 1.4389 \left( \frac{\rho_w g D}{\rho} \right)^{1/2}
 \Deqlab{shuutan5}
\end{equation}
が得られる. 

\section{雲水の衝突併合}

雲水の衝突併合による雨水混合比の変化率 $CL_{cr}$ は, 直径 $D$ の単一の
雨粒の衝突併合による質量変化率 $(dm(D)/dt)_{cr}$ と $D$ から $D+dD$ の
範囲の直径を持つ雨粒の数 $N_{D}$ を用いて
\begin{equation}
  CL_{cr} = \Dinv{\rho}\int_{0}^{\infty}\left(\DD{m}{t}\right)_{cr}
   N_{D}\Dd D
  \Deqlab{衝突併合による雨水混合比の変化率}
\end{equation}
と表される. $(dm(D)/dt)_{cr}$ は,
\begin{equation}
  \left(\DD{m}{t}\right)_{cr} = \frac{\pi}{4}D^{2} V_D E \rho q_{c}
  \Deqlab{衝突併合による雨粒の成長方程式}
\end{equation}
と表される. ここで $V_D$ は雨粒の落下速度, $E$ は雨粒と衝突した雲粒
のうち雨粒に併合される割合を表す係数(捕捉係数)である.

雨粒のサイズ分布関数と雨粒の落下速度 $V_D$ を以下のように仮定する.
\begin{eqnarray}
  N_{D} &=& N_{0}\exp (-\lambda D), \Deqlab{マーシャル・パルマー型分布関
   数1} \\
  V_D &=& 1.4389 \left( \frac{\rho_w g D}{\rho} \right)^{1/2}
   \Deqlab{マーシャル・パルマー型分布関数2}.
\end{eqnarray}
ここで $N_{0}, \lambda $ はパラメータである.  式\Deqref{マーシャル・パ
ルマー型分布関数2}の分布は一般にマーシャル・パルマー型分布 (Marshall
and Palmer, 1948) と呼ばれる.  Kessler (1969) では $N_{0}=10^{7}$ とする. 
これを式\Deqref{衝突併合による雨水混合比の変化率}に代入すると,
\begin{eqnarray}
  CL_{cr} &=& \frac{1.4389 \pi}{4} \left( \frac{\rho_w g D}{\rho}
				   \right)^{1/2} E N_{0} q_{c} 
            \int_{0}^{\infty}D^{2.5}\exp(-\lambda D)\Dd D \\ 
          &=& \frac{1.4389 \pi}{4} 
              \left( \frac{\rho_w g D}{\rho} \right)^{1/2} E N_{0} q_{c}
	      \frac{3.75}{\lambda ^{3}}\frac{1}{2} 
             \sqrt{\frac{\pi}{\lambda}}
              \nonumber \\
          &=& \frac{1.4389 \times 3.75 \pi^{3/2}}{8} 
              \left( \frac{\rho_w g D}{\rho} \right)^{1/2}
              E N_{0} q_{c} \lambda ^{-3.5}
  \Deqlab{CL_cr 項} 
\end{eqnarray}
を得る. ここで $E$ は $D$ によらないと仮定した. 
Kessler (1969) では $E=1$ とする. 

雨粒のサイズ分布曲線の傾きを表すパラメータ $\lambda$ は, 以下の式を用
いて雨水比湿 $q_{r}$ で置き換える.
\begin{eqnarray}
  q_{r} &=& \frac{1}{\rho}\int _{0}^{\infty} \rho_{w}
            \frac{\pi}{6} D^{3} 
            N_{D}\Dd D
        \nonumber \\
        &=& \frac{\pi N_{0}\rho_{w}}{6 \rho}\int _{0}^{\infty}
            D^{3}\exp(-\lambda D)\Dd D 
       \nonumber \\
        &=& \frac{\pi N_{0}\rho_{w}}{6 \rho}\frac{6}{\lambda ^{4}}
       \nonumber \\
        &=& \frac{\pi N_{0}\rho_{w}}{\rho}\lambda ^{-4}.
   \Deqlab{λの式}
\end{eqnarray}
ここで $\rho_{w}$ は液相の密度である. 
これを $\lambda$ について解き, 式\Deqref{CL_cr 項}に代入すると,
\begin{eqnarray}
  CL_{cr} 
   &=& \frac{1.4389 \times 3.75 \pi^{3/2}}{8} 
              \left( \frac{\rho_w g D}{\rho} \right)^{1/2}
              E N_{0} q_{c} 
              \left( \frac{\rho q_r}{\pi N_0 \rho_w} \right)^{0.875}
              \nonumber \\
   &=& \frac{1.4389 \times 3.75 \pi^{5/8}}{8} 
        E N_{0}^{1/8} 
        g^{1/2} \left( \frac{\rho}{\rho_w} \right)^{0.375}
        q_{c} q_r^{0.875}
              \nonumber \\
   &=&  10.344
        g^{1/2} \left( \frac{\rho}{\rho_w} \right)^{0.375}
        q_{c} q_r^{0.875}
\end{eqnarray}
となる. 最後の式変形では, $N_{0}=10^{7}, E=1$ を代入した.


\section{平均終端速度}

平均終端速度 $V_{term}$ は雨滴の鉛直フラックス $F_r$, 雨滴密度 $\rho_r$ 
により
\begin{eqnarray}
 V_{term} = \frac{F_r}{\rho_r} \Deqlab{TermVel1}
\end{eqnarray}
と表される. 
$\rho_r$, $F_r$ はそれぞれ以下のように表される. 
\begin{eqnarray}
\rho_r &=& \int_0^{\infty} N_D m dD, \Deqlab{TermVel2} \\
F_r &=& \int_0^{\infty} N_D m V_D dD. \Deqlab{TermVel3}
\end{eqnarray}
ここで $m$ は直径 $D$ の雨滴の質量であり, 
\begin{eqnarray}
m = \frac{\pi}{6} D^3 \rho_w \Deqlab{TermVel4}
\end{eqnarray}
と書ける. 
\Deqref{マーシャル・パルマー型分布関数1}, \Deqref{マーシャル・パルマー型
分布関数2}, \Deqref{TermVel4} を \Deqref{TermVel2}, \Deqref{TermVel3} に
適用すると, 
\begin{eqnarray}
\rho_r 
 &=& \int_0^{\infty} N_0 \exp(- \lambda D) \frac{\pi}{6} \rho_w
  D^3 dD \nonumber \\
 &=& \frac{\pi \rho_w N_0}{\lambda^4} \Deqlab{TermVel5} \\
F_r 
 &=& \int_0^{\infty} N_0 \exp(- \lambda D) \frac{\pi}{6} \rho_w
  D^3 1.4389 \left( \frac{\rho_w g D}{\rho} \right)^{1/2} dD \nonumber \\
 &=& \frac{1.4389 \times 105 \pi^{3/2}}{96}
   \rho_w^{3/2} g^{1/2} \rho^{-1/2} N_0 \lambda^{- 9/2}
  \Deqlab{TermVel6} 
\end{eqnarray}
となる. 
\Deqref{TermVel5} を \Deqref{TermVel6} に代入して $\lambda$ を消去すると, 
\begin{eqnarray}
F_r 
 &=& \frac{1.4389 \times 105 \pi^{3/2}}{96}
   \rho_w^{3/2} g^{1/2} \rho^{-1/2} N_0 
   \left( \frac{\rho_r}{\pi N_0 \rho_w} \right)^{9/8} \nonumber \\
 &=& \frac{1.4389 \times 105 \pi^{3/8}}{96}
   \rho_w^{3/8} g^{1/2} \rho^{-1/2} N_0^{-1/8} \rho_r^{9/8}
      \Deqlab{TermVel7} 
\end{eqnarray}
となる. 
\Deqref{TermVel1} に \Deqref{TermVel7} を代入すると
\begin{eqnarray}
V_{term}
&=& \frac{1.4389 \times 105 \pi^{3/8}}{96}
   \rho_w^{3/8} g^{1/2} \rho^{-1/2} N_0^{-1/8} \rho_r^{1/8} \nonumber \\
&\approx& 0.3224 g^{1/2} \left( \frac{\rho_w}{\rho} \right)^{0.375}
   q_r^{0.125}
 \Deqlab{TermVel8}
\end{eqnarray}
が得られる. 


\section{雨水の蒸発}

蒸発による雨水混合比の変化率 $EV_{rv}$ は, 式\Deqref{衝突併合による雨
水混合比の変化率} と同様に
\begin{equation}
  EV_{rv} = \frac{1}{\rho_{d}}\int _{0}^{\infty}\left(\DD{m}{t}\right)_{ev}
  N_{D}\Dd D
  \Deqlab{蒸発による雨水混合比の増加率}
\end{equation}
と表される. ここで $(dm(D)/dt)_{ev}$ は直径 $D$ の単一の雨粒の蒸発によ
る質量変化率である.

雨水の蒸発は雨粒の表面からの水蒸気の拡散によって律速されると仮定する. 
雨粒周囲の水蒸気フラックスを $F$ とすると, 雨粒の質量の変化率は
\begin{equation}
  \left(\DD{m}{t}\right)_{ev} = - 4\pi r_{d}^{2}F(r_{d})
\end{equation}
と表される. ここで $r$ は雨粒中心からの距離, $r_{d}$ は雨粒の半径で, 
$F$ は
\[
   F = - K_{d}\DD{\rho_{v}}{r}
\]
と表される. $\rho_{v}$ は水蒸気の密度, $K_{d}$ は水蒸気の拡散係数であ
る.  雨粒の周囲では水蒸気フラックスの収束発散はないと仮定すると,
\[
  \frac{1}{r^{2}}\DP{}{r}\left(r^{2}F\right) = 0
\]
が成り立つ. これを積分し 
\[
   \rho _{v} = -\frac{C_{1}}{r} + C_{2}
\]
境界条件 $r=r_{d}$ で $\rho_{v} = \rho_{v,s}$, $r=\infty$ で 
$\rho_{v}=\rho_{v,\infty}$ を適用すると, 
\[
  C_{1} = (\rho _{v,\infty}-\rho_{v,s})r_{d}, 
  \quad C_{2} = \rho _{v,\infty}
\]
これより, 雨粒表面での拡散による水蒸気フラックスは
\begin{eqnarray}
  F(r_{d}) &=& - K_{d}\left.\DD{\rho_{v}}{r}\right|_{r=r_{d}} \nonumber \\
    &=& K_{d}\frac{\rho_{v,s} - \rho _{v,\infty}}{r_{d}}
\end{eqnarray}
よって,
\begin{equation}
  \left(\DD{m}{t}\right)_{ev} = - 4\pi r_{d}K_{d}(\rho_{v,s}-\rho_{v,\infty})
\end{equation}
と表される. 雨粒が落下しながら蒸発する場合には, $K_{d}$ に補正項のついた
\begin{equation}
  \left(\DD{m}{t}\right)_{ev} = - 4\pi r_{d}\left(1+\frac{Fr}{s}\right)
                K_{d}(\rho_{v,s}-\rho_{v,\infty})
  \Deqlab{蒸発による雨粒の成長方程式}
\end{equation}
が用いられる. ここで $F$ は換気因子, $s$ は雨粒表面でのクヌーセン層の
厚さである\footnote{
\Deqref{蒸発による雨粒の成長方程式} は Kinzer and Gunn(1951) で導出され
ている. 
導出方法については要確認である. 
}.

Kessler (1969) では, \Deqref{蒸発による雨粒の成長方程式}の右辺の項を
以下のように近似する.
\begin{eqnarray*}
  4\pi r_{d}\left(1+\frac{Fr}{s}\right) &\sim& 2.24\times 10^{3}D^{1.6}, \\
  K_{d}(\rho _{v,s} - \rho_{v,\infty}) &\sim& 10^{-5}(\rho_{v,s}-\rho_{v,\infty}).
\end{eqnarray*}
このとき\Deqref{蒸発による雨粒の成長方程式}は
\begin{equation}
  \left(\DD{m}{t}\right)_{ev} \sim 
     - 2.24\times 10^{-2}(\rho_{v,s}-\rho_{v,\infty})D^{1.6}
  \Deqlab{近似された蒸発による雨粒の成長方程式}
\end{equation}
となる. これを式\Deqref{蒸発による雨水混合比の増加率}に代入し, 
雨粒のサイズ分布として\Deqref{マーシャル・パルマー型分布関数1}を
仮定すると, 
\begin{eqnarray}
  EV_{rv} &=& - \frac{1}{\rho_{d}}2.24\times 10^{-2}(q_{v,s} -q_{v,\infty})\int _{0}^{\infty}
              D^{1.6}N_{D}\Dd D, \nonumber \\
          &=& - 2.24\times 10^{-2}(q_{v,s} -q_{v,\infty})\int _{0}^{\infty}
              D^{1.6}N_{0}\exp(-\lambda D)\Dd D \nonumber \\
          &=& - 2.24\times 10^{-2}(q_{v,s} -q_{v,\infty})
              N_{0}\frac{\Gamma(2.6)}{\lambda ^{13/5}}  \nonumber \\
          &=& - 2.24\times 10^{-2}
              \Gamma(2.6) (\pi \rho_{w})^{-0.65}
	      N_{0}^{0.35}
	      (q_{v,s} -q_{v,\infty})
	      (\rho_{d}q_{r})^{0.65}
              \nonumber \\
          &=& - 1.7\times 10^{-4}N_{0}^{0.35}(q_{v,s} -q_{v,\infty})
                (\rho_{d}q_{r})^{0.65} \nonumber \\
          &=& - 4.81\times 10^{-2}(q_{v,s} -q_{v,\infty})
                (\rho_{d}q_{r})^{0.65} 
\end{eqnarray}
最後の式変形を行う際には\Deqref{λの式}式の関係を用いて $\lambda$ を消去し,
$\Gamma(2.6)=1.4296245$, $N_{0}=10^{7}$ とした
\footnote{Kessler (1969) では最終的には
\[
  EV_{rv} = 4.85\times 10^{-2}(q_{v,s} -q_{v,\infty})(\rho_{d}q_{r})^{0.65}
\]
としている. }.  


