はじめに
季節〜気候スケールの大気の水平の流れ場を表現する場合には、流線関数(stream function)や速度ポテンシャル(velocity potential)を用いることがある。
これらの量を東西・南北風速から計算するには、後述するようにPoisson方程式を解く必要があり、localな差分計算のみからは求まらない。このため、例えばMetPyのmetpy.calcには関数が用意されておらず、別のパッケージを用いる必要がある。
ここでは、xarrayで読み込んだデータに対して、windspharmというPythonパッケージを用いて計算する方法を説明する。windspharmは、全球の風速データに対して球面調和関数展開を適用し、流線関数・速度ポテンシャルなどの風に関わる諸量を計算することができる。
PythonのwindspharmパッケージではなくCDOを使う方法については、以下の記事を参照されたい。
なお、詳しくない読者のために流線関数や速度ポテンシャルについての数学的基礎や気象学的応用について記したが、必要のない読者は読み飛ばしてもらって構わない。
数学的基礎
まず、流線関数・速度ポテンシャルに関する数学的な基礎を説明する。気象データを扱う場合は球面上での計算になるが、簡単のためにデカルト座標での計算を考える。
Helmholtz分解定理
水平風ベクトル ${\bf V}$ は、東西風 $u$ と南北風 $v$ を用いて ${\bf V} \equiv (u, v, 0)^\mathrm{T}$ と成分表示される。
気象学でよく使われる2次元のHelmholtz分解定理は、任意のベクトル場（ここでは風の場）を発散なし成分 $\mathbf{V_ψ}$ と渦なし成分 $\mathbf{V_χ}$ に一意に分解されることを示す。
\displaylines{
\mathbf{V} = \mathbf{V}_ψ + \mathbf{V}_χ \\
∇ \cdot \mathbf{V}_ψ = 0 \\
∇ \times \mathbf{V}_χ = \mathbf{0}
}
流線関数・速度ポテンシャル
水平発散なしの流れ場に対しては、一般に以下のような流線関数が導入できる。
\mathbf{V}_ψ = (u_ψ, v_ψ, 0)^\mathrm{T} = \left(-\frac{∂ψ}{∂y}, \frac{∂ψ}{∂x}, 0 \right)^\mathrm{T}
流線関数の等値線は、その名の通り流れのベクトル $\mathbf{V_ψ}$ に平行な流線になる。
水平渦なしの流れ場に対しては、一般に以下のような速度ポテンシャルが導入できる。
\mathbf{V}_χ = (u_χ, v_χ, 0)^\mathrm{T} = \left(\frac{∂χ}{∂x}, \frac{∂χ}{∂y}, 0 \right)^\mathrm{T}
速度ポテンシャルは、流れのベクトル $\mathbf{V_χ}$ に対して直交する。
なお、流線関数と速度ポテンシャルいずれも定数項の任意性があるが、通常は領域平均が0になるように定数項を決めることが多いと思われる。
以上より、水平風は（発散なし成分の）流線関数 $ψ$ と（渦なし成分の）速度ポテンシャル $χ$ の2つのスカラー量で表現できることが分かる。
鉛直渦度・水平発散との関係
鉛直渦度 $ζ$ は流線関数 $ψ$ のLaplacianである。
ζ
\equiv \mathbf{k} \cdot ∇ \times \mathbf{V}
= \mathbf{k} \cdot ∇ \times \mathbf{V}_ψ
= \frac{∂^2 ψ}{∂x^2} + \frac{∂^2 ψ}{∂y^2}
= \nabla^2 ψ
ここで $\mathbf{k}$ は鉛直($z$)方向の単位ベクトルである。
水平発散 $δ$ は速度ポテンシャル $χ$ のLaplacianである。
δ
\equiv ∇ \cdot \mathbf{V}
= ∇ \cdot \mathbf{V}_χ
= \frac{∂^2 χ}{∂x^2} + \frac{∂^2 χ}{∂y^2}
= \nabla^2 χ
このことから、水平風 $(u,v)$ から、流線関数・速度ポテンシャル $(ψ,χ)$ を求めるためには、まず渦度・発散 $(ζ,δ)$ を求めた後、2つのPoisson方程式
\displaylines{
\nabla^2 ψ = ζ \\
\nabla^2 χ = δ
}
を解く必要があることが分かる。
気象学的応用
流線関数や速度ポテンシャルが季節予報で好まれる理由
流線関数と鉛直渦度はどちらも渦を表す量であるが、鉛直渦度のほうが水平スケールの細かい構造を表現する。というのも微分演算子Laplacianはハイパスフィルタとして機能するからである。例えば、$ψ = \sin(kx)$ と東西波数 $k$ の正弦波で表現される場合、$ζ = -k^2 \sin(kx)$ となるため、より高波数(短波長)成分が強調されるためである。
このような性質から、短期予報では細かい擾乱を把握するために鉛直渦度や水平発散に着目することが多いのに対して、季節予報ではより大きな循環場を把握するために流線関数や速度ポテンシャルを用いることが多い。
速度ポテンシャルと対流活発域
速度ポテンシャルは、熱帯の対流活動等に伴う大規模な発散・収束等を監視するのに利用される。熱帯域では通常下層850hPaで収束(発散)、上層200hPaで発散(収束)、その間の高度で上昇流(下降流)というパターンになることが多いため、200hPaの速度ポテンシャルが小さい発散域では対流活動が活発で上昇流域になっていると解釈される。
なお、前述のように水平発散と速度ポテンシャルは表現する空間スケールが異なり、水平発散は個々の対流(あるいは降水)のパターンに対応するのに対し、速度ポテンシャルはMJOやENSOのようなある程度の広さを持った対流活発域と対応する。
流線関数と水平風の流線
流線関数 $ψ$ は発散なし成分 $\mathbf{V_ψ}$ に対する流線関数であるが、気象学ではあたかも流線関数の等値線が元の水平風 $\mathbf{V}$ の流線であるかのように解釈される。つまり値が大きいほうを右手に向いて風が吹き、流線関数の等値線の間隔が狭いほど風が強いと解釈される。
この背景には、発散なし成分と渦なし成分の大きさについて
|\mathbf{V}_χ| \ll |\mathbf{V}_ψ|
という関係が経験的に成り立つことに基づいている。この関係が成り立つ理由については、筆者の推測だが、水平収束は連続の式より（浮力が必要な）上昇流を伴うのに対し、渦を作る場合は上昇流を伴う必要がないことが関係していると思われる。
流線関数を使うメリットの一つは、熱帯域も含めて等値線を描くだけで風速・風向を把握できる点にある。中高緯度では地衡風平衡が概ね成り立つため、等ジオポテンシャル高度線から風を把握できる。しかし、熱帯域では地衡風平衡が成り立たないため、等ジオポテンシャル高度線から風を推測することができない（そもそもジオポテンシャル高度勾配はほとんど存在しない）。流線関数はジオポテンシャル高度と比べて、熱帯域の風の把握に使えるという利点がある。
Pythonでの計算
windspharmのインストールはcondaを用いて行うことができる。このあとxarrayを使うため、インストールしていない場合はあらかじめインストールしておく。
conda install -c conda-forge windspharm
(追記) Intel Mac の場合はこれでインストールできるが、Apple Siliconの場合はインストールすることができない。windspharmの最終更新は2018年であり、メンテナンスがなされていないようである。
windspharmは入力に対する様々なインターフェースが用意されている。
本パッケージは、numpy と pyspharm が利用可能であることを最低限必要とし、インストールには setuptools が必要です。windspharm.iris インターフェースは、iris パッケージが利用可能な場合にのみ使用できます (iris のドキュメントを参照)。windspharm.cdms インターフェースは、cdms2 モジュールが利用可能な場合にのみ使用できます。このモジュールは、UVCDATプロジェクトの一部として配布されています。windspharm.xarray インターフェースは、xarray パッケージが利用可能な場合にのみ使用できます（xarray のドキュメントを参照）。
ここではxarrayを用いた方法を説明する。xarrayインターフェースについての公式ドキュメントは以下のページにある。
次に示すサンプルコードは、風速の適当なGRIBファイル（６か月アンサンブル数値予報モデル統計GPV（全球域）３か月統計値）から、流線関数と速度ポテンシャルを計算する。
import matplotlib.pyplot as plt
import xarray as xr
from windspharm.xarray import VectorWind
# ファイル読み込み
file_u = "grib/Z__C_RJTD_20211122000000_EPSC_GPV_Rgl_Gll2p5deg_Lp200_Pwu_E3em_grib2.bin"
file_v = "grib/Z__C_RJTD_20211122000000_EPSC_GPV_Rgl_Gll2p5deg_Lp200_Pwv_E3em_grib2.bin"
ds_u = xr.open_dataset(file_u, engine="cfgrib")
ds_v = xr.open_dataset(file_v, engine="cfgrib")
# 時刻0の2次元配列を変数に代入
u = ds_u["u"][0]
v = ds_v["v"][0]
# 流線関数と速度ポテンシャルを計算
# 流線関数のみの場合は w.streamfunction()
# 速度ポテンシャルのみの場合は w.velocitypotential()
w = VectorWind(u, v)
sf, vp = w.sfvp()
# 描画
sf.plot.contour(levels=20, colors="black")
# vp.plot.contour(levels=20, colors="black")
plt.show()
plt.clf()
plt.close()
実行すると、以下のように流線関数や速度ポテンシャルが描画される。
流線関数
速度ポテンシャル
