第1章 はじめに
世界中にはさまざまな形態の地表面状態（水面（海面、湖面、川面））、雪面、裸地面、草地、森林、人工物（コンクリート、アスファルトなど）があり、それぞれの表面における熱・物質（水、二酸化炭素）交換の特徴の違いによって局地的スケールから全球スケールのエネルギー・物質循環が形成されている。こうした様々な地表面上において熱・物質循環輸送過程を調べることは気候形成の理解、環境問題の改善処置などのために非常に重要だと考えられている。過去においてさまざまな地表面上の熱・物質循環フラックスについて研究が行われてきて、近年では地球温暖化問題に関連して草地、森林などの陸域生態系が熱・物質循環過程に与える影響は大きいと考えられている。
熱・物質循環過程を計算で求める陸面過程モデルには現在までに様々なものがある。Sellers et al.(1997)によると、陸面過程モデルはその年代やモデルの構造から大きく分けて第一世代から第三世代に分けられる。
第一世代は陸面過程モデルとして最初に考えられたバケツモデル（Manabe,1969）である。バケツモデルは陸水貯留量の上限値を全球にて同じ深さとされている。しかし簡便さと引き替えに、植生の作用が考慮されていないので、それを考慮した結果と比較して日変化がおかしい、複雑さを表現できていない等の批判がある。また流出が簡単すぎて、水文過程を再現できていない。
第二世代にはBATS（Dickinson et al. 1993）、SiB（Sellers et al. 1986）がある。これらは植生の作用が考慮され、地表面植生と背の高い植生（以下キャノピーとする）が分離されている。また水文過程が改良され、放射や積雪の影響も考慮されている。その結果、日変化や降水の再現性は向上した。しかし、その改良で第一世代に比べて気候シュミレーションがどれぐらい良くなったかは依然として疑問点である。またこの段階では実際の植生で行われている光合成などの生物過程は取り入れられていなかった。
その後の第三世代にあたるSiB2（Sellers et al. 1996）、LSM（Bonan，1996）、MATSIRO（Takata et al. 2003）では、光合成の影響を考慮した炭素循環過程が導入されている。光合成を行う反応経路が異なる植物を分けて扱い、それぞれについて生物化学的な光合成-気候モデルを適用している。
本研究では先に挙げたモデルの中で、SiBを取り上げる。SiBはモデルが作られてから長い間使用されてきた。今でも多くの場所で用いられていて、最近まで気象庁でもSiBモデルを使用していた。また最近の第三世代のモデルの多くはSiBモデルをもとに作成されている。よって陸面過程モデルの中心であるSiBについてどのようなモデルの構成であるかをSellers and Mintz(1985)で詳細に調べる。そしてSiBモデルを実際計算できるように実装し、計算機器として使用するのが最も簡易であるExcelを用いて計算して、熱・物質循環を表すことができるか調べる。
第2章 SiBの概要
2.1 モデルの概念
太陽エネルギー（以下 可視光）と微量の気体の出す目に見えない長波放射（以下 近赤外域の放射）の2つの放射は植生空間を通って減衰した放射を植生の葉や地表面が吸収しエネルギーを得る。そのエネルギーは、水蒸気を蒸発したり植生の気孔から蒸散することによって使われるエネルギー（潜熱）と大気を直接暖める熱（顕熱）によって費やされる。
降水は植生に遮断され、遮断されなかった水分は地表面に到達し土壌に潜入していく。また植生や土壌の表面からの直接蒸発や植生の根から土壌の水分を根から吸い上げ気孔から蒸散することによって、大気中の水蒸気となる。
このような熱、水分の循環過程をSiBモデルによって得ることができる。
植生が吸収した放射量と顕熱、潜熱の作用によって、植生の温度は変化する。植生の水分量は葉による降水の遮断、排水と植生表面からの直接蒸発により変化する。土壌の含水率は、土壌内の水分移動、土壌表面からの直接蒸発、植生の根による水分の抽出により変化する。
これらの時間変化を得る事で、熱・物質循環過程を知ることができる。
また、SiBモデルでは背の高い樹木やかん木をまとめてキャノピ－と呼ぶ。
2.2モデルの構造
2.2.1 放射
樹葉の光学的特性と樹葉同士の反射により、樹木は0.4～0.72μｍの可視光をよく吸収し、赤外領域に近い0.72～4.0μｍに対しては適度な反射特性を示す。これに対し、植生の無い地表面では0.4～4.0μｍの近赤外域で反射率は波長とともに少しずつ増加する（図１）。
図１：植生と乾燥、湿潤裸地面のアルベド（地面の反射率）
そこで入射放射フラックスは、波長0.72μｍと4μｍそれぞれで区切られ分類される。また、直達光による放射と、拡散による放射でも分類され、全部で5つの入射放射フラックスに分けられる。
可視光（PAR）（＜0.72μｍ） 直達光による放射 Fs,b(0)
可視光（PAR）（＜0.72μｍ） 拡散による放射 Fs,d(0)
近赤外域（0.72～4μｍ） 直達光による放射 Fn,b(0)
近赤外域（0.72～4μｍ） 拡散による放射 Fn,d(0)
熱赤外線（＞4μｍ） 拡散による放射 Ft,d(0)
放射フラックスについては3.2章で示す。キャノピー、地表面によって吸収された放射フラックスではこの5通りに分けて計算する。
2.2.2 潜熱、顕熱輸送の構成
顕熱フラックス、潜熱フラックスを得る方法として、顕熱、潜熱輸送の流れを電気回路のように考える（図2）。よって、電気回路での基本的な公式のオームの法則（電流＝電圧÷抵抗）が用いられる。電流が顕熱、潜熱フラックスに、抵抗は熱が葉の気孔や土壌の表面や大気中の空間を通る時の抵抗に、電圧は顕熱の場合は温度差で潜熱の場合は水蒸気圧の差になる。
図2 潜熱、顕熱輸送の概念図。キャノピーの
右が顕熱、左が潜熱に関する図である。
図2で用いる文字は次の意味である(添字のｃはキャノピー、ｇは草、ｓは土壌表面、gsは土壌表面と草の両方を表す)。
H：顕熱フラックス Wｍ－２
λE：潜熱フラックス Wｍ－２
T：温度 K
ｆh：土壌表面での空気の相対湿度
ｒa：キャノピー空間と参照レベル間の抵抗 ｍｓ
：キャノピー空間とキャノピー葉面間の抵抗 ｓｍ－１
：キャノピー全体の気孔抵抗 ｓｍ－１
ｒｄ：キャノピー空間と地表面の抵抗 ｓｍ－１
ｒｇ：草全体としての気孔抵抗 ｓｍ－１
ｒsurf：裸地面の蒸発に対する抵抗 ｓｍ－１
：植生の脈管による面積平均の抵抗 ｓｍ－１
：土壌と根の面積平均の抵抗 ｓｍ－１
Ψl：葉の水分ポテンシャル ｍ
Ψr：根の地帯での土壌水分ポテンシャル ｍ
Ta：キャノピーsourceの高さの空気の温度 K
ea：キャノピ－空間の大気の水蒸気圧 ｍｂ
：Tc、Tgsでの飽和水蒸気圧 ｍｂ
Tｒ：参照レベルでの大気温度 K
eｒ：参照レベルでの水蒸気圧 ｍｂ
潜熱、顕熱フラックスを求める計算式は3.3章に示す。
2.2.3植生、土壌の構成
SiBで植生と土壌をどのようなモデルに構成するかは、水分の循環過程にとっては大事な要素である。
世界の植生は二つのグループ、キャノピーと地表面植生（草や他の草本性の植物）に分けられる（図3）。また計算の際に用いられるパラメータは、植生の種類や組み合わせから12パターンある。
図3：SiBでの植生の形態
土壌の構造は土壌温度を計算するときは1層で考えられる。地表面の植生と土壌を合わせて一つの温度としている。土壌水分を計算するときは3層に分けられる（図4）。土壌の表面約２cmが第１層、植生の根がある部分が第２層、根の届かない層が第３層となる。土壌に関するパラメータは土壌の種類によって6パターンある
図４：SiBでの根の状態と土壌の層
植生水分の詳細は3.6章で、土壌水分の詳細は3.7章で示す。
2.2.4 予報変数と診断変数
SiBモデルは7つの時間変化する値を得ることができる。それは次の2つの温度予報変数と５つの水分予報変数である。
Tc：キャノピーの温度 K
Tgs：地表面植生と土壌の温度 K
Mｃ：キャノピーの含水量 m
Mg：地表面植生の含水量 m
W1：飽和状態に対する第1層の土壌の水分含水率（無次元）
W2：飽和状態に対する第2層の土壌の水分含水率（無次元）
W3：飽和状態に対する第3層の土壌の水分含水率（無次元）
またキャノピーsourceの高さの空気の温度Ta とキャノピー空間の大気の水蒸気圧eの時間変化も得ることができる。この2つ変数を以下では診断変数と呼ぶ。
第3章 SiBの詳細
3.1 温度変化
3章のSiBの詳細はSellers and Mintz(1985)の論文を用いる。
予報変数である温度の時間変化ΔTc/Δt、ΔTgs/Δtは次式で表す。
(2) (1)
C：熱容量 Jm―２K－１
Rn：吸収された放射量 Wｍ－２
H：顕熱フラックス Wｍ－２
λE：潜熱フラックス Wｍ－２
T：温度 K
しかし、今回は時間ステップを短く考えている為、2式の右辺の温度によるエネルギー変化の項は微小なものと考えられるため、0とみなす。すると、
(3)
(4)
が得られる。この式は吸収された放射量から出て行った顕熱、潜熱を引いたものが温度変化になることを表す単純な熱エネルギー保存の式である。
3.2 放射フラックス
植生内の放射伝達式はDickinson(1983)の提案により次式である。
(5)
(6)
：入射フラックスに対する上向き、下向き拡散放射フラックスの割合
μ：直達入射光の天頂角の余弦
K：単位葉面積あたりの直達光の光学的深さ 
G(μ)：μ方向に投影した葉面の相対的面積
：散乱光の単位面積あたりの光学的深さ
：散乱フラックスの方向
β、β0：散乱光と直達光に対する上方散乱パラメータ
ω：葉の散乱係数
L：累積的な葉面積 ｍ２ｍ－２
ωβの値はNomanとJavis(1975)の解析より、
(7)
α：葉の反射係数
δ：葉の透過係数
：水平面からの葉の傾斜角
である。直達光に対する上方散乱パラメータはDickinson(1983)に従う。
(8)
散乱光の単位面積あたりの光学的深さ
は
(9)
である。
放射では高さ方向を放射フラックスが通過した葉面積で表す。つまり、キャノピー空間の上端ではL=0、キャノピー空間の下端から地表面植生に到達する前まではL=Ltc(キャノピーの葉面積(ｍ2ｍ－２))、地表面ではL=Lt(キャノピーと地表面植生の葉面積(ｍ2ｍ－２))となる。
そこで
を求めるための境界条件としてL=0でI↓=1（キャノピー空間に入る前だから）、L=LtでI↑=αdI↓が適用される(αd：拡散放射の地表面の反射率)。
キャノピー、地表面に吸収された放射量は
(10)
(11)
FΛ,μ(c,gs)：キャノピー、地表面に吸収されたFΛ,μ(0)の量 Wｍ－２
FΛ,μ(0)：入射放射エネルギー（Λ：波長 μ：放射角度） Wｍ－２
：キャノピーの上端を通過する拡散フラックス（L＝0のときのI↑）
：キャノピーの下端を通過する拡散フラックス（L＝LtcのときのI↓）
：キャノピーを貫通する直接flux
Ltc：キャノピーの葉面積 ｍ2ｍ－２ 
：面積平均のキャノピーの葉面積 ｍ2ｍ－２
AΛ,d、AΛ,b：拡散、直接フラックスに対する地表面スペクトル反射率
である。また地表面のスペクトル反射率AΛ,bは
(12)
αｓ(Λ) ：等方性と仮定した土壌スペクトル反射率
：草の上の上向き拡散フラックス（L=LtcのときのI↑）
である。
最後に吸収されたネットの放射フラックスは、吸収された放射量から放出された熱放射を引いて次式より得られる。
(13)
(14)
<Fc>、<Fgs>：キャノピーと地表面植生に吸収された放射の5つの合計 Wｍ－２
σｓ：ステファンボルツマン係数 5.67×10－８Wｍ－２K－４
Vc、Vg：キャノピー、草に覆われている割合
δt：熱赤外線のキャノピー透過率 
3.3 顕熱、潜熱フラックス
顕熱、潜熱フラックスは2.2章にあるように、電気回路に類似していると考えフラックス＝電位差÷抵抗で表す。電位差をここでは、顕熱の場合は温度差で、潜熱の場合は蒸気圧の差で決まる。（図2を参照）
・キャノピー植生からキャノピー内の空気へ
(15)
(16)
・地表面からキャノピー内の空気へ
(17)
(18)
γ：乾湿定数 
ｈs：土壌表面での空気の相対温度を調節する要素 
：土壌表面での空気の相対湿度 
Ψ1：土壌の最上位層の水分ポテンシャル 
ｇ：重力加速度 9.8ｍｓ－１
R：気体定数 287JkgK－１
ｃp：空気の定圧比熱 1004Jkg－１K－１
ρ：空気の密度 kgm－３
Wc、Wg：濡れているキャノピー、草の葉の割合 
Mc,、Mg：キャノピー植生、草の含水量 ｍ
Sc、Sg：Mc、Mgの最大値 ｍ
3.4 表面抵抗
キャノピー全体の気孔抵抗
は、キャノピーの葉それぞれの抵抗を並列と考え、さらに葉一枚の抵抗は全て等しいと考える。
個々の葉の気孔抵抗ｒsはJarvis(1976)が次式のように提案する。
(19)
Qp：光子flux密度 μ Ein m－１s－１
ｈ１：Qｐが無限大のときの最大コンダクタンス ｍｓ－１ 
ｒsmin：気孔抵抗の最小値 sm－1
：
（Qp＝0のとき）
、μEin m－２ｓ－１
ｒsmax：気孔抵抗の最大値 sm－1
κ：消去係数 
ωs：PARの散乱係数
Fs,μ(0)：キャノピーの上のPRA Fs,b(0)+ Fs,d(0)
：葉面温度T、葉の水分ポテンシャル
、水蒸気圧 
の
影響による調整要素 0～1
ここで、
(20)
Tc：キャノピーの葉の温度 K
：最適温度 (f(T)=1) K
：上限温度 (f(T)=0) K
：下限温度 (f(T)=0) K
(21)
ｈ5：種類依存定数 ｍｂ－１
：Taでの飽和水蒸気圧 ｍｂ
(22)
Ψ0：気孔が完全に閉じたときの葉の水分ポテンシャル ｍ
ｈ６：種類依存定数 ｍ－１
Ψlc：キャノピーの葉の水分ポテンシャル ｍ
(23)
Ψs：飽和時の土壌水分ポテンシャル m
B：経験定数
ｚd：根の深さ ｍ
Wi：i層の土壌水分 
：第i層の土壌中の水の体積割合 m3m－３
：飽和時の
の値 m3m－３
Di：i層の土壌の深さ ｍ
ha：キャノピーsourceの高さ ｍ
Ed：蒸発率 kgm－２s－１
ρｗ：水の密度 kgm－３
：植生の脈管による面積平均の抵抗 ｓ
(24)
R：根の単位長さあたりの抵抗 ｓｍ－１
Dd：根の密度 mm－３
Vr：土壌の単位体積あたりの根の体積 m３m－３
Kr：根の地帯での水のコンダクタンス ｍｓ－１
Ks：飽和状態の水のコンダクタンス ｍｓ－１
式(23)の右辺の第一項は土壌層の水分ポテンシャルの合計である。
そしてキャノピー全体の気孔抵抗
は個々の葉が並列に作用するため、
(25)
Nc：生きている葉からなるLtcの割合
：葉面角度分布関数（葉の方位ξ、傾角
）
となる。
地表面植生全体の気孔抵抗
もキャノピーのそれと同じ過程である。
裸地面の蒸発に対する抵抗ｒsurfはShu Fen Sun(1982)に従い、
(26)
ｄ１、ｄ２、ｄ３：経験定数
である。
3.5 空間抵抗
3.5.1 キャノピー空間とキャノピー葉面間の抵抗
高さとせん断力、風速の関係には次の関係式がある。
ｚ＞ｚ2にて
(27)
ｚ1＜ｚ＜ｚ2にて
(28)
あるいは
(29)
(30)
(31)
ｚ1＞ｚにて
(32)
：キャノピー空間の上限、下限の高さ ｍ
Kｍ：運動量輸送係数 ｍ2ｓ－１
u：風速 ms－１
Cd：葉の抵抗係数
τ：せん断力 kgm－１s－２
u*：摩擦速度 ms－１
ρ：空気の密度 kgm－３
ｄ：ゼロ面変位 ｍ
：粗度長 ｍ
：transition
height m
ｚ1からｚ2のせん断力については、式(28)と(29)のどちらを用いてもよい。そして式(27)から(32)から高さｚ1，ｚ2でのせん断力、風速について解き、β1からβ3を得る。
キャノピー空間とキャノピー葉面間の抵抗
は、次式となる。
(33)
C１：定数
u2：キャノピー空間の上限高さｚ2での風速 ｍｓ－１
ｌ：葉の長さ ｍ
ここで定数C１はβ1、β2、β3を使い、次の条件式より得る。
と
<0、あるいはβ１＜0のとき
(34)
と
>0、あるいはβ2＜0のとき
(35)
Cｓ：輸送係数
：面積平均の茎と葉の密度 m2ｍ－３
Ps：葉による遮断要素
式(33)の右辺の第一項は自由対流、第二項は強制対流を表している。つまり、自由対流と強制対流の影響を並列関係とみなしたものである。物理的に考えた場合に非現実的だが、風速が弱く放射量が大きい時にキャノピーと大気の温度差が非現実的になるのを防ぐような関数である。
3.5.2 キャノピー空間と地表面の抵抗
キャノピー空間と地表面の抵抗は運動量輸送係数Kｍを考慮して、
(36)
ｚgs：有効な地面の粗度長 ｍ
rd’：運動量輸送に対する抵抗 ｍ2ｓ－１
と表せる。高さの範囲はｚgsからhaだが、キャノピー空間の下端z1で分けて考える。まずｚgsからz1では次式で表す。
(37)
次にz1からhaでは次の条件式より得る。
あるいは β１＜0 のとき
(38)
あるいは β2＜0
のとき
(39)
式(36)から(39)によって得られたrd’を用いて、キャノピー空間と地表面の抵抗rdは次のように表す。
(40)
式(40)はGoudriaan(1977)、Businger(1971)によって使われた修正要素をもとにしている。ΦHで割ることによって、地表面とキャノピー空間の温度差が考慮される。
3.5.3 キャノピー空間と参照レベル間の抵抗
キャノピー空間と参照レベル間の抵抗はｒｄと同じように運動量輸送が考慮される。高さの範囲はhaからzrだが、キャノピー空間の上端z2とtransition height zmの２点で区切る。
(41)
haからz2は、同じキャノピー空間であるz1からhaの範囲で求められる式(38)、(39)と同じ方法で、範囲をhaからz2にする事によって得られる。
そしてz2からzm、zmからzrそれぞれの範囲で次の式より得る。
(42)
(43)
ここでΦ2はPaolson(1970)の修正要素で、顕熱フラックスによって変化する。
3.6 植生水分増加
SiBでは降水の遮断と排水については単純にモデル化されている。キャノピー、地表面植生による降水遮断の割合Pc、Pg(ｍｓ－１)とキャノピー、地表面植生からの排水の割合Dc、Dg(ｍｓ－１)は、
(44)
Mc<Sc のときDc=0
Mc=Sc のとき Dc=Pc
(45)
Mg<Sgのとき Dg=0
Mg= SgのときDg =Pg
P：降水率 ｍｓ－１
Sc、Sg：Mc、Mgの最大値 ｍ
K：消去係数
によって与えられる。またキャノピー、地表面植生の湿った部分の蒸発率Ewc、Ewg(kgm－２s－１)は、潜熱フラックスなので電気回路と考えることによって、
(46)
(47)
Wc、Wg：濡れているキャノピー、草の葉のfraction 
となる。よって植生水分増加ΔMc、ΔMgは遮断された水分による次の2つの支配方程式より得る。
(48)
(49)
3.7 土壌水分増加
2.2章のモデルの構造に記したように、SiBでは土壌温度は１層だが、土壌水分を考える場合３層に分けられる。まず土壌のi層からi+1層への水の移動Qi,i+1(ｍｓ－１)は、
(50)
によって与えられる。また重力排水のみによる底の層の流出Q3(ｍｓ－１)は、
(51)
ｘ：地表の傾斜角
となる。
土壌表面からの直接蒸発率Es(kgm－２s－１)は、
(52)
となる。またI層の根から吸い上げ蒸散した水の割合Edc,i、Edg,i(kgm－２s－１)は式(23)のキャノピーと地表面それぞれの葉の水分ポテンシャルを用いて次式より得る。
(53)
(54)
よって、３層の土壌水分の支配方程式は、
(55)
(56)
(57)
P１：最上位層の土壌への降水の潜入 ｍｓ－１
となり、土壌水分増加ΔW1、ΔW2、ΔW3が得られる。
3.8 診断変数
診断変数とは、ここではキャノピーsourceの高さの空気の温度Taとキャノピーの高さでの大気の水蒸気圧eaを表す。
キャノピー内の空気から大気境界層の参照高度への顕熱、潜熱スラックスを考えると次式になる。
(58)
(59)
式(58)、(59)の左辺の顕熱、潜熱フラックスと右辺のキャノピー空間と参照レベル間の抵抗を計算した上で診断変数が得られる。
3.9 雪の影響
SiBでは参照高度の温度Trが融点Tfを下回ったら降水が降雪になると定義する。雪による影響はいくつかある。
地表面植生に遮断される雪の深さzfは、
(60)
Mg：地表面植生の含水量 m
によって計算される。つまり雪の場合は水の5倍と仮定される。
放射フラックスに対する散乱係数は遮断された雪の量によって調整される。
(61)
：キャノピーと遮断された雪による散乱係数
ω：葉の散乱係数
ωf：雪の散乱係数
地表面の反射率は散乱係数と同じように調整される。
(62)
(63)
：雪のある地表面のアルベド
Ags：雪のない地表面でのアルベド
Af：雪のアルベド
雪は放射フラックスに対しては反射のみで透過はしないと仮定する。よってωf＝Afとなる。それらの値は可視光に対しては0.8、近赤外域に対しては0.4である。
ゼロ面変位d、粗度長z0、抵抗
の計算過程で用いる定数C1、キャノピーの遮断容量の最大値Scは次のように調整できる。
(64)
(65)
(66)
(67)
(68)
‘～’は調整値を表す。
エネルギーでは潜熱フラックスが次式のように調整される。
(69)
λf：結晶の潜熱 JKg－１
第4章 SiBモデルの実行
4.1 計算の流れ
前章のSiBモデルの詳細を使って計算ができるかどうかを試みる。そこである1日分の予報変数、診断変数の変化を計算してみる。計算に使うものとして、今回は使用するのがもっとも簡易なExcelを用いて計算する。
計算の流れは、まず植生が吸収する放射量を求める。そして植生の気孔抵抗、キャノピー、地表面とキャノピー内の大気の間の空間抵抗を求め、植生が費やす顕熱フラックス、潜熱フラックスが求められる。吸収された放射量、費やされた顕熱、潜熱フラックスにより、その時間ステップでの植生の温度変化を得る。
植生の水分変化は、植生の葉が遮断した降水の割合、排水した降水の割合、そして植生の湿った部分からの蒸発率から得る。
土壌層の水分変化は、各層間の水分移動、土壌表面からの直接蒸発率、土壌内から植生によって抽出された水分から得る。
そして７つの予報変数の変化量が求められ、ここまでで最初の時間ステップは終了し、次の時間ステップで変化した値を加えた予報変数を使い同じ計算を進めていく。また次の時間ステップから、前の時間ステップで求めた顕熱、潜熱フラックスを用いてその次の診断変数を得る。
計算の流れの概略図を図5で示す。
予報変数 Tc、Tgs、Mｃ、Mg、W1、W2、W3
放射の吸収 Rnｃ、Rngs (5)～(14)
表面抵抗 ｒｃ、ｒｇ、ｒsurf (19)～(26)
次の時間ステップへ
空間抵抗 ｒb、ｒd (27)～(43)
潜熱、顕熱 λEｃ、λEgs、Hｃ、Hgs (15)～(18)
診断変数 Ta、e a (58)～(59)
温度変化 ΔTｃ、ΔTgs (3)～(4)
水分増加 ΔMｃ、ΔMgs、ΔW１、ΔW2、ΔW3 (44)～(57)
図５ SiBの計算の概略図
4.2 観測値、パラメータ
SiBモデルで使用する観測値は岐阜県高山市にある高山測候所で、2004年8月1日の正午から2004年8月2日の正午まで観測された現地気圧、気温、相対湿度、風速を用いる。
パラメータはSellers et al(1989,1995)とWalter Larcher(2004)を用いる。植生に関するパラメータは地表面を覆う植生で12パターンに、土壌に関するパラメータは土壌の種類から6パターンある。そこで高山の場合、植生はスギの木を考えて常緑針葉樹、土壌は有機土壌のパラメータを用いる。またスギの木の場合地表面植生はほとんど生えていないため、全く生えていないとみなす。
4.3 仮定した数値
SiBモデルを計算するときの必要な値として分からない値が数多くあった。そこでそれらの値を存在しうるように仮定した。
キャノピーの熱容量Ccは、1ｍ2あたりのキャノピーの体積を考える。観測地である高山のスギの木の状態から、直径0.4ｍ、高さ12ｍのスギの木が５m間隔に１本生えていると仮定する。すると体積が分かり、スギの木の比重300kg/m3と比熱1.3kJ/kgKを用いると、熱容量は0.023Jm－２K－１となる。
予報変数は計算する事によって次の時間ステップの値を得るが、初めの時間ステップでは値が必要である。よって、キャノピーの温度Tcは正午ということで気温＋５℃で、キャノピーの含水量Mcはキャノピーの体積の５割と考えて0.03ｍで、i層の飽和状態に対する土壌水分の割合Wiは３層とも0.8と仮定する。
入射放射エネルギーFΛ,μ(0)の直達光、拡散光の分け方は拡散光が全体の1/4と仮定する。また、可視光、近赤外域の分け方は4.4章で計算によって得られる。
単位葉面積あたりの直達光の光学的深さKはMonsi and Saeki(1953)、Kira et al(1969)の葉面積指数に対する放射減衰率のグラフから、葉面積指数が1のときの放射減衰率は50％と仮定する。すると放射の式(5)と(6)の右辺のe―KLは減衰率を表しているので、e―K×１＝0.5からK＝0.7となる。
生きている葉からなるLtcの割合 Ncは全てが生きていると仮定して1とする。
散乱光の単位面積あたりの光学的深さ
と散乱光と直達光に対する上方散乱パラメータβ0はどのような範囲の値かもわからない為、１と仮定して計算する。
4.4 計算が必要な値
３章で述べたSiBモデルから計算するのだが、その中にある値で、３章に述べていない計算によって得る値がある。その方法を以下で述べる。
温度Tでの飽和水蒸気圧
はTetens(1930)の式から、
(70)
より求める。また温度Tの大気の水蒸気圧eは、
に相対湿度を掛けて求める。
空気の密度ρは理想気体の状態方程式から気圧ｐ、大気の気体定数Rとして、
(71)
である。
i層の土壌の深さDiはSellsrs et al(1995)から、３層合計の土壌の深さDTを2ｍ、根の深さDrが1.5ｍより、D1は0.02ｍ、D2はDr－D1＝1.48ｍ、D3はDT－Dr＝0.5ｍとなる。
輸送係数Cｓ、葉による遮断要素Ps、transition height
はSellers et
al(1995)より次のように得る。
(72)
(73)
(74)
lw：葉の幅の長さ ｍ
lwが0.001ｍ、ｚ2が17ｍ、ｚ0が0.05ｍ、
が0.94ｍ2ｍ－２よりCｓは9、zmは20.59ｍ、Psは1.96となる。
土壌の単位体積あたりの根の体積Vrは、まずMi-sun Lee(2003)より1ｍ2あたりの根の炭素Cの量が673ｇである。根の成分はセルロースCH2Oなので炭素が673ｇのときセルロースは1682.5ｇである。そして比重400kgと根の深さ1.5ｍで割れば、Vrが2.8×10－３ｍ３ｍ－３となる。
根の密度Ddは単位体積あたりの根の長さのことである。根の半径をｒとすると体積は体積Vrはπ×ｒ2×Dd、根面積指数は2×π×ｒ×Ddであるが、Vrは2.8×10－３ｍ３ｍ－３で根面積指数は10ｍ2ｍ－２なので、根の密度Ddは2842ｍｍ－３となる。
根の単位長さあたりの抵抗RはWalter Larcher(2004)から比水分通導度LSCが針葉樹の場合4×10－４ｍ２ｓ－１Mpa－１より、単位を直して逆数にするとRは2.5×10５ｓｍ－１になる。
植生の脈管による面積平均の抵抗
に関しては、まずWalter
Larcher(2004)より流量Q、圧力勾配dp/dzとすると比水分通導度LSCは、
(75)
と表すことができる。そこで式(22)の右辺の第３項の（Ed/ρｗ）×
は水頭差である。圧力勾配を用いても水頭差は表せるので、
(76)
となる。（Ed/ρｗ）は流量なので式(75)(76)から、
(77)
となる。よって、
は4×105ｓとなる。
入射放射エネルギーFΛ,μ(0)の可視光、近赤外域は近藤純正(1994)より得る。まず大気上端の水平面に入射する可視光の量S0↓は、
(78)
Φ1：緯度（北は+、南は－）
ｈ：太陽の南中からの時角
M：月数 DAY：各月の日にち
である。それを用いて地上における可視光の量S↓は、
(79)
ｐ：地上の気圧 ｍｂ
e：水蒸気圧 ｍｂ
ref：地表面の反射率
より得る。近赤外域の放射量の日平均値Lf↓は、
(80)
σ：ステファンボルツマン定数 5.67×10－８Wｍ－２K－４
T1：日平均気温 K
より得る。なお式(79)(80)は推定式である。
4.5 展開が必要な式
３章でのSiBモデルの詳細はSellers and Mintz(1985)をもとにした内容であるが、SiBモデルをExcelで計算するときに３章の式のままでは計算できない式がいくつかある。それらをExcelで計算できるようにした式と導いた方法を以下で述べる。
4.5.1 放射フラックス
式(5)と(6)の植生内の放射伝達方程式は、上向き放射フラックスI↑、下向き放射フラックス↓を求める事ができるように変形する必要がある。
放射伝達方程式には微分形が含まれているので、式(5)(6)を葉面積Lで積分する。しかし不定積分なので、積分定数D1、D2が発生する。積分定数を解く為に、3.1章で述べたように2つの境界条件L=0でI↓=1、L=LtでI↑=αdI↓を用いて積分定数を与え、I↑とI↓を求める事ができる式を導く。その導いた式は、
(81)
(82)
である。そして、式(81)の葉面積にキャノピーの葉面積Ltcを代入すると式(10)で用いるキャノピーの下端を通過する拡散フラックス
に、式(82)の葉面積に0を代入すると式(10)で用いるキャノピーの上端を通過する拡散フラックス
に、式(82)の葉面積にキャノピーの葉面積Ltcを代入すると式(12)で用いる草の上の上向き拡散フラックス
になる。
4.5.2 気孔抵抗
今回計算する気孔抵抗は、草はないものとして考えているのでキャノピーの気孔抵抗
のみを求める。そのために式(24)を用いるが、葉の向いている方位や角度を表す葉面角度分布関数
を得る事が出来なかった。しかし、葉はある程度太陽の方向を向いているはずである。そこで葉は天頂方向から90°を向いていると仮定する。すると式(24)は、
(83)
と表すことができる。ｒsは式(19)からわかるので、ｒsの逆数を積分範囲0からLtcで積分すると、
(84)
となる。
4.5.3 空間抵抗
3.5章にある3つの空間抵抗を得るには、式(29)の中にあるパラメータβ1、β2、β3の値を得る必要がある。そのために式(26)から(32)の、キャノピー空間の下端z1と上端ｚ2で区切ってｚ2以上、ｚ1からｚ2、ｚ1以下それぞれの風速、せん断力を表す式を用いる。
まずz1からｚ2の風速は式(29)だが、右辺に風速u2が含まれている。よって、式(29)で高さｚ2での風速を表すとき次に条件を導くことができる。
(85)
次にz1とｚ2それぞれでのせん断力τは同じであるという条件を考える。z1では、ｚ1からｚ2のτを表す式(27)の右辺と、ｚ1以下のτを表す式(32)の右辺を高さで微分した式は等しい。つまり、
(86)
と表せる。ｚ2は、ｚ1からｚ2のτを表す式(28)の右辺と、ｚ2以上のτを表す式(26)の右辺は等しい。また、ｚ2以上のせん断力は空気の密度と摩擦速度で表すことができるのでキャノピー上端より上空はせん断力は一定である。よって式(26)の高さと風速は、参照レベルの高さzrとその高さでの風速urを用いて、
(87)
と表せる。
式(85)、(86)、(87)の3式からβ1、β2、β3を導く。その結果、
(88)
(89)
(90)
となる。未知数3つに対して式(88)(89)(90)の非線形の式が3つなのでSellers and Mintz(1985)に従い、β3に関する繰り返し計算を行う。初めにβ3を適当に与えてβ2、β1と求め式(90)のβ3を求めてそのβ3を使ってまたβ2、β1を求める。これを繰り返し、収束するまで行いβ1、β2、β3の値を得る。
4.6 計算方法および結果
今回はExcelで計算を行う。その方法として、行を一つの時間ステップにして予報変数が求められたら、次の行で前の時間ステップで得た予報変数を用いて次の時間ステップを計算する。しかし、4.3.5章で述べたβ1、β2、β3の計算や、式(34)、(35)、(38)、(39)にあるnが1から無限大まで合計する計算は収束するまで繰り返し計算が必要なので、1行で1時間ステップではできず、β1、β2、β3の計算は10行で1時間ステップ、式(34)、(35)、(38)、(39)にある計算は3回で収束するようだったので3行で1時間ステップとした。
計算過程は以上の方法で、進めていけば予報変数の値が求められ、その値を次に時間ステップで使って新しい予報変数が求められるはずである。しかし、4.5.3章のβ1、β2、β3の値が収束すると思われたが、一定の値に定まらなかった。そこで式(88)、(89)、(90)からさらに違う方法を試みた。
まず式(89)、(90)からβ1を消しβ2、β3のみの式となる。その式と式(88)で未知数β2、β3による2つの非線形の式を得る。そこで線形の式にするため、テーラー展開を使う。式(88)を右辺が0の形にしてその式をF(β2、β3)＝0とおく。また式(89)、(90)から得た式を同じようにG(β2、β3)＝0とおく。その式をテーラー展開より、
(91)
(92)
となる。ここで初期値β2、β3を与えると式(91)、(92)はβ2、β3の変化分であるδβ2、δβ3のみとなり、行列で表すと、
(93)
よりδβ2、δβ3の値がわかり、それを初期値に足し変化分を求める。これを繰り返し、変化分が0に限りなく近づいたところのβ2、β3が正しい値である。しかし、この方法も変化分が0に近づかないという結果になった。
最後に正確な値を得ることはできないが、テーラー展開をするときに得たF(β2、β3)＝0 とG(β2、β3)＝0をグラフに描いて、その交点を調べた。すると2本のグラフは近づいていくが交わらないという結果になった（図6）。
図6：F(β2、β3)＝0 とG(β2、β3)＝0のグラフ
β1、β2、β3の値を得られなかったため、空間抵抗が得られず顕熱、潜熱フラックス、予報変数のキャノピーの温度変化が得られなかった。そして、次の時間ステップから先では当然、変化した分を足したキャノピーの温度は必要なので計算できなかった。
そこで、キャノピーに吸収された放射量Rnｃを計算するためにキャノピーの温度Tｃ、地表面の温度Tgsの値が必要であるが、計算する過程において式(13)、(14)で用いるだけで、大きな影響を与える値ではないと考えられる。よってTｃ、Tgsを観測値である観測値と同じとし、またTｃ、Tgsは大きく変化する値ではないので時間変化はしないと仮定する。そして2004年8月1日の正午から2004年8月2日の正午までのRnｃを計算した結果が図7である。
(W/m2)
図7：キャノピーに吸収された放射量の時間変化のグラフ
4.7 考察
計算内容に関しては、収束した値を得る計算や式(34)、(35)、(38)、(39)にあるnが1から無限大まで合計し収束するまで繰り返し行う計算等はかなり複雑で難しい内容となっている上に、Excelで行う場合は毎時間ステップで繰り返し計算を行う作業が必要である。またExcelにてSiBモデルを実装するために、式を直したり必要なパラメータを他の資料で調べて値を得る等の作業も必要である。
計算結果では空間抵抗を求めるために必要なパラメータβ1、β2、β3の値を得るために幾つかの方法で試みたが、結局得られなかったためキャノピーの温度変化、キャノピーの含水量、土壌水分の含水率の時間変化を得られなかった。しかし得られなかったキャノピー、地表面の温度の時間変化を仮定してキャノピーに吸収された放射量の時間変化を得ることができ、正午で高くなり太陽が沈んだ夜の時間帯で熱を放出する為、負になるという現実的な値を得ることができた。
β1、β2、β3の値を得ることができなかったのはExcelが原因で得られなかった値ではないので、そのパラメータさえ求めることが出来れば、空間抵抗を求めてキャノピーの温度の時間変化分を得ることが出来る。そして予報変数の時間変化を計算していくことはでき、ExcelにてSiBモデルを実装できると考えられる。
第5章 まとめ
本研究では陸面過程モデルの中心であるSiBが放射フラックス、温度変化、植生、土壌の水分変化についてどのように構成されているかを詳細に調べた。そしてSiBモデルを使って実際計算できる状態にするために改良し、使用するのが最も簡易なExcelを用いて計算を行った。
SiBの詳細では、地球規模に対応させるために計算量を減らすよう簡略化している部分がいくつかあることがわかる。例えば、顕熱フラックスを求めるときに植生とキャノピー空間の温度差と、潜熱フラックスを求めるときに植生とキャノピー空間の水蒸気圧の差を一定と捉えているのがわかる。しかし、葉と大気の間の温度あるいは水蒸気圧の差は葉一つ一つ違うはずであるが、それを考慮して地球規模で計算すると膨大な計算になるため簡略化されているのがわかる。
またSiBを用いて計算するには第4章全体で述べたように、値を得るためにさらに計算が必要な所や、SiBモデルとして示されている式では計算が困難な式を実際計算できる式に展開する必要がある所などが存在することがわかった。
そして実際Excelで計算した結果は計算することができなかったパラメータがある為予報変数の時間変化を計算することができなかったが、キャノピーに吸収された放射量の1日の時間変化を得ることができたので、得られなかったパラメータさえ求められれば全体の計算ができるので、ExcelでSiBを実装できると考えられる。
今後の課題としては4.5.3章で述べた空間抵抗を求めるためのパラメータの値を得る事ができない原因を調べ、SiBモデルをExcelによって計算する。そして、一番新しい第3世代のモデルは今回調べたSiBモデルと同じように地球規模に対応させるために作成されている。よってそれらを詳細に調べれば、狭い範囲で考えた場合にSiBと同じように計算量を減らすために幾つか簡略化していると思われる。よってその部分を改良し、改良したモデルが狭い範囲で再現できるかを検討していきたい。
謝辞
本研究を進めるにあたり、ご多忙であるにもかかわらず、終始懇切丁寧な助言、ご指導をしていただきました、玉川一郎助教授には心から感謝の意を表します。
また、当研究室の学部生、院生の方々にも数々のご協力、ご支援を頂きました。
この場を借りてお礼申し上げます。
参考文献
高山測候所wwwページ：http://www.tokyo-jma.go.jp/home/gifu/takayama/
近藤純正：水環境の気象学 －地表面の水収支・熱収支－、朝倉書店、pp55-92, 1994
Kira, T., Shinozaki, K., Hozumi, K.: Structure of forest canopies as related to their primary productivity. Plant Cell Physiol. 10, pp129-142, 1969
Mi-sun Lee, Kaneyasu Nakane, Takayuki Nakatsubo, Hiroshi Koizumi : Seasonal changes in the contribution of root respiration to total soil respiration in a cool-temperate deciduous forest, Plant and soil, 255, pp311-318, 2003
P. J.Sellers, Y. Mintz : A simple biosphere model (SiB) for use within general circulation model , Journal of the atmospheric sciences, 43, pp505-531, 1985
P. J. Sellers et al : A revised land surface parameterization (SiB2) for atmospheric GCMs. PartⅠ: model formulation , Jaurnal of climate , 9, pp676-705, 1995
P. J. Sellers et al : A revised land surface parameterization (SiB2) for atmospheric GCMs. PartⅡ: The generation of global fields of terrestrial biophysical parameters from satellite date, Journal of climate , 9, pp706-730, 1995