(大学院) 熱収支水収支特論

1学期 月曜2限 4月14日から

このウェブページは http://macroscope.world.coocan.jp/ja/edu/ris/energy_balance/ です。

2025年度の授業予定

• 地表面エネルギー収支の概念を学ぶ (松山・増田 (2021) 第5章)
• 世界の地表面エネルギー収支・水収支各項の地理的分布や季節変化をみる。例、JRA3Q または JRA55 再解析
• 地表面エネルギー収支・水収支の基本と、エネルギー収支項を地上気象観測値からもとめる方法を学ぶ。 近藤 (1994) 『水環境の気象学』 または 近藤 (2000) 『地表面に近い大気の科学』 のその関連の部分を読む。
• 地上気象観測値から地表面エネルギー収支項を計算してみる。例、熊谷キャンパスでの観測値による。

おおまかには2024年度と同様です。 くわしい内容や日程などは、受講生との相談で変えていきます。

• 第1回 4月14日(月) 2限 (10:55から)、1号館 1306
  • この授業の予定
  • 地表面エネルギー収支の基本
  • Linux マシン hamana にユーザー登録
  • 課題: 松山・増田 第5章を読む
• 第2回 4月21日(月) 2限 3号館 3410
  • 地表面エネルギー収支の基本 つづき
    • [松山・増田 (2021) 第5章]
  • 世界の地表面熱収支各項 (月ごとの平年値) の分布 (再解析データによる)
    • 「再解析データ」
      • [JRA-55]
      • [JRA-3Q]
      • (用語補足) [再解析 (reanalysis)] (2016-05-18)
      • (用語補足) [データ同化 (data assimilation)] (2022-06-25)
    • まず、(データが旧版だが)「JRA-55アトラス」(1981-2010 平年値) の図を見る。
      • [気象庁 「JRA-55 アトラス」]
        • 表示は「カラー1」「カラー2」「白黒」を選択可能、画像ファイルまたは数字テキストファイルをダウンロード可能
        • 地表面および大気上端の気候値
          • 地表面における正味の熱フラックス ... G (陸では地中熱伝導)
          • 地表面顕熱フラックス ... H
          • 地表面潜熱フラックス ... L E
          • 地表面における下向き太陽放射 ... R_S^↓
          • 地表面における下向き長波放射 ... R_L^↓
          • 地表面における正味の太陽放射 ... R_S^↓ - R_S^↑
          • 地表面における正味の長波放射 ... R_L^↓ - R_L^↑
        • 地表面の気候値 ... 気温 (地上 2 m)、降水量、風 (地上 10 m でのベクトル平均風) などがある。
    • JRA-3Q の平年値 (1991-2020) を、Linux 上の Python で作図してみる。(次回につづく。)
      • hamana の /home/masudako/ex/jra3q_clim9120/ に例文がある。
      • これまでに用意したプログラム
        • 1地点の季節変化の折れ線グラフ
          • graph_srb1.py ... 地表面の放射4成分
          • graph_seb1.py ... 地表面エネルギー収支項 Rnet, H, LE, G
          • (水収支項 準備中)
        • 各月の平年値の空間分布の地図表示
          • map_psl.py ... 海面気圧。エネルギー収支項への変更はあとで説明する。
      • 作業手順の例
        • cp -r /home/masudako/ex/jra3q_clim9120
        • mkdir image (一度だけ、画像ファイルがつくられる場所を確保しておく)
        • cd jra3q_clim9120
        • (ここからプログラムのファイル名や内容をかえてくりかえす)
        • emacs graph_seb1.py (編集する)
        • python graph_seb1.py
        • ls ../image
        • display image/seb40.00N110.00E.png (例)
      • 変更すべきところ (graphのプログラムについて)
        • ほしい緯度経度 plon, plat に代入する値
        • 図を画像ファイルに保存するか (savefig=True)、画面表示か (savefig=False)
        • 凡例をかく縦軸上の位置 ypos, dypos に代入する値
• 第3回 4月28日(月) 2限 3号館 3410
  • 画面表示用のコマンドについて [display と animate]
  • JRA-3Q の平年値 (1991-2020) を、Linux 上の Python で作図してみる。(つづき。さらに次回につづく。)
    • 例文プログラムを修正したので、コピーしなおしてください。
      • cd jra3q_clim9120
      • cp /home/masudako/ex/jra3q_clim9120/map_psl.py . 【スペース ピリオド】
    • 月ごとの地表面エネルギー収支各項の地理的分布を作図する。map_psl.py のプログラムを、つぎのところを変更して実行する。
      • プログラムファイル名は、気象要素名略号 (増田が便宜上きめたもの) をつかって、map_rsdown.py のようにつけよう。
      • 気象要素名略号は、[JRA-3Q の 平年値 (1991-2020) の月平均値のデータ] のページの表中の「ファイル名中の略号」。
        • rsdown, rsup, rldown, rlup, shf, lhf。
      • 図を画像ファイルに保存するか (savefig=True)、画面表示か (savefig=False)。
      • 読まれるファイルが置かれたところ (datapath)。例文は気圧なので「surf125」となっているところを、 エネルギー収支項では「phy2m125」と変更。
      • 「バイナリデータを読む」ところで、「psl」(海面気圧) を「rsdown」(下向き短波放射) に変更。 さらに、順次、見たい気象要素名略号に変更していく。
      • データを読んだあとの単位換算。 pslは気圧なので Pa → hPa の換算のために 100 でわっているが、他の気象要素ではこの操作をやめる。 (気象要素によってはそれに即した換算をする。たとえば気温ならば K → ℃の換算のために 273.15 をひく。)
      • 「等値線によってぬりわける」のところで、配列 levels (および colors) の内容。 【ここは、各気象要素ごとに数値の見当をつけたうえで、試行錯誤する。】
        • (cmap.set_under, cmap.set_over が無効になっているようであり、 対策として levels, colors を上下にひろげているので) mappable = の行の「extend='both',」をとりやめる。
      • 「カラーバーをかく」の cb.set_label で 単位の文字列を変更。'hPa' を 'W/m2' とする。【'W/m$^2$'として上つきがでるならばそうしたいが、未確認。】
      • 「タイトルをかく」の stitle に代入する文字列の 'Sea level pressure' を変更。 表の「Item」の列を参照。長すぎるときは省略形をくふうする。
      • 「図を画像ファイルに保存」の imagefilename の psl の部分を変更 (気象要素名略号にあわせるとよい)。
      • 表示する緯度経度範囲をかぎって図を拡大したいばあいは、ax.set_extent のコメントをはずし、 図の左端、右端、下端、上端の経度・緯度の数値を書きかえる。
    • 地表面の太陽放射反射率 (Surface albedo) alphas = rsup / rsdown を作図してみる。 (分母ゼロの対策が必要。) 2つのファイルを読み、演算してから作図する。
    • 地表面の正味放射 (Net radiation) rnet = rsdown - rsup + rldown - rlup を作図してみる。 4つのファイルを読み、演算してから作図する。
    • 地表面エネルギー収支の G の項 (陸では地中熱伝導) (暫定的に Ground heat flux とよぶ) ghf = rnet - shf - lhf を作図してみる。 6つのファイルを読み、演算してから作図する。
    • 「ボーエン比」(Bowen ratio) bowenr = shf / lhf を作図してみる。 (分母ゼロの対策が必要。) 2つのファイルを読み、演算してから作図する。
• 第4回 5月12日(月) 2限 3号館 3410
  • JRA3Q によるグローバルな量の作図 つづき
    • 地表面水収支
      • 各項は、単位時間・単位面積あたりの水の質量の流れ、SI単位は kg/(m² s)。
      • 液体の水の標準の密度をつかうと、1 kg/m² は 1 mm にあたる。
      • 降水量 precip, 蒸発量 evap, 降雪 snowfall (降水のうちわけ)
      • 降水量と蒸発量の差 (仮記号 pminuse)
      • 大気中の水蒸気の水平輸送の (水平) 収束量 column_qconv も見たいが今回は見送る。
      • 流出量は JRA3Q では陸面の予報値の変数グループ fcst_land125 にあるので今回は見送る。
    • 大気上端 (top_of_atmosphere, 英文中では TOA、略号中では toa と略す) の放射収支
      • toa_rsdown, toa_rsup, toa_rlup (toa_rldown は 0 なのでデータファイルはない)。
      • toa_rsnet = toa_rsdown - toa_rsup
      • toa_rnet = toa_rsnet - toa_rlup
      • 大気上端でみた太陽放射反射率 (planetary albedo) alphap = toa_rsup / toa_rsdown (alphas 用のプログラムから修正)
    • 月値から年平均値 (annual mean) を計算し、その水平分布を作図する。 [第4回には降水量と蒸発量。次回につづく。]
      • 月から年への平均操作と、 作図とを、別のプログラムでやる (平均操作の結果をファイルに書き出す) か、 1つのプログラム中でやるか、方針をきめてから、プログラムをつくる。
        • 1つのプログラム中でやることにした。
        • 例文 (海面気圧) hamana の /home/masudako/ex/jra3q_clim9120/map_psl_ann.py
      • 月の日数が一定でないことをどう処理するか?
        • 一定とみなしてしまう人もいる。
        • ここではグレゴリオ暦の日数にしたがう。ただし2月は 28.25 日とする。
• 第5回 (5月19日(月) 2限 3号館 3410)
  • JRA3Q のエネルギー収支・水収支各項の図 (とくに年平均) を見て (必要に応じてつくりなおして)、世界の分布の特徴をつかむ。
    • 地表面水収支の各項について、 年平均値の東西全経度平均値 (zonal mean) を計算し、その南北分布を作図する。
      • 例文 hamana の /home/masudako/ex/jra3q_clim9120/lat_evap_ann_zon.py (蒸発量)
    • 格子点の地表面状態分類のファイルも読んで、海と陸にわけて平均し、作図してみる。
      • 例文 hamana の /home/masudako/ex/jra3q_clim9120/lat_evap_ann_zon_landsea.py (蒸発量)
    • プログラムを改造して、蒸発量と降水量を同じわくに作図してみる。
    • 大気上端の放射収支の各項について、 年平均値の東西全経度平均値 (zonal mean) を計算し、その南北分布を作図する。
      • 例文 hamana の /home/masudako/ex/jra3q_clim9120/lat_toa_rad_ann_zon.py (下向き短波、上向き短波、上向き長波)
      • 同 lat_toa_rad_ann_zon_2.py ([上記の項から計算した] 正味下向き短波、上向き長波、正味放射)
    • 【今回のプログラム例は、 複数の気象要素をいっしょに表示するばあい、 ほぼ同じでわずかにちがう作業をくりかえす形の、うまくないプログラムになってしまった。 今後の改良としては、年平均の計算と東西平均の計算をそれぞれ下請け関数にして、 引き数となる変数だけ変えてよびだす形にしたほうがよい。】
    • 教科書的文献と対応づける。
      • Hartmann (2016) 『Global Physical Climatology, 2nd edition』の 2章 (大気上端の放射収支)、4章 (地表面エネルギー収支)、5章 (水収支) の図を参考にする。
• 第6回 (5月26日(月))、第7回 (6月2日(月))
  • 近藤 (2000) 『地表面に近い大気の科学』
    • まず 1.8 節「熱輸送量の観測」を読む。 (その必要に応じて 1.1～1.5節に書いてあることも確認する。)
    • つづいて、第5章 5.1節「熱収支量の大きさ」と5.2節「地表面の熱収支式」を読む。 とくに「バルク式」は何を入力データとして必要とするかを確認する。
    • 陸の地中熱伝導について、熱伝導式、表4.1の土壌の熱的パラメータの代表的な値、 4.4節後半および5.1節の「参考5.1」に出てくる、境界条件を時間の三角関数であたえたばあいの 地中温度の変化の計算例を見る。計算例を追加する予定。
• 第8回 (6月9日(月))
  • 陸の地中熱伝導について
    • [注] 地面から下をあつかうとき、z軸は下向き (値が正) と、上向き (値が負) の両方の流儀がある。 そのときごとにどちらであるか確認が必要。
    • 熱伝導式と土壌層のエネルギー収支式を確認する。(近藤 2000 の本にはあまり明確に書いてない。)
      • zの正の向き、Gの正の向きをどうとるかきめたうえで、符号を確認すること。
      • 熱伝導のエネルギーの流れ [W/m²]: G = -λ_G ∂T/∂z (p. 113 の (4.3) 式に z = 0 のばあいについて書いてある。) (熱は温度の高いほうから低いほうにむかう。)
      • 厚さΔz の層のエネルギー収支: ∂(ρ_G c_G T Δz)/∂t = -{G(z+Δz) - G(z)}
      • ρ_G, c_G, λ_G を定数とみなせるならば (みなせないこともある)、Δz → 0 で ρ_G c_G ∂T/∂t = -∂G/∂z
      • 熱伝導式から G を消去すると ∂T/∂t = +(λ_G/ρ_G c_G) ∂²T/∂z²
    • 土壌の熱的パラメータ (密度、比熱容量、熱伝導率)
      • いくつかの地表面物質の代表的な値 ... 4.2節 表4.1 p.115
      • 土壌水分 (体積含水率) への依存性をふくむ値 ... 4.4節 (表4.3ほか) p. 124-125
    • 境界条件として地表面温度を時間の三角関数であたえたばあいの日変化 (「解析解 = 微分積分による解」の地温) 4.4節後半 p. 121-122
    • この「解析解」の数値例を計算するプログラム例
      • cp -pr /home/masudako/ex/tg/ .
      • cd tg
      • emacs tg_model_1.py ... プログラムを見る
      • capacg: 土壌の体積あたりの熱容量 ρ_G c_G [J m^-3 K^-1]
      • lambdag: 土壌の熱伝導係数 λ_G [W m^-1 K^-1]
      • (camacg, lambdag の現状は、仮に、黒ボク土、体積含水率 0.2 を想定した概算値を入れてある。)
      • tsmean, tsamp1: 地表面温度の日平均 T_M と1日周期の振幅 A₁ [K]
      • deltathr: 時間ステップ Δt [h]; nstep: ステップ数
      • deltaz: 深さきざみ Δz [m]; ndepth: 鉛直レベル数
      • outfilename: 結果を書き出すテキストファイル名
      • python tg_model_1.py
      • 結果 tg_model_1_out.txt 横が深さ、縦が時間
      • python transpose.py ... 計算結果のファイルの縦横をいれかえる
      • 結果 tg_model_1_timex.txt 横が時間、縦が深さ
      • python to_ztv.py ... 計算結果を 深さ、時間、温度の3列の表にする
      • 結果 tg_model_1_ztv.txt
      • gnuplot で作図
        • load "tg_model_1_bytime.gplt" ... 時間を横軸にとったグラフ
        • load "tg_model_1_bydepth.gplt" ... 深さを縦軸にとったグラフ
        • load "tg_model_1_contour.gplt" ...時間を横軸、深さを縦軸にとった等値線 (contour) 表示 (等値線の数値ラベルが まだ)
  • 今後の発展 (キャンパスの観測値によるエネルギー収支項の計算をしたあと、時間があればやる。)
    • 地温だけでなく、熱伝導によるエネルギーフラックスの量も見る。
    • 解析解でなく、時間差分による「数値解」を計算する。 (きれいな境界条件でないばあいにも応用できる。)
    • 日変化でなく年変化のばあい。空間領域 (deltaz × ndepth) をどのように大きくするか?
    • 境界条件が、地表面温度をあたえるのではなく、 エネルギーフラックスの値があたえられるばあい。
  • 関連する旧教材ページ (プログラム言語が おもに Fortran 77)
    • [ローカルな地表面エネルギー収支を表現した0次元・1次元モデル]
    • (このうちのプログラム kakusan を Python になおしたものが /home/masudako/ex/kumagaya_campus/ にいれてある。)
• 第9回 (6月16日(月))
  • 陸の地中熱伝導について
    • 境界条件として地表面温度を時間の三角関数であたえたばあいの日変化 (「解析解 = 微分積分による解」の地温) 4.4節後半 p. 121-122
      • この「解析解」が解になっていることを確認する。数学の知識として「合成関数の微分」と「積の微分」をつかう。
        • [手書きページ 1]
        • [手書きページ 2]
  • 「ローカルな観測にもとづく熱収支各項の値を見る」へ
    • 手はじめに、放射フラックス (4成分) を時系列グラフで見る。 日変化と、日々の天気によるちがいを見る。 季節によるちがいについてもいくつかの例とって見る。
    • 地表面の太陽放射反射率 (アルベド) を計算してみる。
    • 正味放射を計算してみる。
• 第10回 (6月23日(月))
  • 「ローカルな観測にもとづく熱収支各項の値を見る」つづき
    • 2018年 1月 23-24日のアルベド (計算結果) とその日の気象状況を確認する。
    • 気温と、地表面温度 (放射温度計による) をみる。
    • 地表面温度 (放射温度計による) の黒体放射を計算し、上向き長波放射と比較する。
    • 湿度をみる。湿度は相対湿度であたえられているが、比湿を計算してみる。
    • 風速・風向をみる。
    • 気温の観測高さを推測する。
    • バルク法の中間変数である、地表面温度と気温との差、地表面温度での飽和比湿と大気の比湿との差、を計算する。

ローカルな観測にもとづく熱収支各項の値を見る

熊谷キャンパス露場の「放射・熱収支観測システム」 (渡来ほか 2010, 中川ほか 2010) の観測値によって、計算・作図してみる。

データと Python プログラム例は、Linux マシン "hamana" に置く。 作業はそこでしても別のところでしてもよいが、プログラムは "hamana" でも動くものがのぞましい。

文献

• 渡来 靖、丸本 美紀、福岡 義隆、中川 清隆, 2010: 立正大学総合気象観測システムおよび放射・熱収支観測システムの拡充について。 地球環境, 12: 139-145.
• 中川 清隆, 渡来 靖, 福岡 義隆, 2010: 埼玉県熊谷市南郊における冬季静穏快晴日の 放射 4成分・地表面温度・気温の日変化の事例解析。 地球環境, 12: 77-84.
  • その図を [環境気象学概論「放射」] の教材ページにも引用している。

データファイル (「放射・熱収支観測システム」の観測値)

2021年 9月末に渡来先生からいただいたファイル。 CSV形式で、2015年10月1日 10:40 JST から 2021年 9月29日 15:10 JST までの10分値。 ファイル名 surface_TBL110.dat だった。

hamana に /raid/kumagaya_campus/obsdata/surface_TBL110_20151001_20210929.dat として置いた。

185715行め。レコード番号 187003 の不完全な内容に 187005 の完全な内容がつづいている。 187003 の不完全な内容をとりのぞいた。(とりのぞくまえのものを bak のついたファイル名で残した。)

それぞれの列が何をさすかは 下にしめす head.txt を 渡来ほか (2010) を参照して、みてほしい。 Avg は 10分平均値、Smp は瞬間値だろう。 ただし、列見出しに不適切なものがある。たとえば、CNR_UP_S_Avg は、内容からみて下向き短波放射である。 (上を向いている放射計による観測値だからこうなったか?)

プログラム例

hamana の /home/masudako/ex/kumagaya_campus/ に置いた。

• readhead.py ... CSVファイルの最初の4行 (4行めは最初の観測値データ) を読み、書き出す。 結果は head.txt 。
• readsel1.py ... pandas の データフレーム (df) に読みこむ。 (上記の不完全なレコードを修正しないとエラーになった。)
• readsel2.py ... 対象とする日時の範囲 (dfの行) を限定する。
• readsel3.py ... さらに、観測項目 (dfの列) を限定する。(ひとまず放射だけにした。) 結果を CSV ファイル readsel3out.csv に書き出す。
• graphsel3.py ... readsel3out.csv を読み、放射4成分の時系列グラフをかく。
• graphsel3x.py ... グラフを2段にして、下に地表面アルベドのグラフをかく。 結果を 画像ファイル (例、graphsel3x20151001.png) に書き出す。

3号館1階のモニター画面に表示されているもの

• 放射4成分 CNR_UP_S_Avg, CNR_DN_S_Avg, CNR_UP_LC_Avg, CNR_DN_LC_Avg
• 気温 Ta_Avg(2) [15 cm]、Ta_Avg(3) [90 cm] ←【高さは下にのべたように数値例から推定したもの】
• 地表面温度 (赤外放射温度計) IRT_Avg(1)
• 地温 Ts_Avg(5) [2 cm], Ts_Avg(6) [5 cm], Ts_Avg(7) [10 cm], Ts_Avg(8) [30 cm]
• 土壌水分 [%] SWC_per(3) [5 cm], SWC_per(4) [20 cm]
• 風速 Wind_Spd_WVc(1), 風向 Wind_Spd_WVc(2)

フラックス計算では、このほかに 相対湿度 Rh_Avg をつかう。気温 Ta_Avg と同じ番号が同じ場所と思われる。

今後の計画

• 放射4成分、地表面の反射率 (アルベド)、正味放射をみる。
• 地表面温度をみる。上向き長波放射を地表面温度の黒体放射とくらべる。
• 気温、湿度、風速をみる。湿度は相対湿度であたえられているが、比湿を計算してみる。
• バルク法で、顕熱・潜熱フラックスを計算してみる。
  • 気温とつじつまのあった高さでの風速をつかって計算すべきだが、 暫定的に、観測された風速をそのままいれて計算してみる。
• 土壌の物理パラメータにありそうな値をあたえて、地中熱伝導フラックスを計算してみる。 方法はつぎのいずれか。
  • ある時刻の地温の鉛直プロファイルから、熱伝導式で計算する。
  • 複数の時刻の地温の鉛直プロファイルからエネルギーのたまりの量の時間変化を計算し、そこからフラックスを知る。
• いろいろな季節や総観的状況のもとでのエネルギー収支各項をみる。

気温・湿度の観測位置の確認のこころみ

• readselta.py ... Ta 4列と Rh 4列をとりだす
• readseltaout.csv ... とりだされた数値のCSVファイル
• graphta.py ... Ta と Rh のグラフをつくる。 凡例がまだないが、red, green, blue, grey の順。
• graphta20151001.png ... Ta と Rh のグラフ。

グラフをみると青だけが一定方向にずれているので、(3) が 高さ 90 cm、 (1) (2) (4) が 15 cm での観測と推測される。

過去の教材ページ

(リンクさきがなくなっているばあいや、別のものにおきかわってしまうばあいもあります。)

[2024年度教材]

[2023年度教材]

関連の教材ページ

[圏間相互作用特論 II]

[環境気象学アドバンスト (増田担当分) 2021年度]

[気圏環境学演習 (増田担当分) 2023年度]

[計算機を利用した気象分野の実習の教材]