1．はじめに
前回の記事では，蒸発散比法(仮称)を使って日単位で流域実蒸発散量を推定してみました．
今回は補完法という方法を使って，日単位で流域実蒸発散量を推定してみたいと思います．
ただ，日単位で流域実蒸発散量の推定に補完法を使うのが果たして適切かどうか自信はないです． (ざっと調べた限りではおそらく大丈夫だと思いますが．．．) なので，今回もあくまで学生や若手技術者の方々の学習の支えとなることを目的に記事を書くことにします．
とはいえ，できるだけ正確な内容にはしたいと思うので，もし内容の誤り等がありましたら優しくご指摘いただけると幸いです．
ちなみに，今回の記事はペンマン法について理解していることを前提に書いているので，ペンマン法についてよく知らない方は下記の記事をご覧ください．
【流出解析】シリーズ 一覧
【流出解析】シリーズ 第1回 流域界の決定と入力データの用意 【流出解析】シリーズ 第2回 貯留関数法【前編】 【流出解析】シリーズ 第3回 貯留関数法【後編】 【流出解析】シリーズ 第4回 流域実蒸発散量の推定―蒸発散比法― 【流出解析】シリーズ 第5回 実蒸発散量の推定―補完法― 【流出解析】シリーズ 第6回 タンクモデル 【流出解析】シリーズ 第7回 積雪・融雪モデル ―降雪量・融雪量の推定― 【流出解析】シリーズ 第8回 タンクモデルの最適化 【流出解析】シリーズ 第9回 単位図法 【流出解析】シリーズ 第10回 流出量・流出高・比流量の違い 【流出解析】シリーズ 第11回 表面流モデル【定常降雨：前編】 【流出解析】シリーズ 第12回 表面流モデル【定常降雨：後編】
2．補完法とは
補完法とは，「ペンマン法による蒸発散位は地域の乾湿状態に応じて実蒸発散量に対して補完的に変化し，両者の和は可能蒸発量の2倍に等しい」という観測結果に基づいた方法です．
ここで，
$$ $$ 「補完法は観測結果に基づく方法です」という説明では，さすがに雑すぎるのでもう少し具体的にイメージしてみましょう．
ペンマンの蒸発散位は次式で表されます．
ここで，
ペンマン法で扱う水蒸気圧\( e \)は，一般的な気象観測が行われる高さ(地上数十センチメートル〜数メートル程度)における値です． 地面が乾燥しているとき水蒸気圧\( e \)は低下し，逆に湿っていれば\( e \)は上昇します．
また，ペンマンの蒸発散位\( ET_{p} \)の第2項は飽差に比例します． つまり，乾燥条件下では飽差が大きくなり蒸発散位も増加します．逆に湿潤条件下では，蒸発散位は減少します．
一方で，実蒸発散量は地面が乾燥していれば減少し，湿っていれば増加します． 洗濯物を干した直後は水分がどんどん抜けていくのに対し，乾ききる直前の水分が抜けづらいのと同じですね．
| 乾湿状態 | 蒸発散位ETp | 実蒸発散量ETa | 合計(ETp+ETa) |
|---|---|---|---|
| 乾燥 | 増 | 減 | 一定 |
| 湿潤 | 減 | 増 | 一定 |
蒸発散位\( ET_{p} \)と実蒸発量\( ET_{a} \)が乾湿条件に応じて補完的に変化し，その合計が一定であるという性質から補完法と呼ばれるわけです．
3．計算方法
補完法の概要と具体的なイメージは掴めたでしょうか？
それでは，実際に計算していきたいと思います． 今回は大槻先生らの論文を参考に補完法を用いて流域実蒸発散量を推定してみます．
大槻恭一, 三野徹, 丸山利輔, 1984. 水収支法と補完関係式による流域蒸発散量の比較 実蒸発散量推定に関する研究（II）. 農業土木学会論文集. 112, 17-23.
補完法は次式で表されます．
ここで，
それでは，それぞれの式について細かく見ていきたいと思います．
はじめに，可能蒸発散量\( ET_{pot} \)について． 可能蒸発散量\( ET_{pot} \)は，移流(水平方向の物理量の輸送)のない広大で均一な湿潤面で生じる蒸発散量のことです．
可能蒸発量は平衡蒸発量(飽差が0の時の蒸発散位．ペンマン式の第一項に相当．)の1.26倍になることが知られているそうです． 式(2)の\( {ET_{pot}}^{''} \)は移流の影響を考慮した可能蒸発量です． ちなみに式(3)は，移流の影響を考慮したペンマン式です．
式(4)は，実蒸発散量が蒸発散位を超えることがないようにする制約条件です．
移流項Mの推定
今回の計算で最も意味不明なものは移流項\( M \)です． 幸い，色々と論文を読んだら\( M \)の推定式が見つかりました．
移流項\( M \)の推定式
ここで，
純放射量\( Rn \)の計算には，ペンマンの蒸発散位を計算するときに使った推定式を使いまわしています． 雲量\( C \)の値は0〜1になります．
4．結果
計算用のワークシート↓
実蒸発散量がマイナスということは結露かな？
5．おわりに
最後までお付き合いいただきありがとうございました．
補完法はもともと，十分な気象データが手に入らなかった，あるいは気象データの蓄積が十分でない時代に何とかして実蒸発散量を推定しようとして作られたそうです． 補完法は今はあまり使われることはなく，最近は蒸発散比法(前回の記事を参照)がよく使われるようです．
とはいえ，気象庁が公開している気象データだけで蒸発散量が求められるというのは魅力的ですね．
今回は，学生や若手技術者の学習の支えとなることを目的としているので，補完法の適用範囲など込み入った話には特に触れていません． 実際に使う場合にはその辺のことについてよく調べるようにしましょう．
あくまで，この記事は「講義で補完法って習ったけど使い方もよくわからんし，具体的な計算方法もわからんから，どう勉強していいかわからん」という人を対象としています．
ただ，記事に誤りがあった場合にはすぐに訂正させていただくので，お気づきの点がございましたら優しくご指摘いただけると嬉しいです．
記事の最後に計算用のプログラムを掲載しています． ご自由にお使いください．
【流出解析】シリーズ 一覧
【流出解析】シリーズ 第1回 流域界の決定と入力データの用意 【流出解析】シリーズ 第2回 貯留関数法【前編】 【流出解析】シリーズ 第3回 貯留関数法【後編】 【流出解析】シリーズ 第4回 流域実蒸発散量の推定―蒸発散比法― 【流出解析】シリーズ 第5回 実蒸発散量の推定―補完法― 【流出解析】シリーズ 第6回 タンクモデル 【流出解析】シリーズ 第7回 積雪・融雪モデル ―降雪量・融雪量の推定― 【流出解析】シリーズ 第8回 タンクモデルの最適化 【流出解析】シリーズ 第9回 単位図法 【流出解析】シリーズ 第10回 流出量・流出高・比流量の違い 【流出解析】シリーズ 第11回 表面流モデル【定常降雨：前編】 【流出解析】シリーズ 第12回 表面流モデル【定常降雨：後編】
今後も農業農村工学(水文学，かんがい排水，土壌物理，水理学)を中心に記事を執筆していきたいと思います． リクエスト等も受け付けておりますので，ご遠慮なく連絡ください． Twitterアカウント：エビぐんかん@6LxAi9GCOmRigUI メール：nnCreatorCircle@gmail.com
検索用Key Words： 実蒸発散量，推定，水文学，流域，蒸発散位，ペンマン，補完法，移流
引用・参考文献
(1)中道丈史, 諸泉利嗣, 三浦健志, 2011. 補完関係式を用いた実蒸発散量推定式の改良. 岡山大学環境理工学部研究報告. 16-1, 35-39. (2)三浦健志, 奥野林太郎, 1993. ペンマン式による蒸発散位計算方法の詳細. 農業土木学会論文集. 164, 157-163. (3)佐渡公明, 1992. 広域蒸発散量に及ぼす土地被服の影響について. 土木学会北海道支部 論文報告集 平成4年度. 623-628. (4)松井宏之, 2010. ペンマン型蒸発散量推定式における有効長波放射量推定式の比較. 農業農村工学会論文集. 78-6, 531-536. (5)佐渡公明, 小島正洋, 1996. 補完関係式による竜式実蒸発散量年変化の推定. 水工学論文集. 40, 39-334. (6)菅原正巳, 面積雨量の求め方についての数学的考察. 水利科学. 112, 23-43. (7)大槻恭一, 三野徹, 丸山利輔, 1984. 計器蒸発量, 蒸発散位と実蒸発散量の関係 実蒸発散量推定に関する研究（I）. 農業土木学会論文集. 111, 95-103. (8)大槻恭一, 三野徹, 丸山利輔, 1984. 水収支法と補完関係式による流域蒸発散量の比較 実蒸発散量推定に関する研究（II）. 農業土木学会論文集. 112, 17-23. (9)大槻恭一, 三野徹, 丸山利輔, 1984. 気象資料から推定したわが国の蒸発散量 実蒸発散量推定に関する研究（III）. 農業土木学会論文集. 112, 25-32. (10)F.I. Morton，1978．Estimation evaporation from potential evaporation :practicality of an iconoclastic approach．Journal of Hydrology．38，1-32． (11)丸山利輔，石川重雄，大槻恭一，高瀬恵次，永井明博，田中丸治哉，駒村正治，赤江剛夫，堀野治彦，三野徹，武田育郎，金木亮一，渡辺紹裕：地域環境水文学，朝倉書店，2014，pp.29-39
6．計算用プログラム
ご自由にお使いください．
VBAの使い方はこちら
【VBA】Excelでプログラミングできる環境を構築し簡単なプログラムを組んでみる【プログラミング】 VBA の使い方、開発環境を整える
Sub ActualEvapotranspiration_ComplementaryMethod()
'大槻ら，1984の論文に書かれている実蒸発散量推定式の実装プログラム
'大槻恭一，三野徹，丸山利輔，1984．水収支法と補完関係式による流域蒸発散量の比較―実蒸発散量推定に関する研究(?)―．農土論集．112，17-24．
'ペンマン法により蒸発散位を計算するプログラム(2019/12/19作成)を一部改変して作成した．
'作成日：2020/11/13
'圧力の単位mbarはhPaに等しいため，本プログラムではhPaと表記することにした
'参考文献1：三浦健志, 奥野林太郎, 1993. ペンマン式による蒸発散位の計算方法の詳細. 農業土木学会論文集. 164, 157-163.
'参考文献2：丸山利輔, 三野徹 編, 地域環境水文学, 朝倉書店, p.18, pp.36-37(2014)
'参考文献3：文字信貴, 平野高司, 高見晋一, 堀江武, 桜谷哲夫 編, 農学・生態学のための気象環境学, 丸善株式会社, p.8(2004)
'有効数字については特別考慮していないので実際に利用する場合は注意すること
'画面更新の停止
Application.ScreenUpdating = False
'変数の定義
Dim i As Integer, j As Integer 'カウント用
Dim lat As Double '緯度(°)
Dim latR As Double '緯度(rad)
Dim H As Double '風速計地上高度(m)
Dim A As Double 'アルベド
Dim TD As Integer '積算日数(Total Days) (日)
Dim dd As Double '地球太陽間の距離(天文単位)
Dim Dcl As Double '赤緯(Declination) (°)
Dim DclR As Double '赤緯(Declination)をラジアンで表現
Dim w As Double '日没時の時角(rad)
Dim N As Double '可照時間N(h)
Dim esat As Double '飽和水蒸気圧(hPa)
Dim ea As Double '水蒸気圧(hPa)
Dim dlt As Double '飽和水蒸気圧曲線の勾配(hPa/℃)
Dim l As Double '蒸発潜熱(MJ/kg)
Dim u2 As Double '地上2mにおける風速(m/s)
Dim fu As Double '風速関数
Dim ET0 As Double '蒸発散位(mm/d)
Dim Gamma As Double '乾湿計定数(hPa/℃)
Dim Sigma As Double 'ステファン-ボルツマン定数(W/m2/K4)
'補間法で使用する変数の定義
Dim C As Double '雲量
Dim ε As Double '射出率
Dim ρ As Double '大気長波放射量と快晴時の大気長波放射量の比
Dim B As Double '有効長波放射量(MJ/m2/d)
Dim M As Double '移流項(MJ/m2/d)
'データ出力欄の作成
Range("H1") = "出力データ"
Range("H1:M4").Merge
Range("H1:M4").HorizontalAlignment = xlCenter
Cells(5, 8) = "赤緯(°)"
Cells(5, 9) = "大気圏外日射量(MJ/m2/d)"
Cells(5, 10) = "純放射量(MJ/m2/d)"
Cells(5, 11) = "蒸発散位ET0第1項 (mm/d)"
Cells(5, 12) = "蒸発散位ET0第2項 (mm/d)"
Cells(5, 13) = "蒸発散位ET0 (mm/d)"
'地点パラメータの取得
lat = Cells(3, 3) '緯度(°)
latR = WorksheetFunction.Radians(Cells(3, 3)) '緯度(rad)
H = Cells(3, 4) '風速計地上高度(m)
A = Cells(3, 5) 'アルベド
'物理化学定数(参考文献1の値を使用)
Gamma = 0.66 '乾湿計定数(hPa/℃)
Sigma = 4.9 * 10 ^ (-9) 'ステファン-ボルツマン定数(W/m2/K4)
ε = 0.92 '射出率
'ペンマン法の本体部分(参考文献1のpp.160-162)
i = 6
Do
'************ペンマン法************
'起算日(1/1)からの積算日数(日)
TD = Cells(i, 2) - DateSerial(Year(Cells(6, 2)), 1, 1) + 1
'太陽からの距離(天文単位)
dd = 1 + 0.01676 * Cos(WorksheetFunction.Radians(0.977 * (TD - 186)))
'赤緯(°)
Cells(i, 8) = 23.45 * Cos((TD - 173) * 0.966 * 3.1416 / 180)
'赤緯(rad)
DclR = WorksheetFunction.Radians(Cells(i, 8))
'日没時の時角(rad)
w = WorksheetFunction.Acos((-Tan(latR) * Tan(DclR)))
'可照時間(h)
N = 2 * WorksheetFunction.Degrees(w) / 15
'大気圏外日射量(MJ/m2)
Cells(i, 9) = 1.37 * 10 ^ (-3) / dd ^ 2 * 86400 / (4 * Atn(1)) _
* (w * Sin(latR) * Sin(DclR) + Sin(w) * Cos(latR) * Cos(DclR))
'飽和水蒸気圧(hPa)
esat = 6.1078 * Exp(17.2694 * Cells(i, 3) / (Cells(i, 3) + 237.3))
'水蒸気圧の計算(hPa)
ea = esat * Cells(i, 4) / 100
'純放射量(MJ/m2)(参考文献1のステファン-ボルツマン定数の表現がちょっとおかしい)
Cells(i, 10) = (1 - A) * Cells(i, 9) * (0.18 + 0.55 * Cells(i, 6) / N) - Sigma * (Cells(i, 3) + 273.2) ^ 4 _
* (0.56 - 0.092 * 0.866 * Sqr(ea)) * (0.1 + 0.9 * Cells(i, 6) / N)
'飽和水蒸気圧曲線の勾配(hPa/℃)
dlt = 0.4495 + 0.2721 * 10 ^ (-1) * Cells(i, 3) + 0.9873 * 10 ^ (-3) * Cells(i, 3) ^ 2 + 0.2907 * 10 ^ (-5) * Cells(i, 3) ^ 3 + 0.2538 * 10 ^ (-6) * Cells(i, 3) ^ 4
'蒸発潜熱(MJ/kg)
l = 2.5 - 0.0024 * Cells(i, 3)
'地上2m風速への換算と風速関数fu
u2 = Cells(i, 5) * Log(200) / Log(100 * H)
fu = 0.26 * (1 + 0.54 * u2)
'蒸発散位(mm/d)
Cells(i, 11) = dlt / (dlt + Gamma) * Cells(i, 10) / l
Cells(i, 12) = Gamma / (dlt + Gamma) * fu * (esat - ea)
Cells(i, 13) = Cells(i, 11) + Cells(i, 12)
'************補完法************
'雲量Cの推定
'※雲量は0〜1になるようにする
'雲量がマイナスになるのは可照時間Nの推定誤差による？
If Cells(i, 6) <= N Then
C = (1 - Cells(i, 6) / N)
Else
C = 0
End If
'ρの計算
ρ = 1 + (0.25 - 0.005 * (esat - ea)) * (C ^ 2)
'Bの計算
B = ε * Sigma * ((Cells(i, 3) + 273) ^ 4) * (1 - ρ * (0.707 * ea / 158)) * 86400 / (10 ^ 6)
'Mの計算
M = 0.66 * B - 0.44 * Cells(i, 10)
'流域実蒸発散量
'Ep"
Cells(i, 15) = 1.26 * dlt / (dlt + Gamma) * (Cells(i, 10) + M) / l
Cells(i, 16) = dlt / (dlt + Gamma) * (Cells(i, 10) + M) / l + Gamma / (dlt + Gamma) * fu * (esat - ea)
If (2 * Cells(i, 15) - Cells(i, 16)) <= Cells(i, 16) Then
Cells(i, 17) = 2 * Cells(i, 15) - Cells(i, 16)
Else
Cells(i, 17) = Cells(i, 16)
End If
i = i + 1
Loop Until Cells(i, 2) = ""
'合計と平均
Cells(i, 1) = "合計"
Cells(i + 1, 1) = "平均"
For j = 3 To 15
Cells(i, j) = WorksheetFunction.Sum(Range(Cells(6, j), Cells(i - 1, j)))
Cells(i + 1, j) = WorksheetFunction.Average(Range(Cells(6, j), Cells(i - 1, j)))
Next j
'画面更新の再開
Application.ScreenUpdating = True
End Sub
