macroscope

( はてなダイアリーから移動しました)

大気のエネルギー収支解析のための質量収支の補正について

学術覚え書き

【まだ書きかえます。どこをいつ書きかえたかを必ずしも明示しません。】

- 1 -
2017年5月27日、日本気象学会の大会のなかで、大気によるエネルギー輸送量に関するデータ解析の研究発表があった。わたしがながねん(1980年代-2000年代)かかわってきた研究方法であり(成果を論文にしたものが少ないことがなさけないが)、同じ技術的問題点にぶつかっていることに気づいたので、質問してみた。発表者はその問題に気づいて対策をとろうとしていることがわかったので、ひと安心した。他方、多くの研究者にデータを提供する活動に向けて要望を述べておきたいと思い、その発言もした。

わたしのエネルギー収支解析の結果は、気象学の概論に出てくるような情報であり(実際、ある入門書[読書メモ]に図を提供した)、それがどのように導かれたかの基本は、気象学の学部上級・大学院初級くらいの教科書を見ればわかるだろう。しかし、実際に計算しようとすると、データの制約からくる問題点への対処が必要になる。その問題点と対処方法については、研究論文には書くものの、論文全体に対して小さい部分だから、読んだ人が再現できるほど詳しくはない[注]。多くの場合は、同じ「研究室」などで顔を合わせるうちに先輩から後輩に伝えられていくような、いわば専門的暗黙知になっている。

  • [注 (2017-06-03 追加)] 増田(1988)の日本語の解説(博士論文にも含めた)では、2つのデータセットについての、わたしが出会った問題点を、その雑誌の解説記事のページ数の許す範囲でなるべく詳しく述べた。その重要な部分は質量収支に関するものである。ただし、対処方法について述べたものではない。

わたしはこの数年、大学の非常勤講師として概論的な講義をしてはいるが、研究の指導をする立場になく、研究所でも同じ研究チームにいた同僚が減ってほぼ単独で仕事をしていたから(人を使って研究するリーダーとしての能力をもてなかったのがいけないのだろうが)、データ解析上の問題点と対処方法について、話題にする機会があまりなくなっていた。今になって、わたしが得た知識を記録しておく義務があると感じている。その内容はまだ他人に提供するのに適した形に整理できていないが、ひとまず「たたき台」として書き出しておく。

- 2 -
エネルギー保存の法則に基づく、大気のエネルギー収支解析について、わたしが知っていることは、このブログの「大気のエネルギー、エネルギー収支、熱収支」のシリーズの記事[第1部(2016-12-09)]、[第2部(2016-12-10)]、[第3部(2016-12-11)]、[第4部(2016-12-13)]、[第5部(2016-12-26)]に書いた。

保存則を前提とした枠組みで、観測データや、気象・気候モデルによるシミュレーションの結果や、観測と予測型シミュレーションとを組み合わせたデータ同化プロダクト(古い用語では「客観解析データ」)を材料として、エネルギー収支の計算がおこなわれる。

エネルギー収支の定式化では、対象が質量保存則を満たしていることを前提とすることがある。現実の世界では(適切に定式化すればだが)質量保存が成り立っているはずだ。しかし、データから計算すると、質量保存が成り立っていればゼロになるべきところに、無視できない残差が出ることが多い。そして(おおざっぱな言いかたになるが)、エネルギーの流れの定式化が、空気の質量の流れに単位質量あたりのエネルギーをかけるような形のものだと、エネルギー収支には、質量収支の残差に単位質量あたりのエネルギーをかけたような量のオーダーの残差が出てしまう。使いものになる結果を得るためには、一方で、質量保存則からのずれがなるべく生じないように解析手順をくふうすること、他方で、生じてしまったずれを減らすように物理量(風速など)の値の修正を追加すること、が必要になる。

- 3 -
わたしは、修士論文とそれをもとにしたMasuda (1983)の学術雑誌論文では、「客観解析データ」を使って3次元の「大気の熱収支」の解析をした(2節にあげたリンク先の第3部参照)。そこで鉛直座標は気圧 p をとっている。鉛直速度にあたるものは ω = dp/dt で、鉛直速度 w、空気の密度ρ、重力加速度 g とすれば、ω = -ρ g w である。

多くの場合、風速の水平2次元成分の情報はある(風向と風速の大きさの情報があれば変換できる)が、鉛直成分の観測値はない。しかし、質量保存の式は p 座標で書けば3次元の発散が0という形になるから、水平発散を計算して鉛直に「積み上げる」ことによってωを得ることができる。しかし、現実には風のデータが完全でないので、下端(地表面)と上端(理屈のうえでは「大気上端」)の一方の境界条件としてもっともなものを与えて積み上げると、他方の境界での状態が物理的・気象学的常識から見て変なものになる。

わたしの修士論文のときは、材料が等圧面高度と気温のデータで水平風速もなかった。対象が中高緯度であることを前提として、準地衡風近似を使い、水平風速としては地衡風を使い、水平発散は準地衡風渦度方程式によって求め、それを「積み上げて」ωを求めたのだった。

ωについての地表面での境界条件としては、水平風が地形にあたって押し上げ・押し下げられる効果による上昇・下降流があるとした。このとき、実際に地形にぶつかる境界層内の風は、自由大気(境界層の外)を想定した客観解析データの風速よりも弱いだろう、と考えて補正係数をかけることもある。また、地表面の高さとして何をとるべきかも自明ではない(あとであらためて論じる)が、ひとまず現実の地面の標高をあらわす緯度経度格子データセットをもってきて自分の解析作業で使う格子点に内挿したものを使うことにした。

データ解析対象領域の上端での境界条件をどう考えるかはむずかしい。単純に考えれば「大気上端」では p = 0 で ρ = 0 だから、w があっても ω = 0 としてよいだろうと考えられる(わたしの修士論文ではそのように考えた)。ただし、たとえば対流圏全体のデータ解析ならば、実際の解析対象領域の上端は、対流圏界面の少し上、つまり成層圏下部にある。そこでは平均の鉛直流が無視できず、その大きさを見積もる必要があるかもしれない。(わたしの修士論文では、記載した結果は対流圏の下半分に限ったのだが、計算の枠組みは対流圏全体を扱おうとしたものだった。)

データから風速の水平発散を求め、pについて地表面気圧psから0まで鉛直積分(の近似となる差分計算)をすると、p = 0 でのωは0になるべきところだが0でない残差が出る。この残差を鉛直差分の各層に分配する。分配方法としては、p座標の間隔あたり均一に分配する(ωの補正量は p - psに比例する、つまりpの1次関数になる)のが代表的だ。わたしの修士論文では、地表付近に比べて上層ほど材料データの誤差が大きいだろうと考えて、pの小さいところほど重みが大きくなるpの1次関数で分配し、ωの補正量がpの2次関数になるようにした。これと均一に分配するのとどちらがよいかはいちがいには言えないと思う。なお、補正式をつくる際には、鉛直座標の関数を、p よりも実際の高さに近似的に比例する log p の関数とみたほうがよいという考えもあるが、このときは採用しなかった。

質量保存を考えて鉛直流(風速の鉛直成分)を補正したのならば、それとつじつまのあうように風速の水平成分も補正するべきだ。実際、1980年ごろにも、中小規模の気象の解析では、3次元の風速の場を質量保存を満たすように補正する「佐々木(嘉和)の変分法」を使う研究もあった。しかし、それは当時としては大量の計算機資源を必要とした。大規模の気象では、風の水平成分の内わけは回転成分が主で発散成分はわずかだから、補正の必要はないと考えられることが多かった。(今の計算機事情と観測誤差を前提とすれば補正をする意義があるかどうかは、わたしはよく知らない。)

- 4 -
Masuda (1988)の論文(わたしの博士論文の主要な部分でもある)になった研究では、2節であげたリンク先の第2部で紹介したPeixoto & Oortの本の式と同様な定式化にしたがって、鉛直積算したエネルギーの水平輸送を計算した。ただしOortが観測データだけに基づいたのとちがって、数値予報とともにつくられた「客観解析データ」から計算した。

この研究では、鉛直に層を分けたエネルギー収支を考えることはあきらめた。そうすることによって、ωを求めるというややこしい問題を回避したのだった。

データが与えられているのは水平方向には緯度経度の格子、鉛直方向には p 座標の標準気圧面(1000 hPa, 850 hPa, ...)である。鉛直積算の下端(pの値が大きいほうの端)は地表面で、そこでのpの値は時刻ごとに変わるのだが、地表面気圧の値はデータに含まれていなかったので、地表面の高さと、標準気圧面の等圧面高度や温度などから、内挿・外挿で計算する必要があった。

そこで、地表面の高さとしてどんなデータを使うべきかという問題がある。「客観解析データ」は観測値と予報値(予報モデルの出力)とを重みつきで混合したものだが、重みは明示されていなかった。観測値ならば現実の地表面の標高、予報値ならば予報モデルの地表面の標高を使うのがよさそうだと思われた。予報モデルの地表面は現実の地表面の値をもとにしているが、差分計算のノイズを小さくするために平滑化してあったり、風に対する山による抵抗を適度な大きさにするために実際の標高の升目平均よりも高めの値になっていたりすることもある。わたしは、質量収支・エネルギー収支を計算するうえではモデルの地表面の標高を使うべきだと判断した。地表面の標高は時間とともに変化しない量だが、当時の「客観解析データ」には含まれていなかったので、「客観解析データ」をつくる事業にたずさわっていたかたにお願いして提供していただいた。

地表面の境界条件を適切に与えると、雑に与えた場合に比べて、質量保存の残差はだいぶ小さくなった。しかし、エネルギー収支計算のためには、まだ補正が必要だった。

詳しいことを言えば、大気の質量の総量も時間とともに変化する。水蒸気以外の成分の質量は、1日から数年の時間規模では、変化しないとみなしてよさそうだが、水蒸気の質量が、降水と地表面からの蒸発の差に伴って変化することは、大気の質量の季節変化にも現われている(教材ページ[大気の質量の季節変化]参照)。ただし1980年代の段階ではこのこと自体が先端の研究課題であり、質量保存の補正の際に考慮することは現実的でなかった。

次に、鉛直積算した質量保存の式は、鉛直積算した質量の水平輸送の水平発散と、鉛直積算した質量(地表面気圧に比例する)の局所的時間変化とによっている。モデル出力だとすれば、水平発散と組み合わせるべきものはモデルの時間差分の1ステップの局所的時間変化だと考えられるけれども、その値は提供されていない。提供されているデータの時間間隔は12時間(6時間のこともあったが期間限定)だった。水平発散を計算した時刻の前後それぞれ12時間(合計24時間)の平均の局所的時間変化ならば計算できるけれども、それを考慮に入れることによって計算手順が複雑になるわりには補正が改善されないだろうと思った(ただしその点をきちんと確認しないままで来てしまった)。それで、わたしは、鉛直積算した質量の水平輸送の水平発散を消すという方針で補正を考えた。

1988年に出した論文では、南北方向のエネルギー輸送だけを考えたので、各緯度円を通過する(見かけ上の)質量の重みつきの東西鉛直平均の風速南北成分を求め、それをその緯度円にある東西・鉛直断面の各格子点の風速南北成分からひく、という方法をとった。

それと並行して、水平2次元の(鉛直積算した)エネルギー輸送の計算もしていた。そちらでは、鉛直積算した質量輸送の水平発散が消えるように補正することにした。
2次元の流速(風速)の場からその発散成分を求めるには、まず水平発散を計算し、Poisson方程式をといて(空間座標による2回積分にあたる)、velocity potentialという量を得る。そしてvelocity potentialの水平gradient (1階微分)によって、速度の発散成分を得る。
補正計算の場合は、この「速度」の代わりに鉛直積算した質量輸送ベクトルを入れて、その発散成分を得る。それを鉛直積算した質量(=地表面気圧/重力加速度)で割って速度の次元の量にし、「客観解析」の風速からこれをひく。

球面上のPoisson方程式をとくには、アメリカのNational Center for Atmospheric Research (NCAR)のScientific Computing Divisionの数値解析専門家 SwarztrauberさんたちによるFISHPAKという、差分近似してできた連立1次方程式の解を直接求める方法によるFortranサブルーチン集があったので(先輩の中村一さんがNCARから持ち帰っていたので)それを使うことができた。その後、球面調和関数展開し、展開結果の成分ごとに適切な係数をかけて、また格子点値にもどす、という方法も使った。球面調和関数展開にもSwarztrauberさんたちによるSpherepackというサブルーチン集を使うことができた。

- 5 -
p座標では、地表面の鉛直座標値が一定しないが、気象モデルではたいてい、各格子点の地表面の鉛直座標値が時間とともに変化しないような鉛直座標を使って、数値を与える鉛直位置(レベル)を与えることが多い。そのような鉛直座標としては、まず気圧pを地表面気圧psで割ったσ (シグマ)座標がある。地表面ではσ = 1 である。σ座標は山に伴うでこぼこが大気上層にまで現われるという欠点がある。地表付近ではσに近く、高いところでは p に比例するような混成座標が使われることがふえてきた。

データ同化の結果は、気象モデルの格子で得られるが、多くの利用者が標準気圧面での数値がほしいので、公式プロダクトとしては、モデルの鉛直座標からp座標への内挿をした結果がおもに提供されてきた。しかし、モデル鉛直座標のまま内挿しないで提供してほしいという利用者の声も知られて (わたしのように収支解析をしたい人のほかに、3次元流跡線解析をしたい人、大気・陸面相互作用の計算をしたい人などからも希望があったと思う)、NCEP/NCAR Reanalysis、NCEP/DoE Reanalysis 2、JRA25、JRA55などでは、モデルの格子のままのデータも公式プロダクトに含まれるようになってきた。モデルの鉛直座標面で与えられたデータからそのまま計算すると、質量収支の残差は、なくなるわけではないが、等圧面に内挿されたデータで計算する場合よりもだいぶ小さくなる。

【[2017-05-28補足] なお、エネルギー収支の計算で、質量収支の誤差の影響を小さくするために、いくらかくふうできることがある。位置エネルギーの輸送を計算する際には、ジオポテンシャル高度 z から、対象となる大気全体の(質量で重みづけした) zの平均値に近い定数 z0 (たとえば 6000 m)をひいたものが運ばれるように考えて計算する。同様に、内部エネルギー(の水蒸気によらない部分)の輸送や圧力による仕事の項を計算する際には、気温 T から同様な T の平均値に近い定数 T0 (たとえば 255 K)をひいたものが運ばれるように考えて計算する。これで、質量収支の誤差がエネルギー収支に一方的にきくのを避けることができる。】

近ごろは、CMIPプロジェクトなどの気候モデル出力のデータが公開提供されるようになってきた。その大気部分の3次元データは、鉛直座標として気圧をとったものが主とされている。モデルの鉛直座標のままのデータも提供されることを希望する。ただし、データ総量やファイル数がふえてしまうので、データ提供サービスにのせることがむずかしいのかもしれない。もしそうだとすると、鉛直内挿で誤差が生じるのは残念だが、等圧面のデータを使うしかないのかもしれない。

また、モデルの鉛直座標のデータを使うことには欠点もある。大気によるエネルギー輸送は、東西全経度平均した量で表現できる「子午面循環」による輸送と、その平均からのずれどうしの積による「渦(eddy)による輸送」(「波による輸送」ともいう)とに分けて考えることができる。(さらに、時間平均した量とそれからのずれとに分けることも組み合わせることが多い。[教材ページ「「渦(eddy)輸送」: 「積の平均」と「平均の積」の違い」]参照。) このとき、東西平均の対象としては、ほぼ同じ高さのところの値を使う必要がある。p座標ならばだいたい同じ高さなのでよいのだが、モデルの鉛直座標を使うと、鉛直座標値が同じでも山の上は平地の上よりもだいぶ高いところを見ることになって、実質的に子午面循環と渦輸送とをうまく取り出せないのだ。このように分けたいときは、気圧面で計算するのが順当だろう。

文献