梅雨前線

梅雨前線#

2018年の西日本豪雨に関する力学解析[9]を行います。これまでは取得した変数をそのまま描画していましたが，MetPyの助けを借りて別の物理量を計算してみましょう。 MetPyには気象のデータ解析に関する様々な機能がありますが，ここでは単位metpy.units.unitsと計算metpy.calcを利用します。

import numpy as np
import xarray as xr
from metpy.units import units
import metpy.calc as mpcalc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

相当温位#

理想気体の状態方程式は

(1)#\[ p\alpha = RT \]

と表すことができます。ここで\(p\,\mathrm{Pa}\)は気圧，\(\alpha=1/\rho\,\mathrm{m^3/kg}\)は密度の\(\rho\,\mathrm{kg/m^3}\)の逆数で比容，\(T\) Kは気温です。定数\(R\)は

\[R=\frac{R^*}{M}\]

で\(R^*=8314.3\,\mathrm{JK^{-1}kmol^{-1}}\)，\(M\)は\(\mathrm{kg}\)で表した分子量です。

乾燥空気の平均分子量は\(M_\mathrm{d}=28.97\)であり，乾燥空気に対する気体定数は\( R_\mathrm{d}=287\,\mathrm{JK^{-1}kg^{-1}}\)です。

気体の状態方程式(1)を用いると熱力学第一法則は

(2)#\[ \mathrm{d}q = c_p\mathrm{d}T - \alpha\mathrm{d}p \]

と表すことができます。ここで\(\mathrm{d}q\)は系に与えられた差分加熱，\(\mathrm{d}T,\,\mathrm{d}p\)はそれぞれ温度と気圧の変化を表します。

\[ c_p = \left(\frac{\mathrm{d}q}{\mathrm{d}T}\right)_{p} \]

は定圧比熱です。乾燥空気の定圧比熱は\(c_p = 1004\,\mathrm{JK^{-1}kg^{-1}}\)です。

断熱\(\mathrm{d}q=0\)のとき，参照気圧\(p_0\)での気温を\(\theta\)，気圧\(p\)での気温を\(T\)として (2)を\(p_0\)から\(p\)まで積分すると

\[ \theta = T\left(\frac{p_0}{p}\right)^{R/c_p} \]

が得られます。 \(\theta\)は温位と呼ばれ，気圧\(p\)における温度\(T\)の気塊を断熱過程で\(p_0\)まで移動したときの温度を表します。

凝結を伴う場合は凝結熱が生じます。水に対する飽和混合比が\(w_\mathrm{s}\)であるとき，凝結に伴い乾燥空気に放出される熱は \(\mathrm{d}q = -L\mathrm{d}w_\mathrm{s}\)です。凝結した液体の水が系の外に出るものとすると，不可逆で非断熱となります。これを偽断熱過程と呼びます。

混合比の変化に対する温度変化が相対的に小さい（\(\mathrm{d}T/T\ll \mathrm{d}w_\mathrm{s}/w_\mathrm{s}\)）ことを仮定すると

(3)#\[ \frac{\mathrm{d}q}{T} = c_p\frac{\mathrm{d}\theta}{\theta}=-\frac{L}{T}\mathrm{d}w_\mathrm{s}\approx-\mathrm{d}\left(\frac{Lw_\mathrm{s}}{T}\right) \]

(3)を積分し，積分定数を定めるため\(w_\mathrm{s}/T=0\)となる温位を\(\theta=\theta_\mathrm{e}\)と表すと

(4)#\[ \theta_\mathrm{e} = \theta\exp\left(\frac{Lw_\mathrm{s}}{c_pT}\right) \]

が得られます。混合比\(w_\mathrm{s}=0\)となるときの温位\(\theta_\mathrm{e}\)を相当温位といいます。気塊を全ての水蒸気が凝結するまで偽断熱的に膨張させた後，参照気圧まで断熱的に圧縮すると\(\theta_\mathrm{e}\)が得られます。このとき凝結熱が発生し，凝結した水は落下して気塊から離れるものとします。

相当温位は気温や湿度が高いほど，大きくなります。相当温位が高いことは，暖かくて湿っていることを意味しています。

MetPyを使って相当温位を計算してみましょう。まずmpcalc.dewpoint_from_relative_humidity()を使って，気温と相対湿度から露点温度を計算します。次にmpcalc.equivalent_potential_temperature()を使って，気温と相対湿度から相当温位を計算します。 MetPyでは，相当温位の計算に持ち上げ凝結高度を計算する精度の高い式[10]を用いています。

なお，MetPyの関数の引数には単位が必要です。 levは850という数字なのでunits.hPaを掛けて単位付きにしています。

dewpoint = mpcalc.dewpoint_from_relative_humidity(ds.air, ds.rhum)
thetae = mpcalc.equivalent_potential_temperature(lev*units.hPa, ds.air, dewpoint)

fig = plt.figure(figsize=[8,6])
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
p = ax.contourf(lon, lat, thetae, transform=ccrs.PlateCarree(),
               cmap="coolwarm", levels=np.arange(296, 362, 2))
ax.quiver(lon, lat, ds.uwnd, ds.vwnd)
ax.coastlines()
fig.colorbar(p)
plt.show()

_images/d8299926092b885413f599b2e16ddb9017393670eaceedc0fe2d483a4238c58e.png

大陸から高相当温位の空気が東シナ海を渡って日本列島に伸びつつある様子が分かります。

Qベクトル#

中高緯度では，気圧傾度力とコリオリ力がおよそ釣り合っています。

(5)#\[ u_\mathrm{g} = -\frac{1}{f_0}\frac{\partial\phi}{\partial y},\; v_\mathrm{g} = \frac{1}{f_0}\frac{\partial\phi}{\partial x} \]

これを地衡風平衡といいます。\(f_0\)は注目している緯度におけるコリオリパラメータです。風は地衡風成分と非地衡風成分の和として表すことができます。

\[ u = u_\mathrm{g} + u_\mathrm{a},\; v = v_\mathrm{g} + v_\mathrm{a} \]

中高緯度でよく成り立つ準地衡風近似では，水平風はおよそ地衡風で近似されますが，水平発散を考慮し非地衡風成分により鉛直流\(\omega\)が生じます。

\[ \frac{\partial u_\mathrm{a}}{\partial x} + \frac{\partial v_\mathrm{a}}{\partial y} + \frac{\partial \omega}{\partial p}=0 \]

(5)を気圧\(p\)で微分して静力学平衡

\[ \frac{\partial\phi}{\partial p}=-\alpha=-\frac{RT}{p} \]

を仮定すると温度風の関係

(6)#\[ f_0\frac{\partial u_\mathrm{g}}{\partial p} = \frac{R_\mathrm{d}}{p}\frac{\partial T}{\partial y},\; f_0\frac{\partial v_\mathrm{g}}{\partial p} = -\frac{R_\mathrm{d}}{p}\frac{\partial T}{\partial x},\; \]

が得られます。

準地衡風近似の下で熱力学の式は

(7)#\[ \frac{\mathrm{d_g}T}{\mathrm{d}t} - \left(\frac{p}{R_\mathrm{d}}\right)S_0\omega = 0 \]

と表されます。

ここでラグランジュ微分は

(8)#\[ \frac{\mathrm{d_g}}{\mathrm{d}t} = \frac{\partial}{\partial t} + u_\mathrm{g}\frac{\partial}{\partial x} + v_\mathrm{g}\frac{\partial}{\partial y} \]

地衡風による移流で表されます。 \(S_0=-\alpha_0\mathrm{d}\ln\theta_0/\mathrm{d}p\)は安定度を表します。 \(\alpha_0=\alpha_0(p),\,\theta_0=\theta_0(p)\)は，それぞれ基本場の水平一様な比容と温位を表します。

温度風(6)の両辺のラグランジュ微分(8)を取り，(7)を用いて整理すると鉛直流\(\omega\)の診断式

(9)#\[ S_0\nabla_\mathrm{h}^2\omega + f_0\frac{\partial^2\omega}{\partial p^2} = -2\nabla_\mathrm{h}^2\cdot\mathbf{Q} \]

が得られます。[11]

ここで\(\nabla_\mathrm{h}^2\)は水平ラプラシアンです。ここで

\[ \mathbf{Q} \equiv -\frac{R_\mathrm{d}}{p}\left(\frac{\partial \mathbf{v}_g}{\partial x}\cdot\nabla_\mathrm{h}T,\; \frac{\partial \mathbf{v}_g}{\partial y}\cdot\nabla_\mathrm{h}T\right) \]

をQベクトルといいます。 Qベクトルの収束は上昇流，発散は下降流を表します。

MetPyを使ってQベクトルを計算しましょう。 Qベクトルでは準地衡風近似が仮定されているので，mpcalc.geostrophic_wind()を使ってジオポテンシャル高度から計算した地衡風を用います。次に地衡風と気温からmpcalc.q_vector()を使って，Qベクトルを算出します。 Qベクトルからその発散をmplcalc.divergence()を使って計算しています。

ug, vg = mpcalc.geostrophic_wind(ds.hgt)
qx, qy = mpcalc.q_vector(ug, vg, ds.air, lev*units.hPa)
dq = mpcalc.divergence(qx, qy)

Qベクトルの発散を陰影で，気温を等値線で，Qベクトルを矢印で描画します。

fig = plt.figure(figsize=[8,6])
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
p = ax.contourf(lon, lat, dq, transform=ccrs.PlateCarree())
c = ax.contour(lon, lat, ds.air, levels=np.arange(276, 296, 2), cmap="coolwarm")
ax.clabel(c, inline=True)
ax.quiver(lon, lat, qx, qy)
ax.coastlines()
fig.colorbar(p)
plt.show()

_images/074d4bede9cd4ae4537aaefe89cefadf0aef2e6f1f71147f0715315868e94c6d.png

Qベクトルは北風の中心から発散しています。日本海にある発散は西日本豪雨に先行して接近した台風第7号の北風に対応しています。

南向きのQベクトルは東西に伸びる等温線を横切って，鉛直流を誘導することにより梅雨前線を強化しています。

梅雨前線

Contents

梅雨前線#

データの取得#

相当温位#

Qベクトル#