乱流の運動方程式、エクマン流

　幼少期から気象や気候に関心があり、地図帳で世界各地の気象要素の平年値を面白がって見ていました。大学卒業後、２度目の受験で気象予報士の資格を取得しましたが、その準備として「基礎気象学」（朝倉書店）を１年程熟読し、気象学を概観していました。「基礎」と冠してあるものの、気象学の広範な理論を要約しているだけあり、読みこなすのに苦労しました。当書の中で、理解に時間を要した箇所を数回にわたり解説していこうと思います。
　今回は、前回の記事（三角関数と$${i}$$の関係）に関連した、エクマンスパイラルの導出について解説します。

流体の運動方程式

　地表面を平面と近似した上で、子午線方向を$${y}$$軸（北極方向が正の向き）とし、$${x}$$軸（東向きが正の向き）を$${y}$$軸に垂直にとります。また、鉛直線を$${z}$$軸（上向きが正の向き）とします。この座標系に関する、ある微小空気塊の時刻$${t}$$における速度を$${(u,v,w)}$$とすると、$${u,v,w}$$は、位置座標$${x,y,z}$$の関数ですが、$${x,y,z}$$はともに$${t}$$の関数であるから、何れも2回連続微分可能だとすれば、$${u,v,w}$$も何れも$${t}$$の2回連続微分可能関数であり、例えば空気塊の$${x}$$軸方向の加速度は、

$$
\begin{aligned}
\frac{\,\mathrm{d}u\,}{\mathrm{d}t}&=\lim_{\Delta t\to0}{\frac{\,u(t+\Delta t)-u(t)\,}{\Delta t}}\\
&=\lim_{\Delta t\to0}{\frac{\,u(x+\Delta x,y+\Delta y,z+\Delta z,t+\Delta t)-u(t)\,}{\Delta t}}\\
&=\lim_{\Delta t\to0}{\frac{\,\{u(x+\Delta x,y+\Delta y,z+\Delta z,t+\Delta t)-u(x,y+\Delta y,z+\Delta z,t+\Delta t)\}+\{u(x,y+\Delta y,z+\Delta z,t+\Delta t)-u(t)\}\,}{\Delta t}}\\
&=\lim_{\Delta t\to0}{\frac{\,\Delta x\,}{\Delta t}\cdot\frac{\,u(x+\Delta x,y+\Delta y,z+\Delta z,t+\Delta t)-u(x,y+\Delta y,z+\Delta z,t+\Delta t)\,}{\Delta x}}+\lim_{\Delta t\to0}{\frac{\,u(x,y+\Delta y,z+\Delta z,t+\Delta t)-u(t)\,}{\Delta t}}\\
&=u\frac{\partial u}{\,\partial x\,}+\lim_{\Delta t\to0}{\frac{\,u(x,y+\Delta y,z+\Delta z,t+\Delta t)-u(t)\,}{\Delta t}}\\
&=\cdots=\frac{\partial u}{\,\partial t\,}+u\frac{\partial u}{\,\partial x\,}+v\frac{\partial u}{\,\partial y\,}+w\frac{\partial u}{\,\partial z\,}
\end{aligned}
$$

と表されます。

　$${\frac{\,\mathrm{d}u\,}{\mathrm{d}t}}$$を、時刻$${t}$$に関する$${u}$$の物質微分といいます。

　すると、$${x}$$方向の運動方程式は$${f}$$、$${p}$$、$${\rho}$$を各々コリオリパラメータ（$${f=2\Omega\sin\varphi}$$（$${\Omega}$$は地球の自転角速度、$${\varphi}$$は観測地の緯度））、気圧、空気塊の密度として、

$$
\rho\frac{\,\mathrm{d}u\,}{\mathrm{d}t}=\rho fv-\frac{\partial p}{\,\partial x\,}
$$

即ち

$$
\frac{\partial u}{\,\partial t\,}+u\frac{\partial u}{\,\partial x\,}+v\frac{\partial u}{\,\partial y\,}+w\frac{\partial u}{\,\partial z\,}=fv-\frac{1}{\,\rho\,}\frac{\partial p}{\,\partial x\,}
$$

と表されます。

乱流の運動方程式

　この運動方程式を、実際の大気の運動に即するよう、乱流成分を加味した運動方程式に書き換えます。その際、大気が非圧縮性流体である（つまり、湧き出しが0：$${\frac{\partial u}{\,\partial x\,}+\frac{\partial v}{\,\partial y\,}+\frac{\partial w}{\,\partial z\,}=0}$$）と仮定します。

　風速の各成分を、次のように平均成分$${\overline{u},\overline{w},\overline{w}}$$と乱流成分$${u',v',w'}$$に分解します。
　$${\overline{u}}$$は、時刻$${t}$$を中心とする、気象学的特徴を相殺しない程度の時間間隔$${\Delta t}$$にわたる$${u}$$の平均です：

$$
\overline{u}=\frac{1}{\,\Delta t\,}\int_{t-\frac{\,\Delta t\,}{2}}^{t+\frac{\,\Delta t\,}{2}}u\,dt
$$

　$${u}$$の、この平均からのずれを$${u'}$$と書きます：

$$
u'=u-\overline{u}
$$

　平均成分の定義により、$${\overline{u'}=0}$$が成り立ちます。

　$${x}$$成分$${u}$$のみ例に挙げましたが、他の成分に関しても同様に定義します。すると、乱流成分の平均が0であることより、例えば

$$
\begin{aligned}
\overline{uv}&=\overline{(\overline{u}+u')(\overline{v}+v')}\\
&=\overline{u}\,\overline{v}+\overline{u'v'}
\end{aligned}
$$

です。同様に、$${\overline{uw}=\overline{u}\,\overline{w}+\overline{u'w'},\overline{u^2}=\overline{u}^2+\overline{u'^2}}$$が成り立ちます。

　前項で導出した運動方程式の両辺の平均をとることで、乱流の運動方程式を得られるのですが、直前に示した3式を使えるよう、運動方程式の左辺各項の係数を偏微分の中に繰り込みます。そのために、本項冒頭で仮定した大気の非圧縮性を用います。

$$
\frac{\partial u}{\,\partial t\,}+u\frac{\partial u}{\,\partial x\,}+v\frac{\partial u}{\,\partial y\,}+w\frac{\partial u}{\,\partial z\,}=fv-\frac{1}{\,\rho\,}\frac{\partial p}{\,\partial x\,}
$$

に$${u\left(\frac{\partial u}{\,\partial x\,}+\frac{\partial v}{\,\partial y\,}+\frac{\partial w}{\,\partial z\,}\right)=0}$$を辺々加え、積の微分による展開を逆向きに見ることで、

$$
\frac{\partial u}{\,\partial t\,}+\frac{\partial (u^2)}{\,\partial x\,}+\frac{\partial (uv)}{\,\partial y\,}+\frac{\partial (uw)}{\,\partial z\,}=fv-\frac{1}{\,\rho\,}\frac{\partial p}{\,\partial x\,}
$$

を得ます。

　両辺の平均をとり、

$$
\frac{\partial \overline{u}}{\,\partial t\,}+\frac{\partial \overline{u^2}}{\,\partial x\,}+\frac{\partial \overline{uv}}{\,\partial y\,}+\frac{\partial \overline{uw}}{\,\partial z\,}=f\overline{v}-\frac{1}{\,\rho\,}\frac{\partial \overline{p}}{\,\partial x\,}
$$

となります。

　$${\overline{uv}=\overline{u}\,\overline{v}+\overline{u'v'}}$$などを代入して得られる

$$
\frac{\partial \overline{u}}{\,\partial t\,}+\frac{\partial \overline{u}^2}{\,\partial x\,}+\frac{\partial (\overline{u}\,\overline{v})}{\,\partial y\,}+\frac{\partial (\overline{u}\,\overline{w})}{\,\partial z\,}=f\overline{v}-\frac{1}{\,\rho\,}\frac{\partial \overline{p}}{\,\partial x\,}-\left(\frac{\partial \overline{u'^2}}{\,\partial x\,}+\frac{\partial \overline{u'v'}}{\,\partial y\,}+\frac{\partial \overline{u'w'}}{\,\partial z\,}\right)
$$

の左辺を展開し、$${\frac{\partial \overline{u}}{\,\partial x\,}+\frac{\partial \overline{v}}{\,\partial y\,}+\frac{\partial \overline{w}}{\,\partial z\,}=0}$$を用いると、乱流の運動方程式

$$
\frac{\partial \overline{u}}{\,\partial t\,}+\overline{u}\frac{\partial \overline{u}}{\,\partial x\,}+\overline{v}\frac{\partial \overline{u}}{\,\partial y\,}+\overline{w}\frac{\partial \overline{u}}{\,\partial z\,}-f\overline{v}=\frac{1}{\,\rho\,}\left(-\frac{\partial \overline{p}}{\,\partial x\,}+\frac{\partial \left(-\rho\overline{u'u'}\right)
}{\,\partial x\,}+\frac{\partial \left(-\rho\overline{u'v'}\right)
}{\,\partial y\,}+\frac{\partial \left(-\rho\overline{u'w'}\right)
}{\,\partial z\,}\right)
$$

が得られます。右辺で偏微分されている、$${-\rho}$$で始まる3項は、平均からのずれにより生じる応力テンソルで、乱渦応力やレイノルズ応力と呼ばれます。これらを順に$${\tau_{xx},\tau_{yx},\tau_{zx}}$$と書きます。

　乱流の運動方程式の、レイノルズ応力に起因する項は、次のように書き換えられます：
　例えば、

$$
\tau_{zx}=\mu\frac{\partial \overline{u}}{\,\partial z\,}
$$

のように、レイノルズ応力は平均流のシアーに比例することが知られています。乱渦の等方性を仮定し、$${\mu}$$が添字に依らないとすると、この関係を用いて

$$
\frac{\partial \overline{u}}{\,\partial t\,}+\overline{u}\frac{\partial \overline{u}}{\,\partial x\,}+\overline{v}\frac{\partial \overline{u}}{\,\partial y\,}+\overline{w}\frac{\partial \overline{u}}{\,\partial z\,}-f\overline{v}=-\frac{1}{\,\rho\,}\frac{\partial \overline{p}}{\,\partial x\,}+\frac{\mu}{\,\rho\,}\left(\frac{\partial^2 \overline{u}
}{\,\partial x^2\,}+\frac{\partial^2 \overline{u}
}{\,\partial y^2\,}+\frac{\partial^2 \overline{u}
}{\,\partial z^2\,}\right)
$$

と、乱流成分が現れない形にできます。

エクマン流

　対流圏は、地表からの高度が概ね100m以下の接地境界層、1㎞以上の自由大気層、それらの間のエクマン層に大別され、層毎に大気運動の特徴が異なります。接地境界層では地面の摩擦の影響が支配的で、風速が高度$${z}$$の対数に比例します。地面の摩擦の影響は$${z}$$の増加に伴い徐々に弱まり、自由大気層では気圧傾度と転向力の釣り合いによって生じる地衡風で大気運動が近似されます。地面の摩擦の影響が残るエクマン層では、以下で示すように、風向が$${z}$$の増加に伴い規則的に変化します。

　簡単のため、エクマン層の下面を高度$${z}$$の基準にとり、平均流を考えます（以下、平均流のバーを省略します）。大気場が次の仮定を満たしているものとします：
・等圧線は$${x}$$軸に平行で、気圧傾度は$${y}$$のみの関数である（高度に依らない）
・風速は$${z}$$のみの関数で、$${w=0}$$である（大気運動が水平方向に一様で、時刻に依らない）
・$${z=0}$$において、風向$${(u,v)}$$と乱渦応力の向きが一致する
・$${(u,v)}$$は$${z\to\infty}$$（エクマン層上面に相当する）のとき地衡風$${(u_g,0)}$$に収束する

　下の2条件は、微分方程式の初期条件です。上の2条件を用いると、前項で得た平均流の運動方程式の$${x}$$成分と$${y}$$成分は、順に次のようになります（$${\nu=\frac{\mu}{\,\rho\,}}$$とします）：

$$
\begin{cases}
-fv=\nu\dfrac{\partial^2 u}{\,\partial z^2\,}\\[3pt]
fu=\dfrac{1}{\,\rho\,}\dfrac{\partial p}{\,\partial y\,}+\nu\dfrac{\partial^2 v}{\,\partial z^2\,}
\end{cases}
$$

　定常流を仮定しているので、気塊に働く気圧傾度力と転向力は釣り合います。前者は高度に依らないとしているので、$${z\to\infty}$$の極限を考えると、$${y}$$方向の釣り合いより、

$$
-\dfrac{\partial p}{\,\partial y\,}-\rho fu_g=0
$$

です。これを代入し、$${(u,v)}$$が満たすべき微分方程式

$$
\begin{cases}
\nu\dfrac{\partial^2 u}{\,\partial z^2\,}=-fv\\
\nu\dfrac{\partial^2 v}{\,\partial z^2\,}=f(u-u_g)
\end{cases}
$$

を得ます。
求解の都合上、第1式の$${u}$$を$${u-u_g}$$に置き換えます：

$$
\begin{cases}
\dfrac{\partial^2 (u-u_g)}{\,\partial z^2\,}=-\dfrac{f}{\,\nu\,}v\\
\dfrac{\partial^2 v}{\,\partial z^2\,}=\dfrac{f}{\,\nu\,}(u-u_g)
\end{cases}
\tag{$*$}
$$

　以下、これを初期条件の下で解きます。
　微分方程式の両辺に現れる（実数値）未知関数と、符号の変化に着目すると、

$$
(*)\Longleftrightarrow
\dfrac{\partial^2}{\,\partial z^2\,}(u-u_g+iv)=\dfrac{if}{\,\nu\,}(u-u_g+iv)
$$

です。

　$${\dfrac{if}{\,\nu\,}}$$の平方根は

$$
\pm\sqrt{\dfrac{f}{\,\nu\,}}\exp{\left(\dfrac{\pi}{\,4\,}i\right)}=\pm\sqrt{\dfrac{f}{\,2\nu\,}}(1+i)
$$

なので、$${(*)}$$の解は

$$
\exp{\left(\pm\sqrt{\frac{f}{\,2\nu\,}}(1+i)z\right)}
$$

の$${\mathbb{C}}$$上の1次結合です：

$$
(*)\Longleftrightarrow
\exist A,B\in\mathbb{C}\;\mathrm{s.t.}\;
u-u_g+iv=A\exp{\left(\sqrt{\frac{f}{\,2\nu\,}}(1+i)z\right)}+B\exp{\left(-\sqrt{\frac{f}{\,2\nu\,}}(1+i)z\right)}
$$

　2つの初期条件「$${z=0}$$において、風向$${(u,v)}$$と乱渦応力の向きが一致する」、「$${(u,v)}$$は$${z\to\infty}$$のとき地衡風$${(u_g,0)}$$に収束する」を用いて、$${A,B}$$の値を決定します。

　まず、後者の極限値が存在するためには$${A=0}$$が必要十分です。一方、前者は、$${z=0}$$において、

$$
\begin{aligned}
\tau_{zx}+i\tau_{zy}=\mu\dfrac{\partial}{\,\partial z\,}(u+iv)
\end{aligned}
$$

が$${u+iv}$$の実数倍であることと同値です。

　$${z=0}$$における$${u+iv}$$の値を$${u_0+iv_0}$$とすると、

$$
B=u_0+iv_0-u_g
$$

であり、従って

$$
\begin{aligned}
\left(\dfrac{\partial}{\,\partial z\,}(u+iv)\;\mathrm{at}\;z=0\right)
&=-\sqrt{\frac{f}{\,2\nu\,}}(1+i)B\\
&=-\sqrt{\frac{f}{\,2\nu\,}}(1+i)(u_0+iv_0-u_g)\\
&=-\sqrt{\frac{f}{\,2\nu\,}}\{u_0-v_0-u_g+(u_0+v_0-u_g)i\}
\end{aligned}
$$

を得ます。

　これが$${u_0+iv_0}$$の実数倍であるための必要十分条件は、$${\alpha_0}$$を$${x}$$軸の正の向きから測った、風速$${(u_0,v_0)}$$の方位角（即ち、$${\alpha_0=\arg{(u_0+iv_0)}}$$）として、

$$
\begin{aligned}
&\det
\begin{pmatrix}
u_0-v_0-u_g&u_0+v_0-u_g\\
u_0&v_0
\end{pmatrix}=0\\
\Longleftrightarrow&
u_0^2+v_0^2=u_g(u_0-v_0)\\
\Longleftrightarrow&
\sqrt{u_0^2+v_0^2}=u_g(\cos{\alpha_0}-\sin{\alpha_0})
\end{aligned}
$$

です（最後の式を用いると、接地境界層上端における風向、風速から地衡風速を求められます）。

　最後の式より、

$$
\begin{aligned}
u_0+iv_0-u_g
&=u_g\{(\cos{\alpha_0}-\sin{\alpha_0})e^{i\alpha_0}-1\}\\
&=u_g\{-\sin{\alpha_0}(\cos{\alpha_0}+\sin{\alpha_0})+i\sin{\alpha_0}(\cos{\alpha_0}-\sin{\alpha_0})\}\\
&=\sqrt2 u_g\sin{\alpha_0}\left\{\cos\left(\alpha_0+\dfrac{3}{\,4\,}\pi\right)+i\sin\left(\alpha_0+\dfrac{3}{\,4\,}\pi\right)\right\}\\
&=\sqrt2 u_g\sin{\alpha_0}\exp\left({\left(\alpha_0+\dfrac{3}{\,4\,}\pi\right)i}\right)
\end{aligned}
$$

であり、求める風速は、

$$
\begin{aligned}
u-u_g+iv&=(u_0-u_g+iv_0)\exp{\left(-\sqrt{\dfrac{f}{\,2\nu\,}}(1+i)z\right)}\\
&=\sqrt2 u_g\exp{\left(-\sqrt{\dfrac{f}{\,2\nu\,}}z\right)}\sin{\alpha_0}\exp\left({\left(\dfrac{3}{\,4\,}\pi+\alpha_0-\sqrt{\dfrac{f}{\,2\nu\,}}z\right)i}\right)
\end{aligned}
$$

の両辺の実部、虚部を比較して得られる

$$
\begin{cases}
u=u_g+\sqrt2 u_g\exp{\left(-\sqrt{\dfrac{f}{\,2\nu\,}}z\right)}\sin{\alpha_0}\cos{\left(\dfrac{3}{\,4\,}\pi+\alpha_0-\sqrt{\dfrac{f}{\,2\nu\,}}z\right)}\\
v=\sqrt2 u_g\exp{\left(-\sqrt{\dfrac{f}{\,2\nu\,}}z\right)}\sin{\alpha_0}\sin{\left(\dfrac{3}{\,4\,}\pi+\alpha_0-\sqrt{\dfrac{f}{\,2\nu\,}}z\right)}
\end{cases}
$$

で与えられます。

　各$${z}$$に対する風速ベクトルの終点を$${uv}$$平面上にとると、地衡風$${(u_g,0)}$$に収束する螺旋を描きます。これをエクマンスパイラルといいます。エクマン層内では、本記事で課した条件の下、東寄りの風が高度上昇に伴い、南北成分の符号を交互に変えながら、地衡風に漸近していくことが分かります。