陸水物理学会誌 第 3 巻 第 1 号

2021 年 12 月1 日発行

ISSN 2435-3043

陸水物理学会誌
Journal of the Japanese Society of Physical Hydrology

第 3 巻 第 1号

2021 年12月

陸水物理学会
The Japanese Society of Physical Hydrology

陸水物理学会誌
第3巻 第1 号

目2021 年 12 月

次

巻頭言

浩

...

1

樹液流測定による流域蒸散量の評価：3次元数値実験との比較
鈴木健史，知北和久

...

3

...

21

諏訪

論 文

論 文
亜寒帯湿原における小湖沼の熱的変動 -浮島の浮沈との関係岡田操

会務報告 .................................................................. 41
投稿規程 .................................................................. 43
執筆要領 .................................................................. 45
査読規程 .................................................................. 49
編集後記 .................................................................. 52

表紙画の説明：
青森県・青池の写真。陸水物理学会誌第1巻及び第2巻に掲載された論文著者である花石竜治会員撮影。

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 1(2021)

巻頭言 (Foreword)

天変地異の恵み浩1

諏訪
暴風や豪雨，落雷などは天変と称し，地震や火山噴火，山崩れなどの地異と合わせて天変地異と云う。
“天
変地異は禍をもたらすこともあるが，恵みをももたらしている”と云うと，“何をバカな”と云われるかも知れない。我々は地震や台風で物的・人的被害を受けることがある。災害である。それで被害を避けるあるいは軽減するための対策に取り組んできた。
そもそも天変地異と称する現象自体は災害ではない。
地球の営みである。天変地異を含む地球の営みによって生成・維持される大気・水・土地・生物からなる多様な環境の恵みを我々は享受している，とも云える。
ヒトはこのような環境に適応しつつ進化してきたからでもある。
では，例えば台風はどのような恵みをもたらしているか。地球には太陽からの放射が降りそそぐ。単位面積当たりにもたらされるエネルギーは緯度による偏りが大きい。ところが，この偏りは大気と海洋熱循環の原因であり条件である。そして，この循環はエネルギ
ー配分の偏りを小さくするべく働く。台風はこの熱循環の一端を担う。すなわち，低緯度海域で大量の熱エネルギーを仕入れ，これを高緯度地域へ配り歩く。大量の水を水蒸気の形で運び，雨を降らせると同時に，放出される水蒸気の潜熱で周囲の大気を加熱する。水と熱の宅配便である。台風は毎年のように日本列島に上陸して，あるいは 1000 km 以遠からでさえ水蒸気を送り込み，列島に雨をもたらす。このため我が国土が旱魃に見舞われることは多くない。瑞穂の国の繁栄は
台風によっても支えられてきた。
列島は中緯度高圧帯，すなわち降水量が少ない場所に位置している。この緯度の陸域における年降水量の平均は 0.6 m ほどであるが，日本列島では 1.7 m にも達する。海に囲まれていることに加えて，山がち，すな
わち脊梁山地が発達していて雨を呼ぶ地形に富むためでもある。4 つのプレートがぶつかり合う場所にあるので，地殻変動が静まることは無く，火山活動が絶えることも無い。地震や火山噴火は被害をもたらすが，陸地の起伏を生成・維持して豊かな降水を保証してくれている。
そのような営みのメッカの一つが信濃川を遡り，源流へ向かう地にある。岐阜県と長野県の県境に焼岳がある。1915 年（大正 4 年）噴火で泥流が発生して梓川を堰き止めた。そのため新湖ができて，大正池と呼ばれることになった。1928 年には湖水を利用して水力発電が始まる。
梓川発電所と称した発電所は，
その後1951
年に東京電力霞沢発電所となり，現在に至る。写真左下隅に写るのは霞沢発電所の取水口付近である。大正池は焼岳や槍穂高の山々を湖面に映し，山岳観光地“上高地”景観の立役者を担う。
上高地がある飛騨山脈とその周辺では地震や火山噴火，豪雨が起きて様々な災害が繰り返されている[1]。
しかしこの繰り返しが飛騨山地を我々が暮らし，観光を楽しみ，生物の多様な生存を支える場所にしている。
地すべりや土石流だけでなく，地震や火山噴火などのハザードの素因や誘因に水の存在とその変化が強く関与することも知られるところとなっている。

霞沢発電所取水口そばから望む焼岳。手前は大正池の湖面
参考文献
1. Suwa H (2020): A review of landslide and debris flow hazards
in the Hida Mountains, Japan. Transactions, Japanese
Geomorphological Union 41(2): 127-146.

*1 東京大学空間情報科学研究センターおよび立命館大学歴史都市防災研究所客員研究員，もと京都大学防災研究所
DOI: https://doi.org/10.34502/physhydro.3.1_1

1

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

原著 (Research Article)

樹液流測定による流域蒸散量の評価：3 次元数値実験との比較
鈴木健史 1・知北和久 2

Estimating transpiration from a catchment by sap flow measurements: Comparison with 3D numerical simulation
Takeshi SUZUKI 1 and Kazuhisa A. CHIKITA 2
摘要 本研究では，樹幹の樹液流測定法の一つであるヒートパルス法を用い，流域にある代表的な樹木 50 本を選びその樹液流速を測定した。また，ヒートパルス法の測定原理に対する 3 次元数値実験による補正に基づいて，流域蒸散量を評価した。対象流域は，札幌市の豊平川上流にある冷水沢支流流域（1.16 km2）で，流域面積の 56%は針葉樹のトド
マツ林，他は広葉樹との混交林である。ヒートパルス法による観測期間は，2001 年 7 月 9 日～9 月 20 日である。先ず，ヒートパルス法によって代表的な針葉樹木と広葉樹木の深さ方向の樹液流の流速分布を求め，樹木の辺材部で積分して各樹木の蒸散量 Tr1 (L day-1)とした。この Tr1 値の無降雨日の平均値 Trm は，樹木の胸高直径 DBH(cm)をパラメーター
として Trm=c(DBH)d (c, d は樹種ごとの定数)の関係があることがわかった。この関係式に対する決定係数 R2 は，トドマ
ツで 0.93，広葉樹としてイタヤカエデ，シナノキ，オヒョウなどをまとめた場合 0.71 と求められた。ヒートパルス法
によって，観測期間中の無降雨日での流域全体の平均蒸散量は 0.72 mm day-1 と求まった。他方，ペンマン-モンティース法による蒸発散量は 2.03 mm day-1 と見積もられた。林床からの蒸発散量を考慮しても，前者の値は後者よりかなり小さく，その理由として，センサーの直径が 2 mm であるにも関わらず，樹液流の速度を“センサーの体積が無限に小さい”という仮定の基で求めたためと考えられた。そこで，辺材部での樹液流の阻害パターンとしてセンサーを挿入した時の辺材部のキズの幅を考慮し，樹液流に関する 3 次元数値実験を行った。その結果，キズの幅 2.8 mm のとき流域の平均蒸散量 1.76 mm day-1 と求まり，クマイザサによって覆われている林床からの蒸発散量を考えると比較的妥当な値が得られた。
キーワード ヒートパルス法，樹液流，蒸散量，森林流域，ペンマン-モンティース法
Abstract Sap flow velocity was measured by the heat-pulse method in a forested catchment (1.16 km2 in area), andtranspiration from the catchment was estimated by three-dimensionally simulating actual sap flow blocked by sensors of definite
size. The catchment is mostly forested by 56 % Abies sachalinensis (Sakhalin fir) and 44 % mixture with broad-leaved trees. Radial
distributions of sap flow in the sapwood part for representative fir and broad-leaved trees were obtained by the heat-pulse method, andthen transpiration Tr1 (L day-1) of each tree was calculated by integrating the sapwood area. The Tr1 values averaged for non-rainfalldays proved to be highly correlated to trees’ breast high diameters DBH (cm) with R2=0.93 for Sakhalin fir and R2= 0.71 for broadleaved trees. The catchment transpiration averaged over non-rainfall days was 0.72 mm day-1 for the heat-pulse method and thecatchment evapotranspiration, 2.03 mm day-1 for the Penman-Monteith method. The possible underestimate of transpiration from theheat-pulse method is probably due to the ignorance of sensors’ size originally supposed in the heat-pulse method. Actual sap flow
blocked by sensors of 2.8 mm in diameter was three-dimensionally simulated. Giving the cut size of 2.8 mm by the sensors in thesimulation, the catchment transpiration was estimated at 1.76 mm day-1. This value is likely reasonable to the evapotranspiration at 2.03mm day-1 from the Penman-Monteith method, since the evapotranspiration includes the evapotranspiration from forest floor coveredby striped bamboo (Sasa senanensis).
Keywords Heat-Pulse Method, Sap Flow, Transpiration, Forest Catchment, Pemnan-Monteith Method*1 (株)日立ソリューションズ 〒140-0002 東京都品川区東品川 4-12-7 Hitachi Solutions, Ltd., Tokyo, 140-0002*2 北海道大学北極域研究センター 〒001-0021 札幌市北区北 21 条西 11 丁目 Arctic Research Center, Hokkaido University, Sapporo, 001-0021
DOI: 10.34502/physhydro.3.1_3

3陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

1. はじめに森林流域における蒸発散は，陸域の水収支・熱収支を考える上で重要な過程であり，水文学・気象学・生態学など幅広い分野で扱われている。このため，蒸発散量とその成分（蒸散量・遮断蒸発量・林床からの蒸発散量など）について，多くの観測方法が利用されてきた。その代表的な例は，流域水収支法[1]，渦相関法[2]，ボーエン比法[3]，ペンマン-モンティース法[4]，樹液流測定法[5-13]である。これらの方法には，それぞれの基本原理に基づく利点と欠点がある。例えば，流域
水収支法では大きな流域スケールでの蒸発散量が評価できるが，時間スケールは 1 年単位である。また，渦相関法やボーエン比法では30 分～1 時間ごとの測定も可能であるが，測定可能な領域には制約がある[2,3]。
現在最も多く使われているのが，流域熱収支に基づ
くペンマン-モンティース法である。その理由は，流域を代表する 2 つのパラメーターである空気力学的抵抗と群落抵抗の決定法についての検討が進み，これらの値を得るのが比較的容易になった点にある。
樹液流測定法は，単木，林分および流域の蒸散量を求める手法として，今日まで多くの改良型が考案されている[10, 13, 14]。この手法の中で幹熱収支法[8]，ヒートパルス法[5]，グラニエ法[11]が知られており，このうちグラニエ法はプローブ製作が容易で安価なこと[15]，および測定・校正の容易さから，わが国で最も広く用いられている[16, 17]。但し，ヒートパルス法とグラニエ法では，樹木にプローブを挿入することに因って辺材部に損傷を与えることになり，これによる樹液流への阻害や長期観測による樹木の損傷反応を常に考慮する必要がある[18, 19]。なお，幹熱収支法は樹木に損傷を与えないが，与えた熱に対する応答遅延で樹液流速に対する誤差が生じやすく，また測定できる樹木の大きさに限界がある[20]。
本研究では，樹液流測定法の中で最も早期に開発され，現在も利用されているヒートパルス法を採用した[14, 21]。この手法は複雑な地形や立木密度の大小に関係なく蒸散量のみを測定できるが，現在，トウモロコ
シやヒマワリといった細い茎内の樹液流の速度を求める方法としても使われている[14]。他方，ヒートパルス法で得た樹液流の速度を用いると流域蒸散量がかなり過小評価されることが以前から指摘されている[18,22-25]。ここでは，この方法で得た樹液流速度の値を数値実験によって補正し，これによって得た流域蒸散量とペンマン-モンティース法による流域蒸発散量とを比較検討した。

図 1. 冷水沢支流域と観測点の位置．斜線域はトドマツ
林．実線上の数値は標高(m)を表す．T 地点：気象観測点，
S1, S2 地点：樹液流速度の観測点，黒点：毎木調査地点
Fig. 1. Location of a catchment for a tributary of the Hiyamizu
Stream, and observation sites in the catchment. Shaded and whiteareas are occupied by Sakhalin fir and its mixture with broadleaved trees, respectively. Numerical values on the contours
show the elevation (m above sea level). Meteorology at site T,
sap flow measurements at sites S1 and S2, tree survey at blackplots

2. 観測流域
観測流域は，北海道・札幌市郊外の豊平川上流にある冷水沢支流域（流域面積1.16 km2；北緯4255
～425608.78，東経1410936.86～1411053.55）である（Fig. 1）。本流域の南約3 kmにある豊平峡ダム管
理所の気象データによると，著者らが樹液流を観測した2001年では，年平均気温4.7 ℃，年降水量1,196 mmである。流域の上流部と中流部右岸には針葉樹のトド
マツ林（0.63 km2）が分布している｡その他は混交林（0.48 km2）で，針葉樹のトドマツ(Sakhalin fir；学名
Abies sachalinensis)，広葉樹のイタヤカエデ(painted
maple；学名Acer pictim)，シナノキ(Japanese linden； 学
名Tilia japonica)，オヒョウ(Manchurian elm；学名Ulmus
laciniata)，ウダイカンバ（monarch birch；学名Betula
maximowicziana）などから成る。トドマツ林の大部分は1957～1958年に植林されたもので，これらは本研究の観測年2001年では樹齢43～44年の一斉林といえる。流域内にある林道とFig. 1のT地点のある草地（合計面積0.05 km2）以外は，全体として密に植生に覆われている。
なお，林床の植生は，トドマツ林，混交林ともに主にクマイザサ（striped bamboo; 学名 Sasa senanensis）である。このため，クマイザサによる遮断蒸発と蒸散，
4陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

地面付近の土壌水の蒸発が起こりうる。

3. 方法
3.1. ヒートパルス法の原理
ヒートパルス法では，熱をトレーサーとして樹幹内辺材部の樹液流の速度を測定する（Fig. 2）[5, 6]。Fig.
2a のように直径 2 mm のヒーターと温度センサーを樹木に垂直に挿入する。Fig. 2b は，その正面図で，樹木の垂直方向を x 軸，接線方向を y 軸とする。ここで，深さ方向 z については樹液流の速度 V と樹幹の熱特性が一様であるとして，xy 平面のみを考える。Fig. 2b のようにヒーターを原点とすると，上流側 xup，下流側
xdown での温度センサーの座標は，単位を cm で表すとそれぞれ(-0.5,0)，(1.0,0)となる。ヒーターからヒートパ
ルスが放出されると，その熱は樹幹内で拡散する。これと同時に，熱は樹幹内の樹液流によって移流拡散する。
Marshall [5]によると，
ヒートパルスの移動速度 HPVと樹液流の速度 V との関係は次の(1)式で表せられる：
𝐻𝑃𝑉

𝜌𝑤 𝑐𝑤
V
𝜌𝑐

図 2. (a) ヒートパルス法におけるヒーターと温度センサ
ーの固定位置と(b)座標系．心材の中心に向かう方向を z
軸（深さ方向）とする．
Fig. 2. (a) Location of a heater and two temperature sensors in atree’s trunk and (b) a coordinate system for the heat-pulse method.
The z axis is set to the center of the core.

(1)

𝐻𝑃𝑉

ここで，，c：それぞれ樹幹の密度と比熱，ρw，cw：それぞれ水の密度と比熱である。
他方，Marshall [5]はヒーターからの放熱量 Q，経過
時間 t，温度上昇 T とし，Fig. 2 の xy 平面に対する 2 次
元の熱移流拡散方程式を，下記の 4 つの仮定のもとで解析的に解くことを考えた：となる。これより，上下センサーの温度が等しくなったとき，与えられた熱はちょうどその中間地点まで移動したと考えられる。つまり，(3)式は，その時間と距離の関係から速度を求めている。本研究で利用した樹液流の流速計は，この t0 の測定から HPV を計算している。最終的には，この HPV 値から(1)式を用いて樹液流の速度 V を求めた。このとき，樹幹内の密度と比熱は，口径 4 mm の成長錐でコアを採取し，その温度，
体積，重量および含水率から求めた。

・熱量 Q は x = y = 0, t = 0 で瞬間的に与えられる・ヒーターとセンサーの体積は無限小である・樹液流の速度 V は xy 平面上で一様に分布する・樹木は熱的に等方均質である(kx = ky)

3.2. 観測方法
観測は，2001 年 7 月～10 月に実施した。Fig. 1 の林
内 S1 地点でイタヤカエデ，シナノキ，オヒョウといった広葉樹の樹液流の速度，S2 地点でトドマツの樹液流の速度をヒートパルス法で測定した。S2 地点のトドマ
ツ林は，上流域のトドマツ林からやや離れた個所に位置するが，樹齢はほぼ同じとみなせるので，ここに観測点を設けた。
トドマツ林の S2 地点に隣接する裸地には，風速計を地上 1.8 m に設置し 10 分間平均風速を得た。樹液流の速度の計測には，林電工業株式会社製 HP-2 型のヒートパルス式蒸散流速計を用いた。蒸散流速計のセンサーは，棒状ヒーター1 本と温度センサー2 本で 1 セ
ットになっており，本研究では 9 セットを用意して観測にあたった。ヒーターとセンサーを設置する際，専これらの仮定の下で，T についての解が次式のように得られる：
𝑇

𝑒𝑥𝑝

(3)

(2)

ここで，熱伝導係数 = kx/c = ky/c である。これから，上流側センサー(-0.5,0)での測定温度 Tup と下流側センサー(1.0,0)での測定温度 Tdown は，(2)式から容易に求まる。このとき，温度センサーTup, Tdown の時間変化から
Tup=Tdown となる時間を t0 とすると，
5

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

用ガイドを樹幹に固定しハンドドリルで樹木に垂直に穴を開けた。これによって，ヒーターと 2 本のセンサ
ー間の間隔が正確に 5 mm，10 mm で一直線に並ぶよ
うに挿入固定された(Fig. 2)。また，外部の気温変化の影響を防ぐため，樹木の外に出ているセンサーとヒー
ターの全体を断熱材の脱脂綿で包み，この上をビニー
ルシートで覆い，最後に遮光用のアルミホイルで包んだ。なお，ヒーターからの放熱時間は基本的に 2.0 秒間である。
本研究では，樹液流の速度の測定を樹木 1 本につきセンサー1～2 セットで行った。このための予備観測として，先ず代表的なトドマツとイタヤカエデを 1 本ずつ選定し，それぞれ同じ深さに 4 方位からセンサーを挿入して同時に樹液流の速度を測定した。
これにより，
方位によって測定流速に差がないことを確かめた。
Fig. 1 の T 地点では，気温・湿度・全天日射量・雨量を計測した。気温・湿度は地上 1.5 m 高で 30 分間隔
で，
全天日射量は地上 2 m 高で 15 分間隔で測定した。
なお，この気温・湿度データを用いて飽差(g m-3)を計算した。
Fig. 1 の黒プロット 8 か所と S1, S2 地点では，20 m
×20 m の範囲で毎木調査を行った。調査項目は，樹種と本数，胸高直径 DBH(cm)である。このとき，蒸散量を考える際，DBH<4 cm の樹木は群落への寄与が小さいと考え，カウントしなかった。なお，トドマツ林については，毎木調査は上流の 1 か所でのみ行った。これは，トドマツ林の大部分が 1957～1958 年の植林で，樹高と胸高直径の大きさから一斉林とみなせると考えたからである。

z に関する 2 次式で表したが，より正確な分布を得る
ため，次の 3 次式を用いた：
𝑉

𝑎′

𝑏′

′

𝑐′

′

′

(4)

ここで，a’, b’, c’：任意定数，sw’：代表木の辺材幅である。(4)式は，トドマツ，イタヤカエデ，オヒョウにそれぞれ異なる a’, b’, c’値を用いて適用し，シナノキについては，同じ散孔材のイタヤカエデの式を使用した。
通常観測では，樹液流の速度を 1～2 点の深さでのみ測った。ここでは，同じ樹種の深さ方向の流速分布
は，DBH の大きさ，つまり sw の大きさに関わらず等しいと仮定した。そして，代表木の流速分布式の(4)式をそれぞれ全ての樹木に適用した。例えば，代表木と異なる辺材幅 sw の樹木について，深さ z で樹液流の速度 V を測定した場合，
V と z の関係式 V = z3+z2+zの係数, , は，
𝛼

′

, 𝛽

′

, 𝛾

′

(5)

として(4)式の a’, b’, c’をそれぞれ変換して, , を求め
流速分布式とした。以上から，樹木 1 本の蒸散フラッ
クス F は，それぞれの流速分布式を辺材全体で積分し次式から求めた：𝐹

3.3. 樹木 1 本の蒸散量算定法
樹液流速度の測定は，樹木1本での1点の値なので，この樹木の蒸散量 Tr1 を求めるには樹液が流れる辺材部で積分する必要がある。これについては，樹木の深さ方向に流速分布を求め，それを辺材部全体で積分した。辺材幅 sw は，既出の成長錐でコアを採取して求めた。先ず，トドマツの場合は，含水率が心材では急激に小さくなることを利用した。広葉樹のオヒョウは辺材と心材で色が明らかに異なるので，肉眼で sw 値を求めた。シナノキとイタヤカエデの場合は，この識別ができないので，Wullschleger et al. [26]による sw～DBH間の経験式を利用した。

2𝜋 𝑟

𝑧 𝑉𝑑𝑧

(6)

ここで，r は樹木の半径(cm)である。ここでは，樹液流の速度を 30 分ごとに測定し，(6)式から 30 分ごとの蒸散フラックスを求めた。さらに，樹液流の速度が観測できた時間帯(7:00～19:00)について F を積分し，樹木
1 本の日蒸散量 Tr1 (L day-1)とした。
3.4. 蒸散量のスケールアップ
前節の日蒸散量 Tr1 (L day-1)を，流域全体の蒸散量 Tr(mm day-1)にスケールアップすることを考える。Hatton
et al. [27]は，観測した樹木の辺材面積，葉面積，胸高
直径，個々の樹木が占める面積(tree domain)をパラメーターとして，流域全体の蒸散量を求めた。その結果，辺材面積または胸高直径が最も適したパラメーターであることを見いだした。一般的に，辺材面積と胸高直径の間には強い相関がある。そのため，本研究では測定の容易さを考え，胸高直径 DBH をパラメーターとして選択した。
本研究では，ヒーター1 本，温度センサー2 本を 1 セ

流速分布については，代表的な樹木として，
DBH=42.5 cm のトドマツ，DBH=18.4 cm のイタヤカエ
デ，DBH=42.8 cm のオヒョウを選び，それぞれ異なる
深さで同時の樹液流の速度を測定した。約 1 週間の観測後に，それぞれの深さにおける樹液流の平均速度𝑉を計算した。Edward and Warwick [22]は，この𝑉を深さ6陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

ットとして，9 セットのみの観測で流域蒸散量を評価する。このため，Vertessy et al. [28]による roaming sensor
technique を利用して，9 セットのセンサーでより多くの樹木の樹液流の速度を測定する方法を考えた。
先ず，基準木を 1 本選び，これを含む 9 本の樹木にそれぞれセンサーを 1 セットずつセットした。これらの樹液流の速度を 1～2 週間観測し，その後，基準木以外のセンサーをそれぞれ別の樹木に移設して 1～2 週間観測した。測定を終了した樹木については，同時に得た基準
木のデータに回帰させた。測定期間外の値は，基準木の測定値を回帰式に代入して求めた。ここで Vertessy
et al. [28]と手法が異なる点は，基準木を 2～3 期間ごとに変えて行ったことである。この理由は，長期間センサーに設置するとドリルで穴を空けた部分にチロースが溜まり，樹液流の速度に影響するからである。
また，
基準木を変えることで，より広範囲の観測が可能になった。
基準木以外のセンサーを移設するまでの期間を 1 期とすると，トドマツ林で計 8 期，混交林で計 9 期の観測を行った。本数としては，トドマツ計 24 本，イタヤ
カエデ計 9 本，シナノキ計 8 本，オヒョウ計 6 本のほか，他の樹種の広葉樹 3 本の計 50 本を計測した。

/

= 0.826B3 1.234B2+1.135B+0.298 (0<B 1)

ここで，流域の日照時間 N は，流域の北東約 20 km にある札幌管区気象台での日照時間と全天日射量との関係を用い，T 地点の全天日射量から求めた。また，流域の可照時間 N0 は，T 地点における快晴時の全天日射量の時間変化から求めた。
(7)式中の空気力学抵抗 ra は，大気の安定度が中立であるとして，次式を使って求めた：

𝜎𝑇

1

1

↓

1

1

exp

𝐶

(14)

ここで，𝐶 𝐶∗ /2𝑘 ，無次元葉面積指数𝐶∗ 𝑐 ∙ 𝐿𝐴𝐼，
葉面抵抗係数 cd (～0.2)である。なお，LAI 値 6 は石井[30]を参考に決定した。(13)式中の zr と ur については，
S2 地点に隣接の裸地で計測した 1.8 m 高の風速 u に次式の速度対数則を適用した：

(8)

𝐴

(13)

ここで，k：von Karman 定数(=0.41)，zr：群落上の代表的な高度(m)，ur：zr における風速(m s-1)，d, zo:それぞれ
ゼロ面変位(m)と粗度長(m)である。ここで，ゼロ面変位 d の値ついては，近藤[29]による次式によって，対象流域の森林群落に対し葉面積指数 LAI=6 として
d/zr=0.73 を求め，北海道森林管理局・担当官（私信）による zr =16 m から d=11.6 m を得た：

𝑢 𝑧

ここで，：アルベド，：短波放射量(W m-2)，L ：下向き長波放射量，：ステファンボルツマン定数
=5.670×10-8 Wm，s：地表面温度(K)，である。
下向き長波放射量 Lは，近藤[29]に従い快晴時の下向
き長波放射量 Lf を求め，次式から算定した：
𝐿↓

/

𝑟

，蒸発散量kg m-2 sここで，：蒸発の潜熱（J kg-1）-3
1)，：大気密度kg m )，cp：定圧比熱(J kg-1 K-1)，：乾湿計定数hPa K-1，：飽和水蒸気圧曲線の傾きhPaK-1，Rn：正味放射量(W m-2)，D：飽差，ra：空気力学的抵抗(s m-1)，rc：群落抵抗(s m-1)，である。ここでは，蒸発散量 E を日単位(mm day-1)で求めるため，地中の貯熱量変化は無視した。(7)式の正味放射量 Rn は，次式で与えられる：
Rn = (1  L - Ts4

(11)
(12)

= 0.2235 (B=0)

(7)

/

(10)

↓ ⁄

ここで，T：日平均気温(K)，である。(10)式のは雲の効果を表す係数で，ここでは日照時間 N(h)と可照時間
N0(h)の比である日照率 B=N/N0 の関数として，次式によって求めた：3.5. ペンマン-モンティース法
本研究では，流域蒸発散量を評価する方法として，
次のペンマン-モンティース式を用いた：

𝜆𝐸

↓⁄

𝐴

∗

ln

(15)

ここで，𝑧 𝑑 𝑧 ，u*：摩擦速度(m s-1)，である。
なお，このとき，地表面の粗度長 z0s として草地の代表
値 0.1 m を用い，これから次式を用いて群落の粗度長
z0 = 1.2 m を求めた[29]：(9)7陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

ln
.

1

exp

𝐶

ln

.

∙ exp

2𝐶

(16)

葉の気孔に関するパラメーターである群落抵抗 rc については，
次式で示すStewart [31]のモデルを適用した：
𝑓 𝐾↓ 𝑓 𝐷 𝑓 𝛹

𝑟

(17)

図 3. 毎木調査に基づく(a)トドマツ林と(b)混交林の胸高
直径 DBH の頻度分布。
Fig. 3. Frequency distributions of DBH (cm) in (a) Sakhalin fir
forest and (b) mixed forest by tree surveys.

ここで，rmin：最小群落抵抗(s m-1)，𝛹 ：土壌水分ポテ
ンシャル(MPa)，である。この式は，短波放射量 K↓，
飽差 D，土壌水分ポテンシャルΨ のそれぞれの関数 fi,
f2, f3 が，全て独立に群落抵抗 rc に影響することを表す。
観測期間では総雨量 477 mm，平均日雨量 6.2 mm で，雨量と降雨イベント数が多かったため，水分ストレスによる蒸散量の減少は全く見られなかった。このことから，Ψ の関数 f3 は 1 とした。なお，(17)式の関数 f1,
f2 は，それぞれ以下の式で表せる：

𝑓 𝐾↓

↓
↓

小さい値に偏った分布を示した。
4.2. 樹木 1 本の蒸散量
Fig. 4 は，トドマツとオヒョウに対して測定された
辺材幅 sw と胸高直径 DBH との関係を示す。sw は，どちらも DBH 増加と共にその増加率は減少したが，全体としてはオヒョウの方が値が大きかった。
一般的に，sw は(20)式のように DBH のべき関数として表せられる[26]：(18)

sw = a(DBH)b

𝑓 𝐷

(20)

(19)
ここで，a, b：樹種によって変化するパラメーター，である。図から，トドマツでは a=0.24, b=0.69，オヒョウ
では a=0.39, b=0.69 と求められた。今回，測定できなかったシオノキとイタヤカエデについては，
Wullschleger et al. [26] が 求 め た そ れ ぞ れ a=2.052,
b=1.654 と a=0.90, b=1.84 を用いて sw を見積もった。

ここで，k1, k2, k3：任意パラメーターである。(18)式と(19)式は，それぞれ Stewart [31]と Zhang et al. [32]によって用いられた関数である。ここでは，Zhang et al.[32]に従い，rmin，k1, k2, k3 の 4 つのパラメーターの値について，実測の樹液流速から求めた蒸散量とフィッティ
ングさせる方法で決定した。
4. 結果
4.1. 毎木調査による胸高直径の分布
Table 1 に，毎木調査によって得た植生の特徴を示す。
なお，立木密度は DBH≥4.0 cm の樹木を対象とした値である。この表より，混交林の DBH はトドマツ林に比べバラつきが大きいこと，混交林でもトドマツが全体の約 23 %を占めることがわかる。また，混交林での平均胸高直径は，トドマツ 20.2 cm，シナノキ 13.9 cm，イタヤカエデ 11.2 cm，オヒョウ 12.8 cm であった。
Fig. 3 は，(a)トドマツ林と(b)混交林における胸高直径 DHB の頻度分布を示す。トドマツ林は左右対称に近い分布を示すが，混交林では DBH がより広範囲で
8

図 4. 毎木調査に基づく(a)トドマツと(b)オヒョウの胸高
直径 DBH と辺材幅 sw との関係。
Fig. 4. Relations between DBH and sw for (a) Sakhalin fir and(b) Manchurian elm by tree surveys.

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)
表 1. 毎木調査で得られた針葉樹林と混交林の分布と特徴．樹木の和名は本文を参照．
Table 1. Distributions and features of needle-leaved and mixed forests obtained by tree surveys. The label “Elm” shows“Manchurian elm” (see text).
Area Standing density DBH Mean
DBH
Sakhalin fir Japanese
2
-2
(km )
(No. km )(cm)
Std. dev. (cm)(%)
linden (%)

Paintedmaple (%)

Elm Others(%)(%)

Needleleaved forest

0.63

1.120×105

20.6

7.6

100

−

−

−

−

Mixed forest

0.48

1.448×105

15.5

10.3

22.7

20.8

17.3

9.2

30

Fig. 5 に，代表木のトドマツとイタヤカエデの無次
元深さ z/sw についての流速分布とその近似曲線を示す。なお，図の白丸にあるように，V=0 (z/sw =1)と仮定している。ここで z は樹皮の下の形成層からの深さを指す。また，図から，流速の極大はトドマツで z/sw～
0.4，イタヤカエデで z/sw ～0.3 の位置にあり，辺材の中心 z/sw= 0.5 より外側にあることがわかる。なお，辺
材幅はトドマツの場合は 3.0 cm，イタヤカエデの場合
は 7.2 cm であった。しかし，センサー長は 5.0 cm であり，測定可能な最大深さは 3.6 cm のためイタヤカエデの場合は深さ 3.6～7.2 cm は測定できなかった。このため，イタヤカエデの場合は，心材では樹液流の速度はゼロとし，z/sw> 0.5 の領域は直線で外挿した（Fig. 5bの点線）。

準となる樹木を選び，ここで固体された 1 セットは動かさずに，他の 8 セットを様々な樹木に固定して同期間の樹液流速を測定し，基準木と他の樹木との間で相関をとってローミングすることを考えた。ここでは，
これを roaming sensor technique と呼ぶこととする。
Fig. 6 は，得られた樹木 1 本の日蒸散量 Tr1(L day-1)
について，基準木のトドマツ E と他のトドマツ F, G, Hとの関係を示す。なお，この観測期間は 7 月 31 日～8
月 10 日であるが，他の期間でもそれぞれの回帰式は全て y=a1x（a1：任意定数）で表すことができた。その決定係数の平均は 0.88 で，最低値は 0.61 であった。一方，広葉樹では決定係数の平均は 0.66 で，トドマツより小さな値をとった｡
ここで，相関の悪かった期間の例として，9 月 23 日
～10 月 9 日の基準木のトドマツ V と他のトドマツ Y,
Z との関係を Fig. 7 に示す。この期間におけるトドマ
ツ V とトドマツ Y, Z との相関は，R2 でそれぞれ 0.89,
0.70 であった。しかし，図にあるように，9 月 23 日～
10 月 2 日（図の〇と□）と 10 月 3 日～9 日（図の●と■）に分けると後者の期間で相関が悪くなり，データに異常性が認められた。なお，図中の回帰式と決定係

図 5. (a)トドマツと(b)イタヤカエデにおける樹液流の速度 V の無次元深さ(z/sw)についての分布．
Fig. 5. Distributions of sap flow velocity, V, for (a) Sakhalin firand (b) painted maple. z: depth from the cambium.

Fig. 5 の近似曲線を用い，(6)式から樹木 1 本の蒸散量を計算した。
なお，
通常は辺材の外側に樹皮があり，この部分は樹液は流れない。DBH と樹皮の厚さの関係を調べると，樹皮の厚さは DBH の約 5%であった。このため，(7)式での樹木の半径 r = 0.95(DBH)/2 としている。
図 6. 日蒸散量に関する基準木のトドマツE とトドマツF,G, H との関係．
Fig. 6. Relations between a reference tree, Sakhalin fir “E”, andSakhalin firs “F”, “G”, and “H” for daily transpiration.

4.3. roaming sensor technique による補間
今回の観測では，１本のヒーターと 2 本の温度センサーを 1 セットとして，全 9 セットでできるだけ多くの樹木の樹液流速を測ることを考えた。このため，基9

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

月 4 日～ 5 日，10 月 7 日～9 日の 10 日間で，この 3
期間前はそれぞれ 9 mm, 17 mm，1 mm の雨量があり，無降雨日の平均気温は 9.2 ℃，10.7 ℃，8.5 ℃であった。つまり，10 月 2 日を境に気象条件が大きく変わる
ことはなかった。また，飽差について，10 月 2 日を境に 10 月 3～5 日でやや減少しているが，10 月 6～9 日の増加時でもトドマツ V の流速は減少傾向にあった。
流れの阻害がある場合，樹液流の速度を回帰式 y=ax
で求めることは困難になる｡このため，10 月 3 日以降のトドマツ V の樹液流の速度について補正を行った。
先ず，Fig. 9 のように，トドマツ Y とトドマツ V の 30
分ごとの樹液流の速度を比較した｡ここでは，
この図における 9 月 23 日～10 月 2 日について，図の回帰式（黒い太線が回帰曲線）を用いて 10 月 3 日以降のトドマ
ツ V の樹液流の速度を計算した。
今回のトドマツ V のように流速が減少する現象は，トドマツで 4 本，広葉
樹で 4 本で認められた。
一方，広葉樹については，y= a1x による近似の決定
係数 R2 が小さかった例を Fig. 10 に示す。サンプル数
は少ないが，決定係数は 0.94, 0.95 で相関がみられた。
このとき，オヒョウ N4 については，y 切片はゼロではなく，回帰式 y= a1x +b1 の適用が妥当であることがわかった。この結果は，基準木の中でオヒョウ 3 本にのみ観測されたが，その理由については不明である。
以上から，広葉樹の蒸散量データの補間には，y=a1xと y=a1x+b1 の 2 通りの回帰式を使用することとした。

図 7. 日蒸散量に関する基準木のトドマツ V とトドマツ
Y, Z との関係．〇，□：観測期間 9 月 23 日～10 月 2 日●，■：観測期間 10 月 3 日～9 日．
Fig. 7. Relations between a reference tree, Sakhalin fir “V”, andSakhalin firs “Y” and “Z” for daily transpiration. 〇, □:observed for 23 Sept – 2 Oct. ●, ■: observed for 3 – 9 Oct.

数は 9 月 23 日～10 月 2 日の結果に対してのものである。この異常性を確かめるため，トドマツ V, Y, Z の 30
分ごとの樹液流の速度(cm h-1)と飽差(g m-3) の変動をFig. 8 に示す。なお，樹液流の測定した深さは，それぞれ 1.8 cm，1.6 cm，1.4 cm で，無次元深さ z/sw = 0.66,
0.53, 0.52 の位置である。ここで，トドマツ V に着目す
ると，10 月 2 日を境に明らかに他の 2 本に比べ樹液流速速が小さくなっている。この原因としては，黒田[19]が指摘する“外傷によって仮道管内で気泡が発生し流路が狭められた”ことに因って樹液流が阻害されたと考えた。なお，9 月 23 日～10 月 9 日の 17 日間について，無降雨日（観測日）は 9 月 27 日～10 月 1 日，10

4.4. 蒸散量と胸高直径との関係
それぞれの樹木の胸高直径 DBH と日蒸散量 Tr1 との関係を調べた。なお，データは無降雨日(49 日間)のみのデータを使用した。Fig. 11 に，トドマツ A～Z（計測できなかった C, N を除く 24 本）の無降雨日での平均

図 8. 2001 年 9 月 23 日～10 月 9 日におけるトドマツ Y, V, Z の樹液流速度と飽差の時間変化（30 分間隔）．
Fig. 8. Semi-hourly series of sap flow velocity for Sakhalin firs “Y”, “V” and “Z” for September 23 – October 9, 2001.

10

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

図 9. 2 つの観測期間の樹液流の速度に対するトドマツ Yとトドマツ V との関係．〇 9 月 23 日～10 月 2 日 ●10
月 3 日～9 日
Fig. 9. Relations between Sakhalin firs “Y” and “V” for sap flow
velocity (cm h-1) of September 23 – October 2(white circles) andOctober 3–9 (black circles).

図 11. 7 月 9 日～9 月 20 日間の無降雨日におけるトドマ
ツ A～Z の平均蒸散量 Trm と DBH との関係．
Fig. 11. Relation between DBH and mean transpiration Trm for
Sakhalin firs A – Z observed on non-rainfall days of July 9 –September 20.

図 12. トドマツ A～Z に対する関係式 Trm=c(DBH)d での係数 c の時間変化．
Fig. 12. Temporal variation of coefficient c in the equation
Trm=c(DBH)d for Sakhalin firs A - Z.

図 10. 蒸散量に対する基準木のウダイカンバ N1 とシナ
ノキ N5，オヒョウ N4 との関係．観測期間 9 月 3 日～8日．
Fig. 10. Relations between a reference tree, monarch birch N1,
and Japanese linden N5 or Manchurian elm N4 for transpiration(L day-1) of September 3 – 8.

日蒸散量 Trm と DBH との関係を示す。両者の関係は，次のべき関数で表せられる：
Trm=c(DBH)d(21)

ここで，c, d：任意定数である。Fig. 11 では，R2=0.93 の高い相関が見られる。前節で示したように，トドマツ
ではそれぞれの樹木間の日蒸散量比が一定，つまり
y=a1x で近似できると考えている。このため，各日の日
11

蒸散量 Tr1 と DBH の回帰曲線は，樹木間で全て等しい
ことになる。つまり，7 月 9 日～9 月 20 日の無降雨日において，平均蒸散量 Trm～胸高直径 DBH 間の関係式の係数 d は，常に一定で 2.864 になる。また，各日のDBH と日蒸散量 Tr1 の回帰式は(ΣDBHi/n, ΣTri/n)を通る。ここで，n は，サンプル樹木の数である。この関係を利用し，係数 d = 2.864 から各日での回帰曲線の係数 c を求めた(Fig. 12）。
次に，広葉樹について説明する。広葉樹は樹種ごとのサンプル数が少ないため，
全樹種をまとめて扱った。
Fig. 13 は，広葉樹に対する 7 月 9 日～9 月 20 日間の無降雨日の平均日蒸散量 Trm と DBH との関係を示す。ト陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

ドマツに比べると回帰曲線の傾き d は小さいが，係数c は大きくなった。また，トドマツに比べバラツキは
大きく，R2=0.71 となった。さらに，Fig. 10 のオヒョウのように，広葉樹ではそれぞれの樹木間の蒸散量比が必ずしも一定ではなかった。このため，トドマツと異なり，係数 c ばかりでなく係数 d も日によって変化する。従って，ここでは Trm～DBH 間の関係式(21)式の係数 c，d を一日ごとに求めた。係数 d については，Fig.
14 のようになった。係数 d は，8 月下旬までほぼ一定であったが，それ以降は大きなバラツキと共に増加した。観測期間中，係数 d は平均 2.305, 標準偏差 0.0181であった。
各日の係数 c の時間変化を Fig. 15 に示す。
7 月 23 日～31 日にかけて最大となり，この時期に樹液流が活発になることを示唆する。その後，9 月にはトドマツの場合と同様，減少傾向が見られた。また，その値はトドマツに比べほぼ 1 オーダー大きな値を示した。

図 13. 7 月 9 日～9 月 20 日間の無降雨日における広葉樹の Trm と DBH との関係．
Fig. 13. Relation between DBH and mean transpiration Trm forbroad-leaved trees observed on non-rainfall days of July 9 –September 20.

図 14. 広葉樹における関係式 Trm=c(DBH)d での係数 d の
時間変化．
Fig. 14. Temporal variations of the coefficient d for broad-leaved
trees.

図 15. 広葉樹に対する関係式 Trm=c(DBH)d での係数 c の
時間変化．
Fig. 15. Temporal variation of coefficient c in the equation
Trm=c(DBH)d for broad-leaved trees.

図 16. (a) トドマツと(b) 広葉樹の流域蒸散量に占める DBH 別蒸散量の割合．
Fig. 16. Ratio of transpiration at each DBH to catchment transpiration for (a) Sakhalin fir and (b) broad-leaved trees.

12陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

図 17. 流域蒸散量(mm day-1)の季節変動．トドマツは混交林に分布するものを加えた値．広葉樹はトドマツ以外の値．
Fig. 17. Seasonal variations of catchment transpiration (mm day-1) for Sakhalin fir and broad-leaved trees. Sakhalin firs in themixed forest are also included.

Tr/Epm 値は常に Tr/Epm<1 で，Fig. 18 では 8 月下旬まで漸減傾向をとりながらも 0.28～0.45 の範囲で比較的一定の値をとった。しかし，9 月に入るとバラつきが大きくなり，9 月 9 日で最大値 0.63 を示した。また，9 月に入ると気温と Tr/Epm 比の変動は類似し，気温上昇の時は Tr/Epm 比も上昇する傾向が見られた。Epm 値は，林床にあるクマイザサの蒸発散量とその下の土壌水分の蒸発量を合計した蒸発散量 Ef を含むため，Epm>Tr となるのは妥当と考える。しかし，Ef の占める割合は一般
的に全蒸発散量の 5～20 %であり[33, 34]，Tr 値は Ef を含まないとしても非常に小さいと言える。
なお，Fig. 18 では，降雨前後に Tr/Epm 比が大きく変化することはなかった。これは，降雨後の遮断蒸発量
による蒸発散量の増加と降雨後の土壌水分増加に伴う

4.5. 流域蒸散量
前節でトドマツと広葉樹に対し，それぞれの Trm～
DBH 関係式として，(21)式中の係数 c, d を求めた。こ
こでは，これらの関係式と Fig. 3 の DBH 分布を用いて流域全体の蒸散量 Tr(mm day-1)を計算した。
Fig. 16 は，(a)トドマツと(b)広葉樹に対する流域蒸散量に占める DBH 別の蒸散量比を示す。
トドマツでは，
DBH=24～28cm の樹木が流域蒸散量に大きく寄与しており，全体として尖度の高い分布となった。一方，
広葉樹では各 DBH でなだらかな分布を示している。
Fig. 3 との比較から，トドマツ，広葉樹ともに DBH の大きな樹木の蒸散量が相対的に大きな割合を占めていることがわかる。
Fig. 17 は，トドマツと広葉樹の流域蒸散量の季節変動を示す。トドマツと広葉樹は，共に夏から秋にかけて蒸散量が減少している。観測期間 7 月 9 日～9 月 20
日の無降雨日における平均蒸散量はトドマツで 0.53mm day-1，広葉樹で 0.27 mm day-1 であった。
4.6. ペンマン-モンティース法との比較
Fig. 18 に，ペンマン-モンティース法で計算した流域
蒸発散量 Epm (mm day-1)とヒートパルス法による流域
蒸散量 Tr (mm day-1)との比，および平均気温と雨量の経日変化を示す。なお，図中のプロットは，全て無降
雨日の値を表す。また，この時に対応する Tr～Epm 間の関係を Fig. 19 に示す。なお，この時，ペンマン-モ
ンティース法で用いた(15)～(17)式での 4 つのパラメ
ーターの最適値は，rmin=120, k1=5.0, k2=0.01, k3=0.11(R2=0.83)と求まった。また，Fig. 19 では，Tr, Epm は共に原点を通ると考え，回帰式はべき関数 Tr=0.44Epm0.66(R2=0.87, p<0.01) で表している。これらの図から，

図 18蒸発散量 Epm と日蒸散量 Tr の比と日平均気温・日
雨量の経日変化．Tr 値は全て無降雨日のもの．
Fig. 18. Temporal variations of Tr/Epm ratio, daily mean air
temperature and daily rainfall. The Tr values were obtained on
non-rainfall days.

13

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

and Whitefield [18]は，xy 平面の 2 次元数値計算によって，従来のヒートパルス法では樹液流の速度を 50～
75 %過小評価することを報告している。本研究では，
深さ方向の流速分布(4)式を利用し，陽的差分法による3 次元数値計算を行った。
5.2. 数値計算での物性に関する条件
数値計算における格子間隔は x, y 軸方向 0.04 cm，z
軸方向 0.08 cm とした。
センサーは全て直径 0.2 cm で，その熱特性については Fig. 20 のような配置で考えた。
ヒーターから放出される熱は，その中心にある格子点（■点）にのみ与えた。温度センサーの測定点は，中心の格子点（●点）の温度を使用した。注意点として，観測では下流側と上流側の温度センサーとヒーターの位置は，それぞれ x= - 0.5 cm, 0 cm, 1.0 cm にあった(Fig.
2)。しかし，計算では格子間隔の都合上 x = - 0.48 cm, 0
cm, 0.96 cm の位置で計算した。樹木の生材の密度g(kgm-3)と比熱 cg(J kg-1 K-1)は，それぞれ次の(22)，(23)式を使って求めた：
図 19.蒸発散量 Epm と日蒸散量 Tr との関係．
Fig. 19. Relation between evapotranspiration Epm andtranspiration Tr.

蒸散量増加とが比較的同期していることを示唆する。
5.ヒートパルス法に対する補正  gmg/Vg
5.1. ヒートパルス法の問題点
これまで，ヒートパルス法では樹液流の速度が過小評価されることが多数報告されている[18, 21-24]。3.2.
節で述べたように，ヒートパルス法によって樹液流の速度を求める場合，2 次元の熱移流拡散方程式の解析
解(2)式が基本となっており，その解析解を得るには，
初期条件以外に次の 3 つの仮定が前提にあった：(22)

cg = 4.1868103(  

ここで，mg, Vg：生材の質量(kg)と体積(m3)，：含水率(i) ヒーターとセンサーの体積は無限小である(ii) 樹液流の速度 V は xy 平面上で一様に分布する(iii) 樹木は熱的に等方均質である(kx = ky)
Swanson and Whitefield [18]は，この過小評価の原因を調べ，実際の樹液流では上の 3 条件が成り立たないことを突き止め，その補正方法を検討した。(i)については，センサーの直径 0.2 cm に対しセンサーの間
隔 0.5 cm，1.0 cm であり，センサーの大きさを無視することはできない。このことから，条件(iii)に関連し，センサーと樹木の熱特性の違いを無視できない。
また，
条件(ii)に関し，センサーが樹液流の流れを阻害することが起こりうる。
こうして，2 次元の熱移流拡散方程式での熱容量や熱伝導率は樹木の値を用いるのではなく，センサーの箇所ではセンサー物質の値を用いるべきである。
また，樹木へのセンサー挿入によるキズが樹液流を阻害するので，これに関する数値計算が必要である。Swanson
14

図センサーの熱特性に関する条件．
a樹幹と計算領域．
中央の矢がヒーター，他の矢は温度センサー．点線で囲んだ部分が計算領域．b数値計算で使用したヒーターとセンサーのグリッド配列．■ヒーターと●温度センサーの中心位置．円筒状のセンサーを計算では格子点の集合としている．
Fig (a) Calculation domain in a trunk. The central arrow andthe others show a heater and sensors, respectively. (b) Central
location of a heater (black square) and sensors (black circles), andtheir arrays on grids in the simulation.

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)(%)である。(23)式は伏谷ほか[35]による経験式で，比熱cg は樹種と密度によらず含水率の関数と考えられた。
熱伝導率 kx, ky, kz について，接線方向の ky と放射方向の kz は，一般的にほぼ等しいとされ，x 軸方向の kx
は，それらの 2 倍程度と考えられている[35,36]。木材は多孔質で熱の絶縁性が高い。その熱伝導率は木材内の空隙に関係するため，密度の関数と考えられる。気乾状態（含水率 12 %）での kx について，次式のような実験式が導かれている[35]：kx(12)= 4.1868102(0.022+1.6810-4



ここで，kx(12)：含水率 12 %の時の x 軸方向熱伝導率，
：含水率のときの密度，である。さらに，木材
よりも水の熱伝導率は大きいため，生材の熱伝導率 kxは含水率の関数と考えられる。なお，ここでの含水率は次式で定義される重量比率である：
・mw/m0)(25)
ここで，mw：木材に含まれる水の重量，m0：木材の乾燥重量である。また，木材で用いられる容積密度(bulkdensity)R(kg m-3)は，次式で定義される：

R=m0/Vg

いた。また，容積密度と含水率については，林業試験場の木材部・木材利用部[37]のデータを使用した。また，センサーの容積比熱（volumetric specific heat）と熱伝導率は，Table 2 に示す SUS304 ステンレスの値を用いた。
5.3. 樹液流の阻害に関する条件
樹液流の阻害に関し，今回の数値計算で用いた 2 つの条件を Fig. 21 に示す。A タイプは Barett et al. [23]を参考にして設定した。彼らはドリルによるキズ部分を薄切りにし，光学顕微鏡と走査型電子顕微鏡でキズの形状を調べた。その結果，x 軸方向に約 0.4 cm，y 軸方
向に約 0.3 cm にわたって導管の損傷があることを見いだした。A タイプの流速ゼロ領域はこのキズの範囲を考慮したものである。B タイプは，Swanson andWhitefield [18]が 2 次元数値実験で用いたもので，x 軸方向に比べ y, z 軸方向の流れが小さいため，センサー
間では樹液流はないとした場合である。ここでは，キズの幅を ww として，センサーの直径と同じ 0.20 cm，Barrett et al. [23]を参考に 0.28 cm，さらに大きなキズとして 0.36 cm といった 3 通りを考えた。z 軸方向の流速分布については，実測に基づく(5)式を使用した。

(26)

ここでは，次の実験式から kx を求め，ky, kz は kx の 1/2の値とした：
kx=kx(12){1 0.011(12   
与えられたトドマツとセンサーの物性値を Table 2に示す。なお，今回の計算では全てトドマツの値を用
表 2. 数値計算で使用されたトドマツとセンサーの物性値．
Table 2. Thermal properties of Sakhalin fir and sensors appliedin the simulation.
Heater or

Sapwood

Heartwood

Bulk density (kg m-3)

362

338

-

Water content (%)

212

76

-

3.69

1.53

4.03

0.716
0.358

0.380
0.190

13.0
13.0

Volumetric specific heat(kJ m-3 K-1)
Thermal conductivity(W m-1 K-1)
kx
ky, kz

図 21. 樹液流の阻害に関するモデル条件．灰色部分は流速ゼロの領域．A タイプ：Barrett et al. [22]を参考にしたゼロ領域，B タイプ：Swanson and Whitefield [18]が与えたゼ
ロ領域．
Fig. 21. Two models related to the obstruction of sap flow. Gray
area shows zero velocity zone. Type A and type B were given byreferring to Barett et al. [23] and Swanson and Whitefield [18],respectively.

sensors15

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

5.4. 計算方法
数値計算では，次の 3 次元移流拡散方程式を陽的差分法で解いた：

𝑐

𝑘

𝑘

𝑘

𝑉𝜌 𝑐

𝑄 (28)

初期条件：T(x,y,z,t)=0 (t=0)
境界条件：
, ,

,

ℎ

𝑇 𝑥, 𝑦, 𝑧 , 𝑡

𝑇 𝑡(z0 平面センサー部分), ,,T(x,y,z,t)=00(z0 平面センサー部分以外)(その他)

計算領域は，Fig. 22 に示すように x,y,z 方向にそれぞれ 4.8 cm2.4 cmcm とした。zo 平面が樹幹と外気の境界面であるとし，z 方向に樹皮 0.8 cm 厚，辺材 1.6
cm 厚，心材 1.6 cm 厚とした。樹液流は辺材領域のみに与え，樹皮と心材の領域は流速ゼロとした。なお，
容積比熱と熱伝導率について，樹皮領域では心材領域と等しい値を与えた。
境界条件は 3 つの部分に分け，z0 平面のセンサー物
質が存在する領域(Fig. 22 の黒い部分)は外気の温度 Ta
によって冷却される条件を与えた。外気温度は，基準温度として Ta(t)=0 で一定とした。熱伝達率 h は，放熱物質の形状や周囲の気体の流れ方によって変化する。
一般的には，空気の自然対流では h=4.2 W m-2 K-1 程度と考えられるが，ここでは，センサーの一部が z0 平面
より外側に存在するため h= 4.2W m-2 K-1 とした。
なお，h = 4.2～4.2 W m-2 K-1 の範囲で h を変えても，計算結果に大きな差はなかった。これは，z0 平面での温度は外気温度 Ta とほぼ等しいことを意味する。
次に，z0 平面のセンサー物質が存在しない領域（Fig. 22の灰色部分）では，断熱条件を与えた。これは，樹皮の最も外側では含水率や熱伝導率が非常に小さいと考えたからである。
数値計算におけるセンサー挿入深さは，ヒーターの熱源が辺材領域の中央(z=1.6 cm)となるように設定した。この z=1.6 cm に適当な樹液流の速度 Vimp を与え，(4)式の流速分布から辺材領域におけるそれぞれの深さでの樹液流の速度を計算して与えた。ここで，与えた Vimp は 5～50 cm h-1 の範囲で 5 cm h-1 刻みとした。また，熱量 Q については，一部の計算を除き，実際に利用した観測機器と同様の 2.0 秒間の矩形波の形で与えた。ヒーターの熱源格子点に熱量を与え，上下の温度センサーの温度が等しくなるまでの時間 t0 を計算し，(3)式から樹液流の速度 Vcal を見積もった。そして Vimpと Vcal を比較し，樹液流の速度がどの程度影響されるかを調べた。
数値計算は Swanson and Whitefield[18]と同様，次の 4
つのケースに分けて行った：
Case 1: センサーの熱特性と樹液流の阻害を共に考慮
しない場合
Case 2: センサーの熱特性のみ考慮する場合
Case 3: 樹液流の阻害のみ考慮する場合
Case 4: センサーの熱特性と樹液流の阻害の両方を考慮する場合
Case 1 では，Vimp と Vcal との間にあまり差はないと推測されるため，(3)式を導いたときの仮定の一つである“瞬時の放熱”について，放熱時間を 0.1～2.0 秒の範囲で変えて，その影響を調べた。
5.5. 計算結果
Fig. 23 に，(a) A タイプ対する Case 1～3 と(b)A,と Bの両タイプに対する Case 4 についての計算結果を示す。なお，放熱時間は全て 2.0 秒である。Case 1 では，計算された樹液流の速度 Vcal は Vimp と 1:1 に対応する
かそれに近かった。しかし，Vimp が大きいとき Vcal はやや過大評価する結果を与えた。これは，樹液流の速度が大きいと上下温度センサーが等温になる時間 t0 は小さくなり，熱の矩形波の影響が出てくることが考えられる。実際，放熱時間を 1.0 秒，0.5 秒，0.1 秒と短くするとこのズレも小さくなった。Case 1 に対し，センサーの熱特性を考慮した Case 2，および樹液流の阻害のみを考慮した Case 3 では，逆に過小評価する結果となった (Fig. 23a)。特に Case 3 は大きく過小評価しており，相対的に阻害の影響が大きいことがわかる。

図 22. 計算領域．z0 平面（灰色の面）は樹幹と外気の境界面．
Fig. 22. Calculation domain. The z0 plane shows an interface
between a trunk and the outer air.

16

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

図 23. 与えた樹液流の速度 Vimp と計算された樹液流の速度 Vcal との関係．(a) A タイプの Case 1～3 の場合，(b) Case
4 の場合（Fig. 21 の阻害条件 A タイプ, B タイプ別の結果）
Fig. 23. Relations between a given flow velocity Vimp and calculated velocity Vcal. (a) Case 1-3 for Type A and (b) Case 4, where
flow velocity was calculated for each obstruction of Type A and Type B in Fig. 21.

Case 4 では，阻害パターン Type A と Type B で影響が大きく異なることがわかる(Fig. 23b)。例えば，樹液流の速度 Vimp=50.0 cm h-1 のとき，Type A では，34.0 cm
h-1，Type B では 22.7 cm h-1 (ww=0.20 cm)，18.5 cm h-1(ww=0.28 cm)，15.8 cm h-1 (ww=0.36 cm)となった。TypeB では，キズの幅 ww が大きいほど 1:1 直線からのズレが大きくなり，過小評価される割合が大きくなった。
なお，樹液流を阻害する要因として，センサーを挿入する深さも関係すると考え，Type B の ww=0.20 cm について挿入深さ 1.2 cm～2.16 cm と変化させて計算した。その結果，Vcal 値はほとんど変わらず，挿入深さは
過小評価の程度に影響しないことがわかった。
5.6. 樹液流の速度に対する補正
前節で得た結果から，樹液流の速度に関する補正式を求めた。補正後の樹液流の速度 Vcor，観測値を Vobs として，補正式は Vcor=Vobs)2+Vobsの二次式で表すこととした。
樹液流の阻害パターンによる係数，を Table 3 に示す。Table 3 の係数を用いた 4 つの補正式により観測
表 3. 樹液流阻害パターン，A タイプと B タイプに対して得られた補正式の係数(Fig. 21 参照)．括弧内は ww 値．
Table 3. Coefficients and in the correction formulaVcor=Vobs)2+VobsValues in the parentheses show
ww values in Fig. 21.
Type AType B(0.20 cm)

Type B(0.28 cm)

Type B(0.36 cm)0.00110.0296

0.06260.1044

1.43731.5341

1.54171.517817

した樹液流の速度を補正し，これによって流域蒸散量
Tr を見積もった。また，この流域蒸散量をペンマン-モ
ンティース法による流域蒸発散量 Epm と比較した。
Fig. 24 は，4 つの補正式による蒸散量 Tr と Epm との関係を示す。(a)の A タイプでは，Tr は補正後にやや大きくなったが，明らかに Epm より小さい。この時，観測期間中の平均蒸散量は補正前で 0.72 mm day-1，補正後で 0.95mm day-1，Epm は 2.03 mm day-1 であった。つまり，補正前で Tr/Epm=0.32，補正後でも Tr/Epm=0.47 で流域蒸散量はかなり小さく見積もられた。B タイプの
ww=0.20 cm，0.28 cm では，補正後の流域蒸散量はそれ
ぞれ 1.46 mm day-1，1.76 mm day-1 であり，それぞれTr/Epm=0.72，0.87 であった。また，ww=0.36 mm では補正後 2.13 mm day-1 となり，Epm を上回った。以上から，林床からの蒸発散量が流域蒸発散量の約 5～20 %を占める[33,34]ことを考えると，ww=0.28 mm の条件が最も妥当と考える。
本研究では，クマイザサを含む林床からの蒸発散量
Ef を見積もっていないため，樹液流の速度の実測に基づく流域蒸散量 Tr の評価には，Epm 評価に含まれる問
題点（特に，
空気力学抵抗 ra と群落抵抗 rc の評価方法）と共に Ef 評価による検証が必要である。
樹液流の測定法として，他にグラニエ法や幹熱収支法がある。しかし，前者では Granier [38]が提案した経験式中の 2 つの係数が樹種によって大きく変化すること，後者では与えた熱量に対する応答遅延による樹液
流速の過小評価など，
これらの手法にも問題点がある。
このため，観測された樹液流速に対し単木でのポトメーター法[13]やライシメーター法[39]などによる検定をしつつ，それぞれの手法の特徴を生かした流域蒸散陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)

図 24. ヒートパルス法による流域蒸散量 Tr とペンマン-モンティース法による流域蒸発散量 Epm との関係．
(a) A タイプ, (b) B タイプ(ww=0.20 cm)，(c) B タイプ(ww=0.28 cm)，(d) B タイプ(ww=0.36 cm)での補正式を使用した場合（実線はその回帰曲線）．〇補正前 ●補正後 点線は 1 対 1 対応直線．
Fig. 24. Relations between transpiration Tr and evapotranspiration Epm from the Penman-Monteith method. Correction
formulae based on (a) Type A, (b) Type B (ww=0.20 cm), (c) Type B (ww=0.28 cm), and (d) Type B (ww=0.36cm) arealso shown. Solid lines and dotted lines indicate regression curves and 1:1 lines, respectively. 〇before correction, ●after
correction

マン-モンティース法による日蒸発散量の28～45 %で，林床からの蒸発散量を考慮してもかなり過小評価して
6. まとめいると判断された。つまり，解析解による樹液流の速度は実際の流速をかなり小さく見積もっていることが樹液流測定法の一つであるヒートパルス法によって， 推察された。そこで，本研究では，実測された樹木の流域蒸散量を見積もることを試みた。従来のヒートパ
深さ方向の流速分布を利用し，有限大のセンサーを考
ルス法では，次の 4 つの仮定のもとで熱の 2 次元移流
慮した 3 次元移流拡散方程式を陽的差分法により数値
拡散方程式を解析的に解くことで樹液流の速度を求め的に解くことを考えた。このとき，センサーと樹木のていた：物性については実測値を与え，樹液流を阻害するセン・熱量 Q は x = y = 0, t = 0 で瞬間的に与えられる
サー挿入時のキズの大きさを変えて，樹液流に与える・ヒーターとセンサーの体積は無限小である阻害効果を調べた。その結果，キズの幅が 0.28 cm の・樹液流の速度 V は xy 平面上で一様に分布する阻害効果を与えたとき，補正された樹液流の速度によ
・樹木は熱的に等方均質である(kx = ky)
る流域蒸散量 1.76 mm day-1 と求まり，この値が最も妥当であると判断された。このとき，ペンマン-モンティ
実際には，熱量 Q は 2 秒間与えられ，ヒーターとセース法による流域蒸発散量は 2.03 mm day-1 と評価された。他方，対象の森林帯の林床にはクマイザサが繁茂
ンサーの直径 0.2 cm のため，流れを阻害することによしており，クマイザサからの蒸発散と地面下の土壌水
り樹液流の速度は一様にならず，また，熱伝導率は分の蒸発を考慮した林床からの蒸発散量を考える必要
kx>ky で等方均質にならない，というように 4 つの仮定はどれも現実的でないことがわかった。実際，解析解がある。この蒸発散量をどのように見積もるのかが今で得た樹液流の速度を用いて日蒸散量を求めるとペン
後の課題である。
量評価の検証が必要である。

18

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)
11. Granier, A (1987): Evaluation of transpiration in a Douglasfir stand by means of sap flow measurements. Tree
Physiology 3: 309–320.
12. Granier A, Biron P, Lemoine D (2000): Water-balance,transpiration and canopy conductance in two beech stands.
Agricultural and Forest Meteorology 100: 291-308.
13. Pasqualotto G, Carraro V, Menardi R, Anfodillo T (2019):
Calibration of Granier-type (TDP) sap flow probes by a highprecision electronic potometer. Sensors 19: 2419
doi:10.3390/s19102419.
14. Miner GL, Hama JM, Kluitenberg GJ (2017): A heat-pulse
method for measuring sap flow in corn and sunflower using
3Dprinted sensor bodies and low-cost electronics.
Agricultural and Forest Meteorology 246: 86-97.
15. 熊谷朝臣（2003）樹液流計測について―グラニエセン
サーの自作法―．http://tomo-kumagai.eco.coocan.jp/granier2.pdf（2020 年 12 月 25 日閲覧）．
16. 篠原慶規・鶴田健二・久米朋宣・大槻恭一(2013)：樹液流計測法を用いた林分蒸散量の計測－森林管理による蒸散量の変化を評価するために－．日林誌 95:321-331.
17. 飯田真一・玉井幸治・野口正二・小林政広(2019)：森
林における水文過程の解明：プロットから流域スケー
ルの手法による知見と今後の展望．日本水文科学会誌
49:135–159.
18. Swanson RH, Whitefield DWA (1981): A numerical analysisof heat pulse velocity theory and practice. Journal of
Experimental Botany 32: 221-239.
19. 黒田慶子(2007)：木部樹液の動きと樹木の健康．樹木
医学研究 11: 83-88.
20. 西岡一洋(2013)：樹液流センサーの農業利用に向けた
課題と展望について．計測と制御 52: 684-689．
21. 竹内真一・高橋理一・飯田真一(2016): 樹液流動計測に基づいた樹木の移植後の生育評価－移植後4年間のタイサンボク(Magnolia grandiflora L.)の事例－．日緑工誌 42: 110-115.
22. Edwards WRN, Warwick NWM (1984): Transpiration froma kiwifruit vine as estimated by the heat pulse technique and
the Penman-Monteith equation. New Zealand Journal of
Agricultural Research 27: 537-543.
23. Barrett DJ, Hatton TJ, Ash JE, Ball MC (1995): Evaluationof the heat pulse velocity technique for measurement of sap
flow in rainforest and eucalypt forest species of south-eastern
Australia. Plant, Cell and Environment 18: 463-469.
24. Wilson KB, Hanson PJ, Mulholland PJ, Baldocchi DD,
Wullschleger SD (2001): A comparison of methods fordetermining forest evapotranspiration and its components:
sap-flow, soil water budget, eddy covariance and catchment

謝 辞
本研究では，ヒートパルス式蒸散流速計の使用にあたって北海道大学農学研究院の玉井裕教授と佐野雄三教
授に大変お世話になった。心から御礼申し上げます。
また，2 名の査読者および編集委員長からは，本稿に対し多くの貴重なご意見をいただいた。ここに厚く御礼申し上げます。
参考文献
1. Bosch JM, Hewlett, JD (1982): A review of catchment
experiments to determine the effect of vegetation
changes on water yield and evapotranspiration. Journalof Hydrology 55: 3-23.
2. Baldocchi DD, Hicks BB, Meyers TP (1988): Measuringbiosphere-atmosphere exchanges of biologically related gaseswith micrometeorologicalmethods. Ecology 69: 13311340.
3. Todd RW, Evett SR, Howell TA (2000): The Bowen
ratio-energy balance methods for estimating latent heat
flux of irrigated alfalfa evaluated in a semi-arid,
advective environment. Agricultural and Forest
Meteorology 103: 335-348.
4. Allen RG, Pereira LS, Raes D, Smith M (1998): Cropevapotranspiration-Guidelines for computing crop waterrequirements. FAO Irrigation and drainage paper 56,
Food and Agriculture Organization of the United Nations
Rome, M-56, ISBN 92-5-104219-5.
http://www.fao.org/3/X0490E/X0490E00.htm (2020 年
11 月 15 日閲覧)
5. Marshall DC (1958): Measurement of sap flow in conifers byheat transport. Plant Physiology 33: 385-396.
6. Swanson RH (1994): Significant historical developments inthermal methods for measuring sap flow in trees. Agriculturaland Forest Meteorology 72: 113-132.
7. Smith DM, Allen SJ (1996): Measurement of sap flow in plant
stems. Journal of Experimental Botany 47: 1833-1844.
8. Sakuratani T (1981): A heat balance method for measuring
water flux in the stem of intact plants. Journal of Agricultural
Meteorology 37: 9-17.
9. Tanaka T, Uchida E, Yokota T, Hagiwara A, Takeda Y (1994):
Comparison between water absorption rate and sap-flow ratemeasured using the improved stemheat-balance method. Journal of the Japanese Forestry
Society 76: 500-505.
10. Lascano RJ, Goebel TS, Booker J, Baker JT, Gitz III DC(2016): The stem heat balance method to measuretranspiration: evaluation of a new sensor. Agricultural
Sciences 7: 604-620.

19陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 3 - 20 (2021)
water balance. Agricultural and Forest Meteorology 106:
153-168.
25. Foster MA (2017): How reliable are heat pulse velocity
methods for estimating tree transpiration? Forests 8: 350
doi:10.3390/f8090350.
26. Wullschleger SD, Hanson PJ, Todd DE (2001): Transpirationfrom a multi-species deciduous forest as estimated by xylem
sap flow techniques. Forest Ecology and Management 143:
205-213.
27. Hatton JT, Moore SJ, Reece PH (1995): Estimating standtranspiration in a Eucalyptus populnea woodland with theheat pulse method: measurement errors and samplingstrategies. Tree Physiology 15: 219-227.
28. Vertessy RA, Benyon RG, O’Sullivan SK, Gribben PR(1995): Relationship between stem diameter, sapwood area,leaf area and transpiration in a young mountain ash forest.
Tree Physiology 15: 559-567.
29. 近藤純正(1994): 水環境の気象学．朝倉書店，東京，350pp．
30. 石井 孝(1998): 衛星データによる葉面積指数 LAI の推定．水文・水資源学会誌 12: 210-220.
31. Stewart JB (1988): Modelling surface conductance of pine
forest. Agricultural and Forest Meteorology 43:19-35.
32. Zhang H, Simmonds LP, Morison JIL, Payne D (1997):

Estimation of transpiration by single trees: comparison of sap
flow measurement with a combination equation.
Agricultural and Forest Meteorology 87: 155-169.
33. 玉井幸治・服部重明(1994)：落葉広葉樹林における林床面蒸発のモデル化と流域への応用．日林誌 76: 233241.
34. Sun X, Onda Y, Kato H, Otsuki K, Gomi T (2013):
Partitioning of the total evapotranspiration in a Japanesecypress plantation during the growing season. Ecohydrology
7: 1042-1053.
35. 伏谷賢美・木方洋二・岡野 健・佐道 健・則元泉・有馬孝禮・堤 壽一・平井信之(1985): 木材の物理．文栄
堂出版，286pp.
36. 石丸 優・古田裕三・杉山真樹(2017)：木材の物理．
海青社，210pp.
37．木材部・木材利用部(1982): 木材の性質一覧表．林業
試験場研究報告 319: 85-126.
38. Grainer A (1985): Une nouvelle methode pour la measure duflux de seve brute dans le tronc des arbres. Annals of Forest
Science 42: 193–200.
39. 手代木純・柳 雅之・内田大貴(2015)：ライシメータ
ー法とポロメーター法を用いた単葉からの実大樹木
蒸散量推定の可能性の検討．日緑工誌 41: 181-186.

受付：2020 年 11 月 19 日
受理：2021 年 2 月 1 日

20陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology) vol. 3: 21 - 40 (2021)

原著 (Research Article)

亜寒帯湿原における小湖沼の熱的変動- 浮島の浮沈との関係 岡田 操 1, 2

Intra- and inter-annual thermal variations of a pond in a subarctic mire
- Relations to the vertical movement of a floating islet Misao OKADA 1, 2
摘要 筆者はサロベツ湿原にある小湖沼である瞳沼（面積約 7,000 m2）で長期間にわたり，沼の凍結過程
や水温特性などの観察・計測を行なった。この沼には浮島がある。調査の過程でこの浮島が冬季には沈んでしまうことが分かった。そこで浮島の沈下過程を理解するために，絶対圧水位計を浮島に取り付け，その浮沈の様子を連続計測した。これと同時に，水深方向に等間隔の深度で水温を連続計測した。この一方で，これまでのシミュレーションモデルに改良を加え，夏冬の湖沼の境界条件を明確に区別することにより，冬期間の水温挙動と湖面の凍結との定量的なメカニズムを説明する新しいモデルを構築した。このモデルを用いて，湖沼の熱的変動を把握し，更に湖面の凍結過程を再現した。これに基づいて瞳沼の凍結過程と浮島が沈下するメカニズムとの関係を明らかにした。
キーワード 浮島，浮沈メカニズム，瞳沼，熱収支，数値シミュレーション
Abstract
The author has observed and measured the physical properties of Hitiomi-numa Pond, a small pond in
the Sarobetsu Mire for about ten years. There is a floating islet in the pond. From the survey, it was found out that thefloating islet virtually sinks every winter. In order to understand this sinking mechanism, the author fixed anabsolute-pressure water level gauge into the floating islet. Thereby, sinking depths of the islet were continuously
measured, while water temperature was monitored at some certain water depths. Meanwhile, the author modified aprototype model for simulating thermal conditions of a pond. Defining the differences between summer and winterboundary conditions, a new model was built up for quantitatively explaining the mechanisms of water temperature
distributions and the freezing processes of the pond. Numerical simulations by the model allowed us to understandthe thermal variations of the pond, since the model well simulated the freezing process of the pond surface.
Furthermore, it was revealed that there is a close relationship between the freezing process and the sinking
mechanism of the floating islet in winter.
Keywords Floating islet, Sinking mechanism, Hitomi-numa Pond, Heat balance, Numerical simulation

1. はじめに視するインターバルカメラを設置した。これらによってこの沼の浮島は冬季に沈むことが確認された。
浮島という言葉が使われる機会は少なくないが，明確な定義とともに使われることは少ない。高層湿原にある池沼である,池溏の水面上にあらわれた泥炭塊を浮島と認識して，浮島沼などと呼ばれている例は全国

筆者はサロベツ湿原に大きな浮島を浮かべる沼があることを確認し，2008 年から 2017 年まで 10 年間にわたって，調査・測量・計測を行なってきた。沼と浮島の形状を把握する測量を手始めに，浮島の動きを監*1 岡田技術士事務所 〒065-0028 札幌市東区北 28 条東 6 丁目 1-11
*2 ㈱水工リサーチ〒062-0933 札幌市豊平区平岸 3 条 7 丁目 1-27
DOI: 10.34502/physhydro.3.1_21

21

Okadax Tech. Engineering Inc., Sapporo 065-0028
Suiko research co., Ltd., Sapporo 062-0933

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)に数多くある。しかし本稿で浮島と定義するのは，泥炭地の池沼に実際に浮かんで漂っている泥炭土壌の塊りであって,上記のようなものは浮島としない。浮島は生活圏から離れた地域に存在する場合が多いため，実在する浮島の数と人間の認識との間には大きな隔たりがあるものと考えられる。そのような浮島に関連する研究は尾瀬ヶ原のように良く知られた湿原や，京都市
内にある深泥ヶ池などでの例がある。
それらの中でも，その浮沈について触れた研究例は多くはない [1, 2]。
国際的に観ても，浮遊する泥炭という特殊な環境に着目して，生態学的・化学的に研究した例はあるが [3- 5]，
物理学的観点から浮島の浮沈に触れた研究はみられない。
瞳沼の観測では，5 年余りの浮島の浮沈と，水深別 5層の連続水温の記録，さらに浮島の挙動と湖面状態を捉えた画像記録を取得した。報告されることが稀な浮
島の浮沈現象や，結氷下の水温変動という現象を，複合的・定量的に把握できた，と考えている。これらを素材として，ここで起こっている現象の詳細なメカニ
ズムを，数値モデルを用いて説明することとした。
コンピュータ・シミュレーションによる解析手法は，
素早い現象や，ゆっくりした現象の解析にその特性を発揮してきた。さらに現象のプロセスが複雑な場合などに手法の長所が利用されてきた。各要素の寄与をプ
ログラムのアルゴリズムを微修正しながら特定していく手法は，現象のメカニズムを解明・理解していく上で有効な方法である。水文・水理学の分野でもこの特性を生かした研究が 1970 年代から行なわれ，実績を積んできている。
ダム湖などにおける水温の鉛直分布は，取水水温・放流水温がこの水を利用する際に大きな影響を与える
ため，関心を持って観察・計測されてきた。この課題に関連し，湖沼における水温鉛直分布の数値解析的研
究は，ダム湖における水温・水質の環境問題に対応するために始まった感がある。古くは MIT モデルと呼ばれる数値モデルが開発され [6, 7]，この系統の研究が続いた [8, 9]。
わが国では貯水池における選択取水水温や，濁水の長期化などの課題に対応して，電力中央研究所 [10,11] や京都大学防災研究所 [12, 13] などによる研究が始められ，鉛直 2 次元の動的なモデル [14, 15] などへと引き継がれている。一方で小湖沼での水温成層の変動についての研究も進み，様々な鉛直 1 次元水温成層の予測モデルが実湖沼に適用されている [16 - 18]。さらに河川の熱収支に適用した研究もある [19]。
水温成層を取り入れた湖沼の熱的挙動の研究は，上記のように様々な研究者によって行われている。その多くは温帯湖・亜寒帯湖で行われ，凍結しない湖沼がvol. 3: 21 - 40 (2021)

対象であったり，あるいは凍結期間を対象としないな
どの理由から，大半の研究は凍結期間を欠いている。
夏季から冬季へと継続する現象は，液体の水と固体の氷との間の相変化を伴なって，熱の移動・温度変化を単一の境界条件で表わせない。そこには凝固エンタ
ルピー（凝固熱・融解熱）という，大きな熱量を介する複雑な過程が存在する。融解～凝固プロセスでは，
結果として生じる水と氷の双方とも，検証の場に存在し続ける一方で，密度が変化し容積も変化するという複雑な過程を示す。
本研究は湖沼における水温構成の変化に着目する。
更に夏季と冬季の熱収支境界条件の相違を，明確に定義することによって，湖面の凍結メカニズムを明らかにする。その上で，水温・熱収支要素・凍結氷厚などの変動と，浮島が冬季に沈没するという特異な現象との関係を，定量的に明らかにする。
2. 研究対象
2-1. サロベツ湿原
研究対象としたのはサロベツ湿原にある瞳沼という小湖沼である。サロベツ湿原は北海道の最北部，豊富町と幌延町にわたって広がっている (Fig. 1)。
辻井ほか[20] によれば，
かつて 15,000 ha 近い広さであったが，1960 年代以降の農地開発を経た現在では約6,700 ha にまで減少している。それでも高層湿原として，わが国最大の面積を有する。湿原の中には多くの河道の痕跡が観られ，激しい河道変遷を経てきた経緯がある [21]。
岡田 [22] によれば，サロベツ湿原のほぼ中央にある環境省のサロベツ湿原センターの南方 3.5 km，ペンケ
沼の北北西 1 km ほどに瞳沼は位置している。Fig. 1 は2003 年4 月と5 月に環境省がレーザースキャナー（オプテック社製 “ALTM1225”，対地高度 800 m，コース
間隔 200 m，レーザー発射回数 25,000 Hz，スキャン角度±20°，サイドラップ 66 %）で計測した DEM デー
タ（平面直角座標第 12 系 (原点北緯 44°0′0″，東経 142°15′0″) に基づいた 1 m グリッドの標高値）
により作図した傾斜量図である。傾斜量図としてのFig. 1 では，傾斜の急なところは濃く，緩やかなとこ
ろは淡色で描かれ，水平な水面は白く表わされている。
図の中央部を占める淡色の部分は，高位泥炭地植生が優占するサロベツ湿原の中核部である。
その縁辺部に，樹枝状の系統がいくつも見られるが，これらは湿地溝と呼ばれるこの湿原特有の地形である。瞳沼は中央部
やや下の部分，縁辺斜面の下部に描かれている。その南方にはペンケ沼がある。
瞳沼はその規模からいえば，湿原にごく普通に存在する小さな池沼の一つに過ぎないが，浮島を浮かべて22

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

図 2. 瞳沼の手前に接岸している浮島．
Fig. 2. Hitomi-numa Pond connected by a floating islet.

図 1. サロベツ湿原の中の瞳沼の位置．
（格子間隔は 500m）
Fig. 1. Location of Hitomi-numa Pond. (Grid interval 500 m)

いるという特徴を持つ。面積約 7,000 m2 の沼の大きさに比べて，浮島の表面積は約 2,000 m2 （上記 DEM データより計測）を占め，相対的に大きい浮島と言える（Fig. 2）
。更に後述の浮島泥炭の厚さ 1 m とから，この浮島の重量は約 2,000 トンほどと見積もられ，国内で知られている，浮いて漂う浮島の中では最大級だと考えられる [22]。
岡田 [22] によれば，Fig. 3 に示す瞳沼の長軸・短軸の長さは約 156 m，71 m で，湖岸はほぼ垂直に切り立っている。
沼の水深の計測結果は 2.05 m から 2.25 m で，ほとんど凹凸のない平らな沼底をなしている。湖岸の一部には，浮島の衝突で崩れた湖岸泥炭の堆積と思われる僅かな高まりもあるが，湖岸を含めて平坦な沼底である。沼の大きさは結氷期に巻き尺を用い，深さは標尺を水中に挿入して測定した。非結氷期には浮島上から，結氷期にはドリルで氷に穴を開けて水深を測定
し，測深位置は簡易 GPS で計測した。
Fig. 4 に示す浮島の長軸・短軸の長さは 74 m，54 m，泥炭の厚さは 0.85 m から 1.25 m と計測され，これも凹凸の少ない平板状をなしている。泥炭の厚さは手製の層厚測定器を用いて計測した。
沼の深さと浮島の厚さとの関係から，夏季に浮島はこの沼のどこにでも浮いた状態で移動できる。しかし

図 3. 浮島を除いた瞳沼の平面図と主要水深.○:結氷期調査地点,●:結氷期調査地点,◎カメラ設置位置.
Fig. 3. Plan view of Hitomi-numa Pond without the floating
islet in Fig. 2 and main water depths. ○: Measured in open
water period,●:Measured in ice covered period,◎: Camera
mounting position.

浮島の長軸が沼の短軸よりも僅かに長いため，浮島の回転は Fig. 4 の状態から時計回りに 90°あまりの範囲に限られている。
2-2. 浮島
浮島は生きている植物を乗せた状態で，泥炭という土壌の塊りが，水面に浮いているものである。瞳沼の浮島のように着底していなければ，風による剪断摩擦力によって水面を漂うものもある。泥炭は植物の遺体

23

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

図 5. 計測機器の設置方法の模式図．
Fig. 5. Schematic diagram of a way for fixing measuring
instruments.
.

福島県蓋沼の約 1,200 年前 [27]，山形県大沼の 1,300
年以上前 [23]などと比較すると，極めて新しい浮島と言える。
浮島の浮沈について，積雪期以外の報告は [1，2] などにみられ，積雪期のものには [28] などがある。

図 4. 瞳沼の浮島泥炭の厚さと計測機器設置位置．
Fig. 4. Peat thickness (m) (black squares) of the floating islet on
Hitomi-numa Pond and locations (gray marks) of measuring
instruments fixed.

3. 方法
3-1. 観測方法
瞳沼には，大きな浮島が浮かんでいるという特異な事情がある。更に，湖岸に設置したインターバルカメ
ラの画像記録と現地観察によって，この浮島が冬期間，
沼の中へ沈没してしまうことも分かった。この沈没現
象と水位の変動を定量的に把握するために 3 基の絶対圧水位計を設置した。沼底に設置した同機により沼の水位を，湖岸と浮島泥炭中の塩ビパイプ中に垂下した
同機により，湖岸と浮島の浮沈を計測した。更に水温構造の推移を計測するために，浮島に係留したフロー
トから，水中に複数の水温ロガーを垂下して，連続計測した。先に示した Fig. 4 には，機器の浮島上の平面的な設置位置を示した。また，Fig. 5 には浮島への機器の設置方法を模式的に表した。
沼底と湖岸に設置していた絶対圧水位計が遺失したため，水位の記録はわずか半年分ほどしか得られていない。その記録は，水位が大きく変動する融雪期を含んでいたにも拘わらず， 0.2 m 程度の小さな変動であった。このことから，この沼では水位が大きくは変動することが無いと判断し，水位を一定として扱った。
浮島に設置した水位計と浮島の相対位置が変化しな
ければ，浮島が沈下すると，水位計は水位が上昇した
ような記録を残す。この記録上の水位変化を，浮島の沈下量の変化と見なすことができる。浮島の浮沈現象の計測は，塩ビパイプに多数の穴を開けたストレーナが倒伏・堆積したもので，土壌中への酸素の供給が湿原の水に阻まれて，十分に行われないため，分解が進
みにくい。もともと植物遺体からできている泥炭の中には，比重が 1 よりも小さいものもある。以下のような経緯で，泥炭本体から切り離された泥炭塊の中には，浮くようになるものもある。それが浮島である。
岡田 [23，24] によれば，浮島の成因は 3 通り考えられる。一つ目は湖沼の湖岸の一部をなしていた泥炭の塊りが，本体から切り離されたり崩れ落ちて，水面に浮かぶようになったもの。二つ目は湖底の堆積物が水面に浮かび上がり，表面に植物が生えたもの [25,26]。浮かび上がるきっかけは，湖底に発生したガスによる浮力増加や，湖底に根を張った浮葉植物の葉や茎が，降雨などによる水位上昇に伴って浮きあがり，湖底泥炭を引き上げる場合などがある。三つ目は元々の傾斜地で，地滑りや土砂崩れが起こって，地上に水が溜まるようになり，そこで元々形成されていた泥炭が浮くようになったものである [27]。
瞳沼の浮島の場合，現在より更に大きかったペンケ
沼の湖岸に形成されていた泥炭塊が，泥炭地開拓の過程で下流の水位が上がり浮くようになった。その後これが泥炭本体から切り離されたという複合的な成因で，1970 年代の前半にできたものと考えられている [22]。
その形成時期は京都府深泥ヶ池の 200～1,000 年前 [2]，24

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)（内径 20 mm）を浮島泥炭に挿入し，中に絶対圧水位計（応用地質社製 “S&DL mini”，測定範囲 4 m，水位
分解能1 mm，
時計精度1 秒/日以下，
電池寿命6～8 年）を吊り下げて，毎正時の圧力を計測した（Fig. 6）。絶対圧水位計の計測値は，大気圧の変化の影響を受ける
ので，約 3.5 km 離れた環境省のサロベツ湿原センターで計測した大気圧データを用いて補正している。
深度別の水温の計測は，フロートから 50 センチ間隔
で吊り下げた，5 基の水中用温度計測データロガー（Onset 社製 “TidbiT V2”，精度±0.2℃，分解能 0.02℃，
時間精度±1 分/月，電池寿命約 5 年）を用いた。最大
深度約 2.2 m の瞳沼で，各水温ロガーを 0.5 m 間隔で，0.0, 0.5, 1.0, 1.5, 2.0 m の 5 深度に沈設し，毎正時の水温を計測記録した（Fig. 7）。
両種の機器とも数年にわたる連続観測が可能であるが，半年から一年以内の間隔で，機器の点検とデータの回収を行なった。その結果 2012 年 5 月から 2017 年10月まで約5年半にわたる連続した観測値が得られた。
この中には欠測期間が含まれている。水温ロガー全 5基のデータは2013 年12 月から2014 年4 月までの冬季
約 5 ヶ月間の欠測を含んでいる。
また，浮島の挙動や湖面の状態を観察するためにイ
ンターバルカメラ（Bicom 社製 time-lapse camera“Garden Watch Cam”， 画角 54°F=33-34 mm に相当，
記録画素 1,280×1,024 画素，撮影間隔 5 秒～12 時間に設定可能，単 3 電池×4 本，電池寿命 3～4 ヶ月）を，
Fig. 3 に示した瞳沼湖岸に設置した。湖面を自由に動き回る浮島を捉えるために，1 台の画角では湖沼全域をカバーしきれないと考え，2 台を並べて設置した（Fig. 8）
。浮島は風向き次第で，弱い風でもその剪断応力を受けて動き出し，途中で止まることなく対岸に着岸することを繰り返している。従って浮島の冬期間の位置は，湖面が凍結する日の風向で決まる。湖面が凍結する時期の気象条件から，浮島は瞳沼の南半分で冬季固定する年が多い。これを踏まえて，カメラは南方と南西方向へ向けて取り付けた (Fig.3, Fig.8)。
一方，数値シミュレーションを実行するためには，
種々の気象観測データを必要とする。最も近い気象観
測所はアメダス豊富地点（地域気象観測システム：北緯 45 度 6 分 6 秒，東経 141 度 46 分 48 秒，標高 14 m，1991～2020 年の統計期間では年平均気温 6.3℃，年降水量 1030.8 mm，降雪の深さ 703 cm）であるが，必要な観測値の全てが観測されていない。不足分は稚内地
方気象台（北緯 45 度 24 分 55 秒，東経 141 度 40 分 43秒，
標高2.8 m，
年平均気温7.0℃，
年降水量1109.2 mm，
降雪の深さ 477 cm）の観測値で補うこととした。瞳沼（北緯 45 度 4 分 52 秒，東経 141 度 41 分 46 秒，標高2.1 m）からこれらの気象官署までの距離は，アメダス

vol. 3: 21 - 40 (2021)

図 6. 浮島に挿入した塩ビパイプ内に絶対圧水位計を吊り下げる様子．
Fig. 6. Absolute-pressure water level gauge “S&DL mini”
is hung in the PVC pipe.

図 7. 水温データロガーを吊り下げた浮き．
Fig. 7. Temperature data loggers “Tidbit V2” are hung froma moored float.

図 8. 湖岸に設置したインターバルカメラ．
Fig. 8. Two interval cameras installed on the shore of pond.

豊富地点までが約 7 km，稚内地方気象台までが約 37
km である。気象観測値の内，降水量 (mm h-1)・気温(℃)・風速 (m s-1)・日照時間 (h) はアメダス豊富地点のものを用い，湿度(%)・全天日射量 (MJ m-2 h-1)・雲
量は稚内の観測値を用いた。
3-2. 観測結果
Fig. 9 は水温変動の観測結果である。左図は 6 年間

25

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

図 9. 瞳沼の実測水温変化 (毎日 14 時の計測値を使用)．左：鉛直水温分布の変動 (白線は浮島の実測沈没量)，右：深度別
水温変動．
Fig. 9. Observed water temperature in Hitomi-numa Pond (First season). Left: Vertical water temperature distribution. The white line
indicates observed sinking depths of the floating islet. Right: Water temperature at five depths.

のうち，最初の一シーズン分の水温鉛直分布の変動を，色分けで示した図である。縦軸は下に向かって水深方向を示し，最大 2.2 m の深さである。水温ロガーを 0.5m 間隔で，
0.0, 0.5, 1.0, 1.5, 2.0 m の 5 つの深度に設置して計測した（Fig. 7）値を用いて，0℃から 25℃までのグラデーション表示している。1 時間ごとの観測値をデータとして取得しているが，煩雑さを避けるため，
14 時のみの値を用いて描いた図である。
Fig. 9 の左図で，水温の低い青で表示されている期間が冬季であり，その中に白い線で浮島の沈下量を示している。
本稿では秋の循環から春の循環までの期間，
即ち湖水の水温成層が逆列をなしている期間を冬季と定義する。Fig. 9 の右図は同じ期間，5 つの深度で計測した深度別水温の変動を示した。この図も煩雑を避けるため毎日の 14 時の観測値のみを表示している。
毎
年表層水温は 7～8 月ころ，30℃前後の最高水温に達し，一方底層水温は少し遅れて 8～9 月ころ，
15℃あまりにまで上昇している。
Fig. 10 は，Fig. 9 の期間を含む浮島の沈下量（黒い実線）の変動を示している。この図には，他に約 7 km離れたアメダス豊富地点の日降雪量（紺色棒グラフ）
・
積雪深（緑実線）
・累加降雪量（赤実線）をも示している。日降雪量は毎日の降雪量を、他は 1 時間ごとの観測値を表わしている。
これらの図を見ると，浮島は 12 月ころ沈み始め，後戻りすることなくほとんど一方的に沈み続けていくが，
1 月末から 2 月にかけて沈下の勢いを失い，4 月に入ると急激に沈下量を減じ浮上する。
この浮島の沈下が始まるのと同じ時期に雪が降り始
め，積雪荷重と浮島沈下との間に関係があることを推察させる。しかし積雪深（緑実線）の変動は浮島の沈下とズレが多く，両者の関連はあまり良くない。積雪図 10. 浮島の沈下量 (黒) と累加降雪量 (赤) の変化．
Fig. 10. Temporal variations of daily snowfall, sinking depth, cumulative snowfall and snow depth.

26

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

別なく凍りつづけている状態と解釈できる。氷の浮力と積雪荷重がバランスを取っている状態が読み取れる。
3-3. 熱環境変動モデル
複数年度の測定結果には，湖面が凍結して雪が積もるなどの周期的に起こる現象と，降雪時期が早かった
り結氷が遅れるなど，各年の特別な条件によって引き起こされる現象とが混じり合って，複雑な挙動を見せる。計測された詳細で複雑な記録を理解するために，普遍的なパラメータを用いた説明は有効な手段となる。
そこで計測結果を定量的な数値解析モデルを用いて，現象へ寄与する要素と境界条件となるパラメータを同定し，現象を再現することにした。湖沼周辺で生起している現象への各要素の寄与を，プログラムのアルゴ
リズムを微修正しながら特定していくことは，現象のメカニズムを理解する上で有効な方法である。
開水面面積があまり大きくなく，また水深があまり深くない瞳沼のような湖沼を対象として，熱収支のモデルを選定するに際して，湖水の流動を含むモデルだと，計算に不安定性が入り込む危惧がある。例えば狭い対象領域に多数のグリッドを設定すると，グリッド間隔が小さくなり，計算が不安定になる。これを防ぐために，安定条件を満たすよう時間差分を小さく取ると，計算に要する時間が増大する。各層が一様であると見做せるような小さな水域では，水平方向の熱の移動を考慮しても意味がない。
そこで本稿では，小規模の湖沼に適用するモデルとして湖水位は変動しないものとした。また風の剪断応力としては，実際の沼水が流動することまでは考慮せず，渦動拡散係数への寄与までを取り込む静的な多層モデルを用いることとした。
ここで用いるモデルは，𝑧 方向にのみ複数の水の層を設定する。したがって各層は，厚さと面積をもった
1 枚のプレートのようなイメージとなる。また水の流動を考慮しないので，計算の安定のための時間差分に配慮する必要が無い。本稿で用いたモデルでは，デー
タ取得間隔 (1 時間) を 2 等分ないし 10 等分した時間を計算間隔とした。熱移動に関する水温の連続方程式
は次のように表わされる。

図 11. 瞳沼の凍結湖面に開けた穴の状態．
Fig. 11. Situation at or around a hole bored on the frozen pond
surface.

密度は積雪の経緯・気象条件などで変化するため，積雪深の変化は荷重の変化を正確に反映していないと考えられる。これに対して，赤実線の累加降雪量の変化と黒実線の浮島の沈下量は，ある時期までほとんど同じ変化を続けていく。Fig. 10 に見られる浮島の沈下量と累加降雪量による荷重の変動の関係，及びその関係が途中から劇的に離れていく現象については解析計算
後に述べる。
Fig. 10 に赤い実線で示したのは，毎日の降雪量（紺色棒グラフ）を積算した累加降雪量に 0.2 を乗じた値である。ここに 0.2 という値は，新雪の比重と見做して用いたもので，積雪初期の積雪密度 [29] や，積雪の上部の密度 [30] から裏付けられた値である。
毛利 [31] によれば，結氷の最初は比較的透明な核氷となり，積雪が無ければ氷は下方へ成長して，核氷のみの層となる。ここで”核氷”とは，水がそのまま凍結して透明氷となったものの意味で使われている。こ
こに降雪があると氷面上に雪と水の混合層ができて凍る。水と雪の混合層が凍結すると，白濁不透明な氷となり，毛利 [31] はこれを”乳氷”と呼んでいる。積雪量が更に多くなると，その荷重によって氷が押し下げられ亀裂が生じ，これに沿って沼水が押し上げられ，氷と積雪の間に浸みだす。この部分の熱収支がマイナスの状態であれば，これも凍結して氷の厚さを増す。
この部分を毛利 [31] は”雪氷”と呼んでいる。
Fig. 11 はこの時期の調査時の写真である。瞳沼の凍結湖面にドリルで開けた穴からは，沼水が噴き出すよ
うなことはなく，被圧されていないことが分かる。また透明ではない氷の状態を見ると，毛利 [31] のいう”雪氷”に相当し，氷面亀裂から浸みだした水が積雪とともに凍結したものと考えられる。このドリル穴の位置
は浮島の真上であり，ドリルの刃に絡んで引き上げられた赤い木の実が穴の左に見える。写真に写った氷面と水面との間には大きな差はなく，湖面も浮島上も区

𝜕𝑇
𝜕𝑇
𝜕𝑇
𝜕2𝑇
+𝑢
+𝑤
= 𝐷𝑀 2
𝜕𝑡
𝜕𝑥
𝜕𝑧
𝜕𝑧
𝜕
𝜕𝑇
1
𝜕𝐻
+ (𝐾𝑧 ) +𝜕𝑧𝜕𝑧
𝜌𝑤 𝐶𝑤 𝜕𝑧

… … (1)

左辺第 1 項は加速度項，2 項は水平移流項，3 項は鉛直移流項，右辺第 1 項は分子拡散項，第 2 項は渦動拡
散項，第 3 項は熱収支項である。

27

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

ここに 𝑇：水温 (℃) ，𝑢：水平 𝑥方向流速 (m s-1),
𝑤：鉛直 𝑧方向流速 (m s-1), 𝐷𝑀：分子拡散係数 (m2 s-1)，
𝐾𝑧：渦動拡散係数 (m2 s-1)，ρw：水の密度 (kg m-3)，𝐶𝑤：
水の比熱 (4.186 kJ kg-1 K-1 )，𝐻：単位面積単位時間当たりの収支熱量 (kW m-2)である。
1) 移流項左辺第 2 項，水平移流項は河川などからの流入や流出のある水域の場合に考慮する必要がある。多層モデ
ルの各層の水は，水温・濁度とそれらに基づいた密度という属性を持っている。河川から流入する水も，同様な属性を持っているので，流入水はそれと同じ密度の層に貫入することになる。静的モデルゆえ動く水の慣性は考慮されない。一方，流出の際も表層水の流出に伴なって，連行される形で上下に隣接する層からも流出水が加わる。これらの流入・流出には，水平移流が伴なう。
水平移流による水の移動を補償するかたちで，左辺
第 3 項，鉛直移流が起こる。したがって水域の容量に比べて大きな流入・流出がなければ，これらの二つの項は無視されることになる。また，水面を通じての日
射や潜熱・顕熱交換などにともなう対流も，鉛直移流という形では考慮されない。
モデルとしては両移流項を考慮しプログラムコードを組んであるが，本研究で対象とした瞳沼では，地表水の流入・流出がないので境界条件として両移流項の寄与はない。

図 12 沼の熱収支模式図．(静的多層モデル)
Fig. 12. Schematic diagram of the heat balance of a pond (static
multi layer model).

𝑆𝑡 = 𝜕𝜌𝑤 ⁄𝜕𝑧で表わし，水温躍層内で最大となる。躍層以深の安定度はこの最大値を用いる。 ρw ：水深 𝑧 における水の密度 (kg m-3 ) で，厳密な値は Chen ほか [33] によった。
一方風の効果 𝐸𝑤 は，∞

𝐸𝑤 = ∑(𝑊𝑑2⁄2𝑑+1 )𝐾𝑧 = 𝛼0

安定度による抑制効果= 𝛼0…… ( 4 )𝑑=0とした。ここで 𝑊𝑑 ： 𝑑 日前の日平均風速 (m s-1)である。ここでは，(4) 式にあるように一日経過するごとに風の影響は半減していくものとした。
分子拡散と渦動拡散を合わせて考慮している例もあるが [34]，本稿では別個に扱う。前者は後者に比べ微小であるが，本モデル中では冬期間の結氷下では分子拡散のみが考慮されている。

2) 拡散項
この静的モデルでは，各層の水が実際に動いて互いに拡散することは考慮していない。しかし，水域の水温構造を決めるのに，水面上を吹く風の影響が大きい
ため，これによる効果を取り入れる。更に，風の影響
は風が吹いている間だけでなく，吹送流となって継続する。
一方で水域では水中に水温躍層が形成されると，これによる安定度が水温伝達のバリアーとして働き，その効果も大きい。ここでは水深 𝑧 における渦動拡
散係数 𝐾𝑧 (m2 s-1) を，風速の履歴と水温構造で決まる安定度の関数として，次の経験式で求めることとした。基本的には風の効果 𝐸𝑤 は拡散現象の促進効果として，安定度 𝑆𝑡 は抑制効果として働く [32]。
吹送流による促進効果………… ( 3 )3) 熱収支項
多層モデルで表わされる沼における，熱の授受のイ
メージは Fig. 12 のように表わせる。沼に対する熱収支式は，沼の単位面積当たりの貯熱量変化 (kW m-2 )として次のように表せる。
𝐻 = 𝑅0 + 𝑅𝐷𝑛 − 𝑅𝑈𝑝 ± 𝐻𝐸𝑣 ± 𝐻𝐴𝑡±𝐻𝑃𝑟 + 𝐻𝑖𝑛 − 𝐻𝑜𝑢𝑡 + 𝐻𝐺𝑛𝑑

√𝐸𝑤
√𝑆𝑡

… … (5)

ここに𝑅0：太陽からの短波放射量のうち水中に進入する部分， 𝑅𝐷𝑛 ：大気・雲からの下向き長波放射量，
𝑅𝑈𝑝 ：湖水からの上向き長波放射量，𝐻𝐸𝑣 ：蒸発・凝結に伴なう潜熱流束， 𝐻𝐴𝑡 ：大気との間の顕熱流束，… (2)ここで安定度 𝑆𝑡 は28

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

𝐻𝑃𝑟 ：降雨・降雪によって持ち込まれる降水熱流束，
𝐻𝑖𝑛 ：河川などからの流入に伴ない持ち込まれる熱流束，𝐻𝑜𝑢𝑡：河川などへの流出に伴ない持ち出される熱流束，𝐻𝐺𝑛𝑑：湖底から供給される地熱流束である。各
項の符号は水塊に加わる量を + とし，奪われる量を－ としている。環境条件次第で入れ替わり得る項は± とした。
通常下向きの長波放射量 RDn と上向きの RUp とは，両者を併せて長波収支として扱われることが多いが，両者の挙動は夏季と冬季とで異なるので本稿では両者を別々に取扱う。同様に潜熱流束 HEv と顕熱流束 HAt
も冬期間の境界条件が異なるので別々に扱う。

m-2 K-4 )=4.90×10-6 (kJ m-2 day-1 K-4 )， 𝑇𝑤𝑘 ：絶対温度で表わした表面水温 (K) である。
一方，上空から水面に入射する長波放射量 𝑅𝐷𝑛(kW m-2) は Swinbank の式で表わされる [37]。
𝑅𝐷𝑛 = 0.97 𝜎
6 }{0.937 × 10−5 (1.0 + 0.17 𝐶𝑙𝑑2 ) 𝑇𝑎𝑘……… (10)
ここで𝐶𝑙𝑑：雲量 (0～1)， 𝑇𝑎𝑘 ：絶対温度で表わした気温 (K) である。
・潜熱流束
水面における蒸発量 𝐸𝑉 (m s-1 ) はRohwerの式で表される [38]。

・短波放射
水中に進入する短波放射量 𝑅0 (kW m-2 )は，水面で太陽からの放射量 𝑆 が，一定の反射率 𝛼（Albedo アルベド）で反射されて (6) 式で表わされる。

𝑅0 = (1 − 𝛼) 𝑆

𝐸𝑉 = (0.231 × 10−3 + 0.139 × 10−3 𝑊)(𝑒𝑤 − ℎ𝑟 𝑒𝑎 ) … … (11)
ここで𝑊：風速 (m s-1)， 𝑒𝑤 ：表面水温に対応する
飽和蒸気圧 (hPa)， 𝑒𝑎 ：気温に対応する飽和蒸気圧(hPa)， ℎ𝑟 ：湿度 (0～1) である。
蒸発量 𝐸𝑉 をもとに蒸発・凝結に伴なう潜熱流束
𝐻𝐸𝑣 (kW m-2) は，… … … (6)

ここに𝑆：太陽からの放射量，𝛼：水面の短波反射率（Albedo アルベド）。水中に入った光エネルギーの多
くは水面近くで吸収され，残りは指数関数的に減衰しながら水中を伝わる。水深 𝑧における光エネルギー𝑅𝑧(kW m-2) は Lambert-Beer の法則によって次式で表わされる。
𝑅𝑧 = (1 − 𝛽) 𝑅0 exp(−𝜂 𝑧)

𝐻𝐸𝑣 = 𝜌𝑤 𝐿𝑉 𝐸𝑉

ここに 𝛽：水面吸収率，𝜂：減衰係数である。更に水深 𝑧1 と 𝑧2 に挟まれた水の層に吸収される放射量𝑅1−2 は次式で表わされる。

𝐿𝑉 = 2500.0 − 2.39 𝑇𝑤𝑐… … … (13)

これらを整理すると蒸発・凝結に伴なう潜熱流束
𝐻𝐸𝑣 は，… … … (8 )𝐻𝐸𝑣 = 𝜌𝑤 𝐿𝑉 (0.231 × 10−3 + 0.139
× 10−3 𝑊)(𝑒𝑤 − ℎ𝑟 𝑒𝑎 )… … (14)

アルベド 𝛼 は Bdyko [35] や Brandt ほか [36] などを参考に夏季の値を 0.05 と与え，また水中減衰係数
𝜂 は大八木ほか [32] などを参考に 2.0 (一定)とした。

・顕熱流束
一方，顕熱輸送と水蒸気輸送に対する拡散係数が等しいものとすると，大気との間の伝導顕熱流束 𝐻𝐴𝑡(kW m-2) も Rohwer の式で表される [38]。

・長波放射
水面での長波放射の収支は以下の通りである。水面から上空に向けた長波放射量 𝑅𝑈𝑝 (kW m-2) は StefanBoltzmann の法則に従う。
4
𝑅𝑈𝑝 = 0.97 𝜎 𝑇𝑤𝑘… … … (12)ここで 𝜌𝑤：水の密度 (kg m-3 )，𝐿𝑉：蒸発の潜熱 (kJ
kg-1 )である。水温に従う蒸発の潜熱 𝐿𝑉 は水面の温度
𝑇𝑤𝑐 (℃) を用いて次式で表わせる [19]。

… … … (7 )𝑅1−2 = (1 − 𝛽) 𝑅0{exp(−𝜂𝑧1 ) − exp(−𝜂𝑧2 )}vol. 3: 21 - 40 (2021)

𝐻𝐴𝑡 =… … (9)

ここで𝜎 ： Stefan- Boltzmann 定数=5.670×10-8 (kW1500 (𝑇𝑤𝑐 − 𝑇𝑎𝑐 )⋅⋅𝐻
𝐿𝑉 (𝑒𝑤 − ℎ𝑟 ⋅ 𝑒𝑎 ) 𝐸𝑣

… (15)

ここで 𝑇𝑤𝑐：表面水温 (℃)， 𝑇𝑎𝑐：気温 (℃) であ29

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

る。
また表面水温に対応する飽和蒸気圧 𝑒𝑤 (hPa)と気温に対応する飽和蒸気圧 𝑒𝑎 (hPa) は2
𝑒𝑤 = 6.11 + 0.4726 𝑇𝑤𝑐 + 0.009 𝑇𝑤𝑐
3+0.0005084 𝑇𝑤𝑐… … (16)2
𝑒𝑎 = 6.11 + 0.4726 𝑇𝑎𝑐 + 0.009 𝑇𝑎𝑐
3+0.0005084 𝑇𝑎𝑐… … (17)

図 13. 沈みゆく瞳沼の浮島（インターバルカメラ画像），
左：湖面の薄氷とその先の浮島 (12 月 8 日)，右：左写真と同じ位置で完全に沈没した浮島 (3 月 1 日)．
Fig. 13. Floating islet sinking in Himomi-numa Pond. (Photos
by an interval camera) Left: Floating islet and Hitomi-numa
Pond of thin ice cover (8-Dec 2013), Right: Islet sunk perfectly(1-Mar 2014).

で表わされる。
・降水熱流束降雨・降雪とともに持ち込まれる熱量 𝐻𝑃𝑟 (kW m-2)
は降雨・降雪の温度を露点温度と見做して次式とした。
𝐻𝑃𝑟 = 𝜌𝑤 𝐶𝑤 (𝑇𝑑𝑒𝑤 − 𝑇𝑤𝑐 ) 𝑃

う重要なチェックポイントを欠いていた。そのためこれまでは冬期間を外した解析期間とならざるを得なかった。
Fig. 13 はインターバルカメラで撮影した浮島の冬季の画像である。左写真では手前の瞳沼湖面の向こう側に浮島が見える。湖面には薄氷が張り始めており，
浮島はこれが融けるまでその位置を変えることができない。右の写真には同じ位置で沈み切り，積雪の上に枯れたヨシの茎が見えている。絶対圧水位計の記録ではこの状態の沈下量は 1.015 m となっている。浮島による遮蔽の状態が変わって，遠方の河畔林の見え方が異なることで判断できる。
基本的には湖沼において前述の熱収支が存在するが，自然界の湖沼では，その境界条件は夏と冬とでは大きな相違を生じる。その相違は水面の凍結・融解という
熱的に不連続な現象に関連して生じるので，湖沼の冬
季結氷の過程を取り扱うには，二つの期間の境界条件の相異や変化の過程を明確にさせる必要がある。

… … (18)

ここで ρw：水の密度 (kg m-3 )，𝐶𝑤：比熱 (kJ kg-1 K-1 )，𝑇𝑑𝑒𝑤 ：気温・相対湿度に対応する露点温度 (℃)，𝑇𝑤𝑐：
表面水温または氷温 (℃)， 𝑃 ：降水量 (m s-1 )である。
・その他
湖底から供給される地熱流束 𝐻𝐺𝑛𝑑 は一般的に他の要素に比べると極めて小さく，計算の趨勢に影響を与えることはない。またこの地域での計測例も知られていないことから考慮しないこととした。
本稿の検証で用いている静的モデルでは，対流混合を水の流動として扱っていないが，これに代わって，
正列か，逆列成層かを決めるアルゴリズムが加えられている。つまり密度が最大となる 4℃の水に比べて，より高い水温の水と，より低い水温の水が同時には存在しない，という処理を行なっている。
湖沼それぞれの形状や深さ・周囲の状況など，湖沼が個別的に持つ特徴については，これらを初期条件として個別に与えるが，それ以外の客観的境界条件，例
えば夏季と冬季とで異なるアルベドの値や分子拡散係
数値などのパラメータについては同一湖沼では値を変えず，また年によっても変えず同一の値を用いることとし，客観性を保つよう配慮した。

・拡散現象
モデルでは，湖面上を吹く風により投入されたエネ
ルギーによる拡散が考慮されている。夏季には水面を通じて入る風に起因する渦動拡散 ((2)式) によって熱が上下に拡散する。しかし冬季，湖面が凍結すると，
このエネルギーが水中に供給されなくなるため，分子拡散のみが働くことになる。

3-4. 水面凍結により変化する境界条件
シミュレーションは，プログラムのアルゴリズムの修正・適正化の過程で，対象現象のプロセスを見極めるという重要な手法と見做せる。著者は本研究に先行するモデルを使って熱収支の検討を行なってきた [39- 41]。冬期間の凍結プロセスについては氷の凍結層厚の変化をとらえたデータがなく，実現象との照合とい・日射の吸収と透過
水面からの短波放射 (日射) エネルギーは，夏冬の間で最も顕著な差異が生じる。夏季には水面から大半の日射エネルギーが水中に入っていく。しかし冬季水面が凍結しさらに積雪に覆われるようになると，より30

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

多くの日射エネルギーが反射されるようになる。ほとんどの日射エネルギーが吸収される夏季と，半分ほどのエネルギーが反射されてしまう冬季との差を，的確に与える必要がある。また後述するように昇華蒸発に消費されるエネルギーは，吸収された日射エネルギーの一部との関連が指摘されており [42]，この関連にも配慮してアルベド値を調整する。
水面または氷面・積雪面から水中に入った日射エネ
ルギーは Lambert-Beer の法則に従って，徐々に減衰しつつ下層に透過していく ((7)式)。

vol. 3: 21 - 40 (2021)

り，露点温度と表面温度との差と，降水量に応じた熱量が水域に供給される。
地熱は他の項目に比べ極めて小さく，計算結果にほとんど影響はないと考え，考慮しないこととした。
・表層水の凍結プロセス
本稿で凍結層厚として扱うのは，その厚さの水が凍結したという意味であり，氷に変化したときの容積変化を表わしていない。0 ℃の氷の密度は約 0.907 で，同じ量の水に比べて容積は 1.09 倍ほど増える。しかしその氷が水の上に浮いている限り，氷の底面の位置
はそれが水であった時と変わらない。
本モデルは水深方向に，複数の水の層を設定した水平多層モデルである。各層における熱収支を計算した
結果，密度構成にアンバランスが生じた場合にはこれを調整する。水面を含む最上層が次第に冷えて 4℃以下になっていくとき，熱の収支ばかりでなく，密度構成に矛盾が生じないような調整が必要である。つまり湖水が循環を起こす時期に，最大密度となる 4℃よりも高い水温の水と，低い水温の水が共存することがな
いような調整を行なう（巻末のフローチャート参照）。
密度構成に矛盾が生じない状態で，未凍結の最も上の層について熱収支を計算し，その総量が氷点を基準としてマイナスに転じたとき，そのマイナス熱量相当
分の厚さの水が氷に変化するものとした。逆に凍結した層での収支がプラスとなった場合に融解が起こる。

・湖面における長波放射
湖面における上空からの長波放射は Swinbank [37]により詳細に与えられ，(10) 式では雲量と気温がパラ
メータとなるから，夏・冬を問わず同じ式で求められる。
しかし，上空に向かう長波放射は，夏季と冬季に起こるプロセスが異なるので，別個に算出考慮しなければならない。上空に向かう長波放射は StefanBoltzmann 法則で求められ，表面温度のみがパラメータとなる。夏季には表層水温がこれに当たるが，冬季の結氷湖面では氷面または雪面の温度が必要になる。
小嶋ほか [43] によれば，積雪表面の乾雪の雪温は気温変化に追随し，また水に浮いた氷の温度は時間によらず氷点であるとしている。そこで，湖面が凍結しているという条件下では，表面温度として氷点の 0 ℃，
または上記のように雪面温度が追随する気温の，いず
れか低いほうの値を用いることとする。

𝛥𝐼 = −・潜熱・顕熱流束(12) 式の潜熱流束と (15) 式の顕熱流束は，Rohwerの式によりセットで求められる。水面が凍結する冬季にも顕熱流束を求めるためには，表面温度として前述の上向き長波放射と同様，氷点または気温のいずれか
低い方を用いることで，この式を準用することとした。
湖面が凍結して蒸発がなくなれば，蒸発潜熱流束はゼロとなるが，水面からの蒸発がなくなっても氷面・雪面からの昇華蒸発現象は続いていく。
東海林 [42] によれば昇華蒸発量は，短波放射の吸収量と風速との関数で表わされる。そこで冬季の結氷期間の昇華蒸発量は，日射の短波放射の吸収量の一部が，昇華蒸発熱流束に変わるとの考え方で求める。すなわち，氷面・雪面のアルベドの値を調整して，試行錯誤的に 0.40 を与えることで表現することとした。

𝐻𝑀
𝜌𝑤 𝐶𝑀

… … … (19 )ここに Δ𝐼 ：凍結厚さの増減 ( m s-1)，𝐻𝑀 ：氷点を基準とした熱量の余剰 (kW m-2 )， 𝜌𝑤 ：水の密度 (kgm-3 )，𝐶𝑀 ：融解熱 (=333.6 kJ kg-1 ) である。
既に凍結した氷の温度は 0℃ のままとし [43]，未凍結の最上層では水も氷も 0℃で，ここにおける熱収支を考える。
4. 計算結果
4-1. 初期条件と境界条件
瞳沼は最深部で約2.2メートルの深さがある (Fig. 3)
ので，深度 0.1 メートルごとに 22 層を設定し，沼の形状から各層の断面積を初期条件として与えた。瞳沼には明瞭な接続河川はなく，流入・流出ともに移流はゼロの設定とした。
境界条件となるパラメータは，それらの値を年毎の傾向に合わせて操作するのではなく，全期間を通して一括して定め，値に客観性と汎用性を持たせることと・その他の熱収支項
降水からの熱流束は，夏冬の差は設けずに扱った。
気温・湿度から決まる露点温度次第で雪にも雨にもな

31

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

図 14. 瞳沼水温の計算結果 (毎日 14 時の計算値)．左：鉛直水温分布の変動 (白線は計算凍結層厚)，右：深度別水温変動．
Fig. 14. Calculated water temperature in Hitomi-numa Pond (First season). Left: Temporal variation of vertical water temperature
distribution. The white line is the calculated thickness of the frozen layer). Right: Temporal variation of water temperature at five
depths.

した。ここで瞳沼の 5 年余りにわたる計算に用いたパラメータの，具体的な値は Table 1 のとおりである。
日射に対する地表面の反射能であるアルベド ((6)式の α) は，夏季については水面の一般的な値 0.05 を用いた [35，36]。水面が凍結した冬期間は，氷表面または雪面からの昇華蒸発熱流束 [42] にも配慮して，上記のように 0.40 を用いた。日射の水面吸収が起こる
範囲としての水面厚は 0.1 m とし，水面吸収率 ((7)式の β) は試行錯誤的に 0.90 を用いた。水中に伝わる
光エネルギーの分布を決めることになる，水中減衰率((7)式の η) は大八木ほか [32] を参考に 2.0 とした。
分子拡散係数値 ((1)式のDM ) は一般的な値7.0 ×10-7とした。渦動拡散係数を算出する (2) 式に用いる α0の値は 1.0 ×10-9 とし，(3) 式で表わされる安定度 𝑆𝑡の下限値は 5.0 ×10-4 とした。

ン分だけ Fig. 14 の左図に，同じく右図には 5 つの深度の計算水温の変動を，深度ごとに色分けして示した。
同じ期間の実測値は先に示した Fig. 9 の左右の図である。計算結果は１時間ごとに求められているが，見やすくするため実測の図と同様に，
毎日の 14 時の計算
値だけを用いて描いた。
計算結果の全般的な傾向として言えることは，浮島
沈下量の実測値に比べると，凍結した氷の厚さはそれよりやや小さく計算され，凍結期間は少し長めな結果となった。どちらか一方を実測に近付けるようにパラメータの操作を行なうと，もう一方の誤差がより顕著になってしまう。
Fig. 14 に示した二つの図は Fig. 9 に示した実測図と良く似た形状となっている。実測した 5 つの各深度(0.0, 0.5, 1.0, 1.5, 2.0 m) の水温を，計算水温とのあいだで客観評価を行なった結果が Table 2 である。両者の相関係数は大半が 0.9 以上となって，強い正の相関があることを示しており，回帰式の係数も 1.0 に近い値となっている。これらから判断して，Table 1 に設定したパラメータの組み合わせは妥当なものと考えられる。

4-2. 計算結果
・水温の変動計算結果の水温鉛直分布の変動を，最初の一シーズ

表 1. 瞳沼の水温変動に対するシミュレーションに用いたパラメータの値．
Table 1. Parameters applied to the thermal simulation of
Hitomi-numa Pond.

表 2. 各深度における計算と実測水温の誤差の指標
Table 2. Correlation coefficients and regressive formulas between
calculated and observed water temperatures at five depths.

32

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

図 15. 熱収支要素の変動（水温は 14 時値，気温と結氷層厚は日平均，それ以外は日総量）．
Fig. 15. Temporal variations of each factor in the heat balance.

ラスとなる日が増え，0℃のままである凍結湖面への顕熱流束はプラスとなって急増する。これは湖面が解氷するとともにマイナスに移行する。
これらの総計である全熱収支 (緑実線) は，盛夏を過ぎるとマイナス傾向に転じ，-0.5 kJ cm-2 day-1 前後で推移し，これが厳冬期まで継続する。全熱収支は氷の凍結層厚が最大となる，2 月から 3 月にかけてプラスに転じ，4 月には年間の最大値となる。
全熱収支の 4 月における最大化の原因は，直接的な日射量の増加と，上記のような，気温上昇に伴なう顕熱流束の増加の寄与が大きい。融雪期におけるこの顕熱流束の増加は，水面がまだ凍結しており，その温度が 0 ℃のままであるため，気温との温度差が大きくなることによる。これと同時に結氷湖面からは積極的な蒸発が起こらず，蒸発熱が奪われないことも影響している。つまり湖面がまだ凍結しているという状況が，湖面解氷の原動力となっている。結氷期に比べて，急激に進む解氷は全熱収支が最大化することが原因である。
完全解氷とともに顕熱流束はプラスからマイナスに転じ，マイナス要素の蒸発潜熱流束も加わって，熱収支は一気にゼロ近くにまで減少していく。瞳沼全体としての熱的挙動は，5 月頃から 6 月頃まで僅かなプラ
ス傾向の受熱期がつづき，8 月頃からマイナスに転じて放熱期に入り，沼の熱がゆっくりと失われていく。
全期間を通した鉛直水温分布の変動状況を，実測値と計算結果を対比して Fig. 16 に示した。この図の上
図 (a) は Fig. 9 の左図を，下図 (b) は Fig. 14 の左図を全期間に拡張したものである。Fig. 16 の上図 (a) で，浮島沈下量を表わしている白い線が，冬季以外の期間

・熱収支要素の変動
Fig. 15 は熱収支の各要素の季節変動を図示したものである。上の図には気象条件である実測値の日平均
気温 (℃ :赤実線) と日降水量 (mm day-1 :紺色棒グラフ)，及び計算結果である蒸発量 (mm day-1 :緑棒グラフ)・凍結厚さ (m :黒棒グラフ)・14 時の計算表層水温(℃ :水色実線) も併せて示した。
Fig. 15 の下図では，実測日射量 (オレンジ点線) 以
外は計算値で，それらの単位は全て (kJ cm-2 day-1 ) とした。日射吸収 (赤紫実線) は日射量の 95 %となる夏
季と，半減する結氷季とが顕著に対比される。下向き
長波 (黒点線) は夏季には日量 3.0 kJ cm-2 day-1 を超える日が多いが，冬季には半分近くにまで減少する。上向き長波 (茶色点線) は夏季には，水温が上がって長波放射量が増えるため，-3.0 kJ cm-2 day-1 以下となる日もあるが，冬季には表面が冷やされるため-2.5 kJ cm-2
day-1 程度となる。両者の和を長波収支 (茶色実線) として示したが，年間を通し夏冬に拘わらず-0.5 kJ cm-2
day-1 程度で推移し，マイナス超過で継続している。
潜熱流束 (青点線) は夏季には-1.0 kJ cm-2 day-1 程度で推移して，秋になると-0.5 kJ cm-2 day-1 前後となり，水面が凍結するとゼロとなる。凍結期間中にも昇華蒸
発・凝結が起こっているが，この分はアルベドの値に反映させているので，日射量と日射吸収量との差に含まれていることになり，グラフには表わされていない。
顕熱流束 (赤点線) は，表層水温が気温を上回ることが多い夏季には，-0.5 kJ cm-2 day-1 近くなることもあるが，雪面と気温との差が小さくなる冬季には目立たなくなる。更に時間が経過し，3 月に入ると気温がプ

33

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

図 16. 鉛直水温分布の変動比較．(a)：14 時の実測値 (白線は浮島の実測沈下量)，(b)：14 時の計算結果 (白線は計算結氷厚)．
Fig. 16. Temporal variations of (a) observed and (b) calculated water temperature distributions. White lines in (a) and (b) indicate
observed sinking depths of the floating islet and calculated thicknesses of the frozen layer, respectively.

図 17. 熱収支要素の変動 (全期間：凡例は図 15 と同じ)．
Fig. 17. Temporal variations of heat balance factors (Whole period: Legend is the same as in Fig. 15).

34

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)にも僅かずつ下がっている。これは夏季の植生成長による植物体の重量増加と，秋季の植物枯死・乾燥の経緯を反映したものと考える。
同じ期間の熱収支要素の変動をFig. 17 に表わした。
この図は Fig. 15 を全期間に拡張したもので，凡例は同図と同じである。熱収支要素の変動は毎年似たパターンで繰り返され，これが水温分布の変動と密接に応答していることが分かる。前述の全面解氷と，全熱収支の最大時期が年によらず一致していることなど，特に湖面凍結と解氷のプロセスは，年毎の量的相違はあるものの，順序とパターンは変わらない。

vol. 3: 21 - 40 (2021)

5. 討論
・インターバルカメラの画像データ
現地では瞳沼での氷の厚さは測定していないが，インターバルカメラの画像（1 時間間隔）がヒントになる。5 シーズンの中で，インターバルカメラの画像によって実際の結氷日・解氷日を特定できたのは (201314) シーズンのみである。この年の結氷日は実際には12 月 17 日であったが，計算氷厚は既に 0.105 m に達していた。
アメダス豊富地点の記録では既に 0.47 m の累加降雪と 0.08 m の積雪があり，絶対圧水位計の記録では，既に 0.182 m 沈下している。
一方，インターバルカメラの画像による実解氷日は4 月 17 日だが，計算では 2 日後の 19 日に解氷している。17 日には計算上まだ 0.077 m の氷厚があり，絶対圧水位計の記録では 0.486 m の沈下量がまだ残っている。
インターバルカメラの画像によって，浮島の沈下現
象は少なくとも，初期と終盤には湖面の凍結とは独立に進行していることが明らかである。

・熱収支要素の変動範囲
Table 3 に，5 年余りの期間にわたる瞳沼の熱収支計
算（Fig. 17）で算出された熱収支各要素の日量（日射・日射吸収・上向長波・下向長波・潜熱流束・顕熱流束）の最大値と最小値を掲げた。単位は (J cm-2 day-1 ) である。長波放射量は上向きと下向きともに，絶対値は日射量を超える日が多く，これは Fig. 17 でも明らかである。一般的には上向き・下向き長波放射はその和を正味長波放射と称して扱われることが多いため，目立たなくなってしまうが，このように分けて扱うとその絶対値の大きさが認識される。
同表には Edinger ほか [44, 45] による米国中緯度地
域の 10 ヶ所あまり (北緯 29 度～40 度) で計測された各要素の変動範囲を示した。彼らの研究は，発電所の冷却水の排水先の池沼を対象としている。本研究の対象とは境界条件が異なっている部分もあるが，多くの項目で近い値が示されている。その観測地域は瞳沼よりもかなり低緯度である上，温排水が溜まった池沼という人工的な環境である。そのため，例えば潜熱流束は全ての期間でマイナスとなっており，わずかではあるがプラスの期間 (凝結環境) を含む瞳沼での自然環
境と異なる。
太陽放射量と下向き長波放射量の相違は緯度の違い
が，その他の要素は温排水池であることが反映された相違であるものと考えられる。

・浮島沈下過程と計算氷厚
Fig. 18 は，浮島の実測沈下量 (横軸) と計算氷厚(縦軸) の関係を図示したものである。5 回の冬期間ごとの変化と，全 5 シーズンを重ねた図を示した。図中の赤い線の部分は各期の増大期，即ち計測された沈下量，計算された氷の厚さともに増加している期間である。黒線部が減少期間に相当する。Fig. 18 各図の左上に示した式は，モデルによる氷厚の計算値 𝑦 への浮島の実測沈下量 𝑥 からの回帰式である。その下に書かれた 𝑟 は，この直線的な関係性の強さを表わす
相関係数の値である。6 枚の図にはそれぞれの回帰式
を青い直線で示した。
これらの線の形を詳細に見ると，黒く描かれた減少
期の線が滑らかなのに対して，増大期の赤い線にはギザギザが多い。それも横軸方向へ短時間でずれる動きを繰り返している。これは実測の浮島沈下量が短時間で増加していることを示し，現地現象としては降雪に

表 3. 熱収支要素の分布範囲
Table 3. Value ranges for the heat balance factors calculated over more than 5 years.

35

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

図 18. 観測された浮島の沈下量と計算された氷厚との関係．赤実線：増大期，黒実線：減少期，緑実線：5 年間履歴曲線の減少期．
Fig. 18. Relations between observed islet's sinking depth and calculated ice thickness in May 2012-October 2017.

当たるものと考えられる (Fig. 10)。
これに対比される計算氷厚（縦軸）に関しては，使っているモデルは熱収支に基づいた，凍結の消長を対象としている。降雪というイベントを熱的な収支として考慮しているものの，その重量を考慮していないた
め，計算による変化が緩慢になるものと考えられる。
増大期を表わす赤い線は極大期に近づくと，急勾配となり，この傾きが 45 度に近くなる年が多い。これは計算氷厚が，浮島の沈下量に近いペースで，増加していることを示している。この時期になると浮島は十分
沈み切り，浮島表面も氷に覆われた状態となって，湖面の氷と一体となった挙動をするようになっている。

浮島沈下は計算氷厚の増加に近いペースで進んでおり，現実に近い計算結果が得られているものと考える。この事実から当モデルで求められる氷厚計算値の信頼性は高いものと考える。
Fig. 18 の右下の図は 5 年間の履歴を緑の実線で表わした図である。赤で描かれた増大期では回帰直線に沿って変動している年が多いが，緑で描かれている減少
期では回帰直線から大きくずれている年が含まれる。
これらはいずれも融雪期に浮島が浮上する過程で生じていることから，浮島と氷との相互作用，例えば潜り込みとか，乗り上げなどの現象が生じた可能性が考えられる。

表 4. 実測された浮島の最大沈下量と計算された最大氷厚．
Table 4. Observed maximum sinking depths, calculated maximum ice thickness and their appearance dates.
.

36

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

Table 4 に浮島の最大計測沈下量と，計算された最大氷厚を年毎に示した。いずれの年も計測沈下量の最大値が計算氷厚の最大値を上回っている。括弧書きしてある生起日を比較すると，大半の年で計算氷厚が最大となる日の方が早い。逆となる年はその差がわずかである。
これらの結果から，氷厚の増加以外に浮島の沈下に寄与する要素が必ず働いており，年による多寡はあるが，それは降雪の荷重であろう。

vol. 3: 21 - 40 (2021)

6. 結論
本研究で対象とした瞳沼には，浮島と呼ばれる泥炭の塊りが浮いて漂っている。この湖沼で計測した水温分布の経年的な変動と，浮島が冬期間に沈下・沈没する現象を連続計測した。一方，このメカニズムを説明するシミュレーションモデルを構築し，複数年にわたって同じ境界条件を与え続けて計算した。
結果として，水温鉛直分布の変動では，計算値は実測値に極めて近く，高い再現性を示した。この中で得られた熱収支要素の変動によれば，解氷の最大要因は水面が氷で覆われていることであることが，定量的に明らかになった。
これは湖面が凍結することによって，冬期間の湖沼からの放熱が防がれていることを示している。ここから類推すれば，湖面凍結の有無は，春季以降の湖沼水温の変動に少なからぬ影響を及ぼすであろう。今後，表面凍結が起こりにくい湖沼の熱収支変動へモデルを適用して，適用限界を探っていきたい。
もう一つの成果は計算氷厚である。実測として計測されている浮島沈下量は，氷の厚さを知る間接的な根拠となる。この両者の変動を比較した結果，浮島の沈下現象は，少なくとも初期と末期には氷の厚さとは無関係に変動することが明らかになった。その一方で，氷厚の極大期には両者は一体となった変動をしている。
これらから，用いたモデルは湖沼の氷の厚さを予測する手段としても合理性を示している。
また浮島の沈下量，ひいては氷厚についても，降雪の荷重に敏感に反応していることが判った。使用した
現モデルは降雪の熱的側面にのみ配慮しているが，荷重については考慮していない。今後のこの降雪荷重の氷厚増加への寄与を，モデルに組み込むべきか否かの検討が課題となる。その確証を得るためにも，氷の厚さのみの計測データの蓄積を図っていきたい。
本稿で用いた熱収支と結氷にかかわる計算モデル，「鉛直一次元静的多層モデル」は水温分布の変動・氷厚の変化について高度な再現性を示した。今後，このモデルが深さや広さの異なる，どのような環境条件の湖沼への適用が可能であるかを確かめていきたい。さらに，このモデルとは別に改良中の「水平二次元動的多層モデル」との機能の差などとも比較しつつ，互いの適用限界を見極めていきたいと考えている。
本稿で使用したプログラムのフローチャートを巻末に掲げた (Appendix 1.)。
なおこの手法を他の湖沼に適用されたい方には，計算プログラム (Fortran) と作図プログラム (VB-net)を提供する用意がある。(okadax@sage.ocn.ne.jp)

・浮島沈下のシナリオ
これまでに得られた知見に基づいて，浮島の沈下・沈没のシナリオを考えていく。先に述べたように，浮島はその上に降雪があると，それによる荷重増加分だ
け沈下する。まだ湖面が十分に凍結する以前に，直ぐに融けてしまわない降雪があれば，湖面の凍結とは独立に，浮島は単独で沈下し始める。
荷重増加で浮島表面が湖面近くにまで沈むと，浮島
表面の低いところに侵入した沼水も凍って，積雪荷重
を氷の荷重に置き換えるプロセスが進んでいく。この時期には日毎に寒さが増していくので，氷が解けて後戻りすることは少ない。これが沈む一方で推移する浮島沈下の理由である。
さらに時間が経過すると，浮島も湖面の氷も一体となって降雪荷重を受け，氷の亀裂から浸みだした沼水が凍る。こうして，氷厚も沈下量も増大を繰り返していく。Fig. 10 における，浮島沈下量（黒実線）と累加
降雪量（赤実線）のグラフが，並行してよく似た推移を示している部分は，この段階に相当するものと考えられる。
このシナリオに沿って進行してきた累加降雪量と，極めてよく似た推移を示していた浮島の沈下量が，1月末から 2 月ころ急速に離れていく。この現象が起こる時期には，全熱収支がマイナス（出）から急激にプラス（入）に転じる (Fig. 15 の緑実線)。このような状態になると，降雪があっても融けるか，昇華蒸発してしまう。氷が厚さを増すことはなくなり，固体としての荷重も増えず，浮島沈下量は横這い状態となる (Fig.
10 の黒実線)。
さらに 4 月に入る頃には，全熱収支量は年間を通じてプラスの極値に達し (Fig. 15)，氷厚は急激に減じ，沈下量は急減，すなわち浮島は急浮上する。しかしこの時点でも，浮島の上には融け残った雪や氷の荷重が存在し続けるため，湖面の氷が解けきっても，沈下量は直ぐにはゼロに戻らない。

37

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

6. Ryan PJ, Harleman DRF (1971): Prediction of theAnnual Cycle of Temperature Changes in a Stratified
Lake or Reservoir. M.I.T. Rep.137.
7. Ryan PJ, Harleman DRF (1973): An Analytical andExperimental Study of Transient Cooling Pond
Behavior. M.I.T. Rep.161.
8. Hurley KA (1977): Vertical heat transport mechanisms
in lakes and reservoirs. M.I.T. Rep.227.
9. Hurley KA, Watanabe M, Adams EE, Jirka GH,
Helfrich KR, Harleman DRF (1979): Mathematical
Predictive Models for Cooling Ponds and Lakes. M.I.T.
Rep.262.
10. 安芸周一・白砂孝夫 (1974): 貯水池水質の調査と解析（その 2）―成層型貯水池の流動形態と選択
取水の数値解析．電力中央研究所報告 74506.
11. 宮永洋一・秋元保・石橋毅・白砂孝夫・斉藤茂・
尾崎幸男 (1979): 貯水池流動形態のシミュレーシ
ョン解析手法．電力中央研究所報告 378022.
12. 芦田和男・江頭進治 (1976): 密度躍層のある場における濁水の挙動に関する研究（2）．京都大学防災研究所年報 19-B: 311-327.
13. 岩佐義朗・松尾直規・井上素行 (1978): 貯水池における濁度解析について．京都大学防災研究所年
報 21-B: 319-329.
14. 小島崇・中村恭志・石川忠晴 (2008): 水温成層を考慮した貯水池内流動解析に向けた CIP-Soroban
法に基づく鉛直 2 次元数値流動モデルの開発．水工学論文集 52: 1039-1044.

謝辞
本研究の対象地域は，利尻礼文サロベツ国立公園の中でも，最も重要度の高い特別保護地区に指定されている。このため，環境省稚内自然保護官事務所の理解と全面的な協力を得て，現地調査を実施することができた。
浜田誠氏からは瞳沼浮島の存在情報の，佐々木伸
宏・村元正己･稲垣紘順・稲垣淳子の各氏からは瞳沼の位置やルートに関する情報の提供を受けた。
瞳沼における現地観測，データ回収には，高橋英紀・
千田智基・辻田香織・賀勢朗子・石橋亮介・片岡澄江・
佐藤大介・伊藤和人・永野雄大・山上佳祐・清原修・
中島信人・嶋崎暁啓・渡邊光晴・原澤翔太・横地譲・
青山留美子・宮永敏明・西村収・川崎正大・茂木京菜の各氏の協力を頂いた。調査に際しては，環境省サロ
ベツ湿原センター及び，NPO 法人サロベツ・エコ・
ネットワークの皆さまに便宜を図って頂いた。
また，北海道大学大学院農学研究院の井上京教授と環境省稚内自然保護官事務所からは計測機材と気圧補正データの提供を受けた。ここに謹んで感謝を申し上
げます。
参考文献
1. 土屋和三・藤田昇・田端英雄・松井淳 (1981): 深
泥池浮島の生態学的研究．深泥池の自然と人．京都市文化観光局文化財保護課，京都: 95-133.
2. 深泥池七人委員会編集部会 (2008): 深泥池の自然と暮らし―生態系管理をめざして―．サンライズ
出版，彦根，247pp.
3. Sasser CE, Gosselink JG, Swenson EM, Swarzenski
CM & Leibowitz NC (1996): Vegetation, substrate andhydrology in floating marshes in the Mississippi riverdelta plain wetlands, USA. Vegetatio 122: 129-142.
4. Zaccone C, Lobianco D, Shotyk W, Ciavatta C,
Appleby PG, Brugiapaglia E, Casella L, Miano TM &D’Orazio V (2017): Highly anomalous accumulation
rates of C and N recorded by a relic, free-floating
peatland in Central Italy. Scientific Report 7: 43040.
https://doi.org/10.1038/srep43040
5. Concheri G, Stevanato P, Zaccone C, Shotyk W,
D’Orazio V, Miano T, Piffanelli P, Rizzi V, Ferrandi C &Squartini A (2017): Rapid peat accumulation favours
the occurrence of both fen and bog microbialcommunities within a Mediterranean, free-floating peatisland. Scientific Report 7: 8511.
https://doi.org/10.1038/s41598-017-08662-y

https://doi.org/10.2208/prohe.52.1039

15. Zhang L, Zhang J, Peng Y, Pan J, Peng Z (2019):
Numerical Simulation of Flow and Temperature Fields
in a Deep Stratified Reservoir Using Water-Separating
Curtain. International Journal of Environmental
Research and Public Health, 16, 5143: 1-16.
https://doi.org/10.3390/ijerph16245143
16. 加藤敬 (1982): 農業用貯水池における水温成層の予測．土木学会第 26 回水理講演会論文集: 543-547.
https://doi.org/10.2208/prohe1975.26.543

17. 高津治・中村由行・早川典生 (1991): 湖沼の水温構造の予測モデル．水工学論文集 35: 179-184.
https://doi.org/10.2208/prohe.35.179

18. 谷慧亮・梅田信・朝岡良浩・山崎剛 (2014): アン
デス高地の貯水池における水温変化特性に関する研究．土木学会論文集 B1 70 (4): I_1627-I_1632.
https://doi.org/10.2208/jscejhe.70.I_1627

19. 佐渡公明 (1983): 河川水の熱収支に関する基礎的研究．土木学会論文報告集 330: 69-79.
https://doi.org/10.2208/jscej1969.1983.330_69

38陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

20. 辻井達一・橘ヒサ子・高橋英樹・梅沢俊・岡田操・冨士田裕子 (2003): 北海道の湿原と植物．北海道
大学図書刊行会，札幌，264pp．
21. 梅田安治・清水雅男 (2003): サロベツ泥炭地形成
図説明書．北海道土地改良設計技術協会，札幌，
12pp.
22. 岡田操 (2010): サロベツ湿原の瞳沼とその形成過程．湿地研究 1 (1): 55-66.

vol. 3: 21 - 40 (2021)

34. Starosolszky Ö (1974): Lake Hydraulics. Hydrological
Sciences Journal XIX (1): 99-114.
35. Budyko MI (1958, 内島善兵衛訳): 地表の熱収支．
河川水温調査会，東京, 181pp.
36. Brandt RE, Warren SG, Worby AP, Grenfell TC (2005):
Surface Albedo of the Antarctic sea Ice Zone. Journal of
Climate, American Meteorological Society 18:
3606-3622.
37. Swinbank WC (1963): Long-wave Radiation from
Clear Skies. Quarterly Journal of the Royal
Meteorological Society 89: 339-348.
38. Rohwer C (1931): Evaporation from Free Water
Surfaces. United States Department of Agriculture,
Technical Bulletins No.271: 1-95.
39. 岡田 操・佐藤大介 (2012): 稚咲内湖沼群における熱収支と水収支．第 15 回日本陸水学会北海道支部
大会予稿集．
40. 岡田操 (2014): 熱収支法を利用した水収支の推定．
日本湿地学会第 6 回大会 予稿集．
41. 岡田操 (2014): 湖沼群における熱収支と水収支．
サロベツ湿原と稚咲内砂丘林帯湖沼群 (冨士田裕子編)，
第Ⅱ部第 3 章，170-177，
北海道大学出版会，
札幌市．
42. 東海林明雄 (1969): 雪面における昇華蒸発と輻射．
雪氷 31 (3): 69-75. https://doi.org/10.5331/seppyo.31.69
43. 小嶋真輔・中村和樹・直木和弘・榎本浩之 (2008):
浸み上がりおよび冠水が海氷上における積雪の温度へ与える影響．雪氷 70 (4): 423-433.

https://doi.org/10.24785/wetlandresearch.WR001006

23. 岡田操 (2013): 浮島が浮かんだ話．モーリーNo.33，
36-39，北海道新聞野生生物基金，札幌．
24. 岡田操 (2017): 浮島はなぜできる．図説日本の湿地（日本湿地学会監修），第Ⅲ部-2，130-131，朝倉
書店，東京．
25. 羽田良禾 (1937): 泥炭地池沼の研究 Ⅱ.擬底による浮島の生成．陸水学雑誌 7 (2): 64-67.
26. 吉岡邦二 (1954): 尾瀬ヶ原濕原植物群落の構造と発達．
「尾瀬ヶ原総合学術調査団研究報告 尾瀬ヶ原」，日本学術振興会: 170-203.
27. 大竹亮作 (1970): 福島県「蓋沼の浮島」の絶対年
代．地球科学 24 (1): 41-42.
28. 岡田操 (2014): 浮島が沈んだ話．モーリーNo.35，
40-43，北海道新聞野生生物基金，札幌．
29. 山崎剛 (1998): 厳寒地に適用可能な積雪多層熱収支モデル．雪氷 60 (2): 131-141.
https://doi.org/10.5331/seppyo.60.131

30. 遠藤八十一・小南裕志・山野井克己・竹内由香里・
村上茂樹・庭野昭二 (2004): 降水量データから積
雪深と密度を推定する方法．雪氷 66 (1): 17-25.

https://doi.org/10.5331/seppyo.70.4_423

44. Edinger JE, Geyer JC (1965): Heat exchange in theenvironment. Cooling water discharge Research Project(RP-49), Edison Electric Institute.
45. Edinger JE, Brady DK, Geyer JC (1974): Heatexchange and Transport in the environment. Cooling
Water Discharge Research Project (RP-49), ElectricPower Research Institute.

https://doi.org/10.5331/seppyo.66.17

31. 毛利茂男 (1946): 松原湖の結氷に就て．雪氷 9 (5):54-62. https://doi.org/10.5331/seppyo.9.54
32. 大八木英夫・市丸康平・中川麻悠子 (2015): 涌池
における水温躍層の形成機構．日本大学文理学部
自然科学研究所研究紀要 50: 125-135.
33. Chen TA, Millero FJ (1986): Precise thermodynamic
properties for natural waters. Limnology andOceanography 31 (3): 657-662.

受付：2020 年 11 月 10 日
受理：2021 年 7 月 11 日

39

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

vol. 3: 21 - 40 (2021)

付録 1. 「鉛直一次元熱収支 静的多層モデル」の計算フローチャート
Appendix 1. Calculation flowchart of “Vertical one-dimensional heat balance model with static multi layers”

Start
Input data① Shape of pond
② Headline of meteorological factors
Output ⑤ Headline of water temperature of each layer
Output ⑥ Headline of heat storage
Output ⑦ Headline of heat balance

Daily loop
Hourly loop
Input data② hourly meteorological factors

Divided loop
I: Calculation of water level and flow
II： Calculation of vertical and horizontal flow
III： Calculation of heat flow
III-① Surface layer heat transfer/ freezing process
1 Surface layer heat diffusion
2・3 Surface layer heat advection
4 Surface layer solar absorption
5 Surface layer heat exchange by evaporation/ conduction/ long wave radiation
6 Surface layer heat transfer from precipitation
7 Surface water temperature change
8 Surface layer freezing algorithm
III-② Middle layers heat transfer
1 Middle layers heat diffusion
2・3 Middle layers heat advection
4 Middle layers solar absorption
III-③ Bottom layer heat transfer
1 Bottom layer heat diffusion
2・3 Bottom layer heat advection
4 Bottom layer solar absorption
5 Bottom layer heat exchange with geothermal
IV： Calculation of vertical water temperature distribution
IV-① Adjustment of water temperature distribution around 4 ℃
IV-② Elimination of negative water temperature
IV-③ Convective mixing/ Density stratification
Calculation of cumulative values

Calculation of deep convection (deep lake)
Daily results additional output⑤ Water temperature of each layer⑥ Heat storage ⑦ Heat balance factors

Final results output⑥ Total heat storage ⑦ Total heat balance

End

40

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

2021 年度会務報告
運営委員長・知北和久
I. 全国大会の開催
2021 年度陸水物理学会第 42 回大会は，2021 年 11 月
20 日（土）～21 日（日）に名古屋市の南山大学で開催された（世話人・大八木英夫氏）。今回は日本陸水学会
東海支部会が共催として参加し，同会員の研究発表があった。新型コロナウイルスの感染拡大が収まらない
ため，感染予防対策を十分に行い，初めての試みとして対面+オンラインのハイブリッド形式で実施した。
参
加人数は計 30 名，このうちオンライン参加者 6 名，対面参加者 24 名であった。また，研究発表の件数は 18件であった。なお，巡検の開催は中止された。

（2）執筆要領などの変更
DOI を記載することの審議と了承があり，執筆要領及びテンプレートを修正し、HP へアップロードした。
（3）論文原稿の受理（受付）に関して
執筆要領及びテンプレートに示している書式などが大幅に異なっている場合には，受理前に編集委員会（編集委員長）の判断で、著者に修正依頼
し、修正原稿を送ってもらうこととした。
（4）検討課題
図の引用など，著作権に関わるイレギュラーな対応が生じた場合，他の学会誌の投稿規定や執筆要領などを調べておきたい。
3）会計報告2020 年度および 2021 年 11 月までの収支決算について，事務局・知北より報告があり，これが了承された。
4）新会員
2020 年度に岡田操氏，2021 年度に苗村晶彦氏，鈴木啓明氏が入会し，2021 年 11 月現在の会員数は 137名となった。
5）会長印，大会実行委員会印の作成
6）地学辞典 改訂版（平凡社）項目提案の学会として
陸水物理学会の参加
2023 年 10 月発行予定
7）MDPI との契約について
スイスのオープンジャーナル学術出版社
MDPI(Multidisciplinary Digital Publisher Institute)と学会パートナーシップを締結した(2021 年 10 月付)。
Water,
Hydrology に受理された論文は、掲載料が 15%OFF。
投稿の際は，当会の学会パートナーシップを申し出る。

II. 運営委員会報告
全国大会の開催に先立つ 2021 年 11 月 19 日（金）に，南山大学で（オンライン+対面）のハイブリッド形
式で運営委員会が開催された。
出席者は 9 名(オンライン 3 名，対面 6 名)であった。
議題：
1）次期会長について
審議の結果，現会長の鈴木啓助氏の再任が了承された。任期は，2022 年 4 月 1 日～2025 年 3 月 31 日の3 年間。
2）2022 年度の開催地について候補地として，北海道・釧路市（世話人，室蘭工業大学・中津川誠氏）
日程（案）：2022 年 11 月 5 日（土）～6 日（日），
巡検先（予定）：釧路湿原国立公園ほか
3）陸水物理学会功績賞の設置と候補者
当会に功績賞を設けることが審議され，徳永英一会員，北岡豪一会員，諏訪浩会員の 3 名に授与することが了承された。今後，この賞を含め学会賞の規程を事務局で作成することとなった。
4）新たな運営委員について新たな運営委員に金沢大学・坂田義隆会員が推薦された。

III. 総会報告
2021年11月20日（土）に南山大学にて総会が開催された。

報告：
1）日本学術会議協力学術研究団体に指定
2021 年 3 月 25 日付で，日本学術会議から正式に協力学術研究団体の指定を受けた。
2）編集委員会からの報告
藤井智康・編集委員長から以下の報告があった（1）第 3 巻の掲載状況について

議題：
1) 次期会長について
2) 2022年度の開催地について
3) 陸水物理学会功績賞の設置と候補者
4) 新たな運営委員について
以上の議題 1)～4)について，事務局・知北より運営委員会での審議内容が説明され承認された。
41

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

報告：1) 日本学術会議協力学術研究団体に指定
2）編集委員会からの報告
3) 会計報告4) 新会員5) 会長印，大会実行委員会印の作成
6) 地学辞典 改訂版（平凡社）項目提案の学会として
陸水物理学会の参加
7) MDPIとの契約について以上の報告事項 1)～7)について，事務局・知北より報告があり，承認された。

42

「陸水物理学会誌」投稿規程
陸水物理学会誌（英文名：Journal of the Japanese Society of Physical Hydrology）は陸水物理学会の機関誌で年 1 回から数回発行される。本誌には，広く陸水学に関する報文，その他を掲載する。本誌の発行は，学会のホームページ上に pdf ファイルの形で公開される。
１．投稿者
第一著者あるいは責任著者は正会員に限る。ただし，編集委員会が依頼する原稿についてはその限りではない。
２．投稿の種類
陸水およびその関連分野の報文（原著，総説，短報，報告）および書評とする。
(a)原著：独創性あるいは新規性のある研究論文で，価値ある結論あるいは事実を含む未発表のものとする。
(b)総説：その分野の研究の進捗状況，現状，展望などを総括した論文で，全体として一つのまとまった主張が展開されている未発表のものとする。
(c)短報：新しい事実や価値ある内容を含む論文，あるいは速報的価値のある内容やデータを含む論文で，未発表の比較的短いものとする。
(d)報告：陸水学に関する野外調査報告，講演記録，シンポジウム報告などの学術的な報告とする。
(e)書評：書籍などの書評。対象書籍等の著者が会員，または書評執筆者が会員とする。
３．用語と制限ページ
投稿原稿で使用する言語は日本語と英語の両方とし，日本語の場合は英語のタイトルとアブストラクトをつける。図表中の言語は英語とし，本文で日本語を使用の時は図表の説明
は英語と日本語の併記とする。
原著および総説の刷上りは 20 ページ程度以内，短報および報告は 10 ページ以内，報告は 5 ページ以内とする。カラーの図・写真は掲載可能で，これによる著者へ負担は生じない。書評の刷上りは 2 ページ以内とする。
４．原稿の作成
学会ホームページから Microsoft Word ファイル形式の投稿用テンプレートと投稿票をダ
ウンロードし，別に定める執筆要領に従い本文，図表，引用文献など作成する。また，必要事項を投稿票に記入する。
５．原稿の投稿と受付
43

(a) 投稿方法
原稿および投稿票を電子メールに添付して編集委員長に送信する。
(b) 受付通知
受付の可否は，原稿が投稿用テンプレートの様式に従って作成されているかどうかの審査後，原則として 1 週間以内に，編集委員長から投稿者に通知される。期間を過ぎても受付通知が到着しない場合は，編集委員長に連絡すること。
６．原稿の査読と受理または却下の通知受付けられた原著，総説，短報，報告の原稿は，複数の専門家による査読を受け，掲載の可否が判断される。査読の過程では原稿ごとに担当編集委員が置かれ，著者との連絡を行う。
査読に要する期間は 3 週間程度以内である。査読の結果，修正を要すると判断された場合，担当編集委員はその内容を著者に伝え，修正を求める。修正原稿または査読結果への反論は，原則として 3 週間以内に再投稿しなければならない。修正原稿提出の際，著者は査読者からの修正意見毎に修正内容，修正箇所，反論などを明示した修正リストを提出すること。特別の理由なく 1 ヶ月を経過した場合は，その原稿は著者が取り下げたものと判断される。陸水物理学会誌への報文の掲載可否は，査読の最終結果を踏まえ，編集委員長が最終的に判断する。著者への受理または却下の通知は，編集委員長が行う。
７．清書原稿の作成と校正原稿が受理された場合，著者は編集委員会による英文校閲結果を含めて速やかに校正を行い，清書原稿を編集幹事に送付する。清書原稿における文章や図表の変更は，編集委員の指示によるもの以外は認められない。印刷原稿の著者校正は，原則として初校のみとし，再校以後は編集委員会が行う。
8．掲載原稿の著作権(a) 本誌に掲載された原稿の著作権は，陸水物理学会に帰属する。
(b) 陸水物理学会による著作権の行使は，電子化および電子ネットワーク上での公開を含む。
(c) 本誌に原稿を投稿する者は，著作権に関する条項を事前に承諾したものとみなされる。
9．附則
本規程は，2020 年 4 月 1 日より施行する。

44

「陸水物理学会誌」執筆要領
１．文章の作成
陸水物理学会ホームページにある投稿用テンプレート（Microsoft Word ファイル）をダウ
ンロードし、その様式に従って原稿を作成する。テンプレートにおいては，日本語摘要は 1
行 48 字で MS 明朝 10pt、英語摘要は Times New Roman の 10pt で作成のこと。本文は 2 段
落１行 21 字で，日本語は MS 明朝 10.5pt，英語は Times New Roman 10.5pt で作成のこと。
原稿における式の使用はテンプレートでの Word ファイルで「挿入」「数式エディタ」から行うこと。
２．原稿の構成(a) 投稿の種別
最初のページの１行目左に投稿の種別（原著論文・短報・報告・総説）を明記する。
(b) 原著論文・短報・報告・総説の種別に続けて，次の順序で記述する。①和文表題，②和
文著者名，③英文表題，④英文著者名，⑤和文摘要，⑥英文摘要，⑦和文の所属と住所，⑧英文の所属と住所
ただし，⑦和文の所属と住所，および⑧英文の所属と住所はテンプレート 1 ページ目下の脚注に記載のこと(c) 摘要とキーワード
和文摘要は 300 字以内，英文摘要は 200 語以内で，いずれも 10 pt を使用のこと。
キーワードは和文・英文とも 5 語以内で，いずれも 10 pt を使用のこと。
(d) 本文は「はじめに」と記して書き始める。原著論文と短報は，「方法（材料と方法，観測，
調査など）
」
，
「結果」
「考察」
（または，
「結果と考察」「結果と討論」
）
「まとめ」
（または
「結論」
）
「謝辞」
「参考文献」の順で記載し，各項目に番号付けを行う。各項目の中には必要に応じて小項目を置くことができる。小項目にも番号を付ける。総説の構成は「はじめに」以外の邪分は特に指定しな いが，内容が理解されやすいように適宜項目や小項目を設定する。
３．活字指定と表示法
テンプレートでは Microsoft Word を使用しており，本文の日本語は MS 明朝体の 10.5
pt，図表は全て英語表現とし，図表中の英数字は Times New Roman 体とする。図は鮮明な画像ファイル(拡張子 jpg, gif, png などのファイル)を使用し，表は字体のサイズが変更できるオリジナルを挿入のこと。句読点は（。
）および（，）とし，和文中は全角，英文中
は半角とする。
45

図表の番号は，“図 1.”や”Fig. 1.”のように MS 明朝体と Times New Roman のボールド表
示で 10pt とし，この後に続く図表の説明文は，普通字体で和文・英文共に 10pt で書くこと。
図表を本文中で引用する場合，Fig. 1, Figs. 2-4, Table 1 などと Times New Roman のボー
ルド表示 10.5 pt で記すこと。なお，Fig., Figs., Table と番号の間には半角空白を挿入のこと。写真は図として扱う。
４．単位や日付などの表記
本文，図表ともに SI 単位を用いる。m s–1, J kg-1 K-1 のような半角空白を挿入した指数表示とし、m/s, J/kg/K のような表現は使わない。数字や欧文表記は半角文字を使用する。
日付の表示は西暦を用い，図表中やその英語説明文では mm/dd/yyyy や mm/dd の形式とする。
５．引用文献の記載(a) 本文中の引用
本文中の文献の引用は，引用順に[1], [1-3]にように記して，最後の「参考文献」にこの順番で記載する。本文中で引用する場合は“奥田[4]によれば….”のような表記とする。著者が多い文献は，“中尾ほか[5]によれば…..”のような表記とし，「参考文献」には著者全員の名を記載すること。
(b) 文献リストの作成｢文献｣の項目には本文中に引用されたものすべてを記載する。記載順序は引用順とする。雑誌名は省略しないで表記し，英語の雑誌名はイタリック体にする。なお，DOI(Digital Object Identifier)が付与されている論文については，下記の「記入方法」に従い必ず記載すること。
記入方法
ア．論文
和文の場合
「著者名（年号）：題目．雑誌名，巻：ページ．DOI」の順に記載する。著者と著者の間は「・」でつなぐ。雑誌名と巻の間は半角空白にし，巻数は半角のボールド表示とする。数字と( )は，Times New Roman を使用のこと。
例）福富孝治・中尾欣四郎・三好日出夫・田上龍一 (1968): 登別温泉大湯沼の水収支および熱収支．北海道大学地球物理学研究報告 19: 1-19. https://doi.org/10.14943/gbhu.19.1
藤井 智康・奥田 節夫 (1995): 中海・宍道湖における連系振動．陸水学雑誌 56： 291296．https://doi.org/10.3739/rikusui.56.291

46

英文の場合「著者名(年号): タイトル. 雑誌名（イタリック体）, 巻（ボールド表示）: ページ. DOI」の順に記載する。著者が 2 名以上の場合は「, 」でつなぐ。first name と middle name はイ
ニシャルのみで「．
」は省略する｡
例) McMahon A, Moore RD (2017): Influence of turbidity and aeration on the albedo of mountainstreams. Hydrological Processes 31: 4477–4491. https://doi.org/10.1002/hyp.11370
Striegl, RG, Dornblaser, MM, Aiken, GR, Wickland, KP, Raymond, PA (2007): Carbon export
and cycling by the Yukon, Tanana, and Porcupine Rivers, Alaska, 2001–2005. Water Resources
Research 43: W02411. https://doi.org/10.1029/2006WR005201
Fujii, T. (1997): Relationship of internal waves with tidal and wind drift effects and propagationof internal Kelvin waves in a brackish lake - Lake Nakaumi. Japanese Journal of Limnology 58:
241-260. https://doi.org/10.3739/rikusui.58.241
イ．単行本の場合
和文の場合「著者名（年号）：書名．出版社，発行地，総ページ」の順に記載する。
例) 吉村信吉（1937）：湖沼学. 三省堂，東京，520p．
英文の場合「著者名（年号）：書名．出版社，発行地，国名，総ページ」の順に記載する。
例) Dingman SL (2002): Physical Hydrology, 2nd ed., Prentice Hall, Upper Saddle River, USA,
646p.
ウ．単行本の章または分冊
和文の場合「著者名（年号）：表題．書名（編著者），引用ページ．出版社，発行地．
」の順に記載する。
例) 奥田節夫 (1985): 土砂災害の危険範囲予測．扇状地の土砂災害（芦田和男編），第 5
章，155-191，古今書院，東京．
英文の場合例 ) Mosley MP, McKercher AI (1993): Stremflow. In Handbook of Hydrology, Chapter 8,
Maidment DR, ed., McGrwa-Hill, New York, USA, pp. 8.1-8.39.
エ．WEB ページ
和文の場合例 ) 気象庁 (2018)：日本活火山総覧 （第 4 版）． URL: https://www.data.jma.go.jp/svd/
vois/data/tokyo/STOCK/souran/menu_jma_hp.html (2019 年 1 月 8 日時点)

47

英文の場合例 ) National Weather Service (2019): Snow depth. URL: https://www.weather.gov/nerfc/
snow_depth_im (2019 年 1 月 7 日時点)

48

「陸水物理学会誌」査読規程
1. 目的
この規程は，陸水物理学会誌の質を高め，陸水物理学および関連分野の一層の発展をはかるために定めるものである。
2. 論文の受付
論文の投稿があった場合，編集委員長は投稿論文の内容や原稿形式等を確認し，受付の可否を判断する。受付可と判断された場合には，担当編集委員を定めて論文の査読を依頼するとともに，投稿者に論文の受付日を連絡する。担当編集委員は査読者 2～3 名に論文の査読を依頼する。
3. 査読内容
査読者は，以下の基準 1）～6）に基づいて査読し，客観的に掲載適否の判定を行うこととする。査読結果は別紙 1 の査読報告書にまとめて担当編集委員に送付する。査読に要する期間は概ね 3 週間とする。
1） 査読者は専門的見地から「論点及び事実の新規性」，「論拠の明示性」，「方法の再現性」，「考察の妥当性」，「全体を通した論理性」について意見及び判定案をまとめる。
2） 査読は内容の批判や討議を行うものではない。したがって見解の相違は掲載不適の理由とはならない。
3） 査読者は，著者に対して研究指導する立場にない。よって，投稿された原稿に対して，新たな調査や実験等の追加要求は極力避けなければならない。
4） 投稿された成果ができるだけ早く公表されるよう，査読者は査読期間を遵守しなければならない。
5） 査読者は査読により知り得た掲載前の原稿の内容を，いかなる形でも利用・剽窃・漏洩してはならない。
6） 個々の投稿原稿についての査読者名は原則として公表しない。また，著者との折衝はすべて編集委員会が行うこととし，査読者が著者と直接に折衝してはならない。
4. 査読結果の報告
担当編集委員は，査読者の査読報告書を参考として投稿論文の審査を行う。この審査で修正等が必要と判断された場合には，担当編集委員は投稿者に修正内容を連絡するとともに査読報告書を送付し修正原稿の提出を依頼する。修正原稿の提出のないまま特別の理由な
く 1 ヶ月を経過した場合は，その原稿は著者が取り下げたものと判断される。

49

5. 再審査
担当編集委員は，投稿者から提出された修正原稿の再審査を行い修正が十分に行われて受理可能と判断された場合には，編集委員長にその審査結果を報告する。再審査の結果，査読者による再査読が必要と判断された場合には，査読者による再査読を行い，その結果を投稿者に連絡する。
6. 論文の受理
編集委員長は，担当編集委員から提出された審査報告書を基に担当編集委員と協議して論文の受理を決定する。論文が受理された場合には，編集委員長は投稿者に受理通知を送付する。論文が受理できない場合には，その理由を投稿者に連絡する。
7. 最終原稿の提出
投稿者は論文が受理された後， 編集委員長に受理論文の最終原稿を電子ファイルで提出する 。なお，編集委員長は最終原稿について，原稿の訂正を投稿者に求めることがある。
8. 印刷原稿の校正
オンライン版原稿の訂正を行う場合，印刷上の誤り以外の字句の訂正，挿入，削除は原則として認めない。
9. 査読規程の制定及び変更
この規程の制定及び変更は運営委員会の決議を経て行うものとする。
2019 年 3 月 1 日制定

50別紙１（査読報告書様式）
年
査読報告書（査読者 A 第１審）
受付番号：XXXX（原稿番号）
投稿分類：■ 原著
□ 総説
□ 短報
□ 報告
□ 書評
表 題 ：
著 者 ：
査読者 ：
判定結果：□ このまま掲載してもよい□ 指摘箇所を編集委員会で修正し掲載してよい□ 著者に戻して修正させ，編集委員会で確認し掲載してよい□ 著者に戻して修正させ，再度査読する必要がある□ 掲載は不適当である
編集委員会への意見（著者に対しては開示されません）

著者へのコメント：
1）全体的意見

2）個別的意見（ページ・行，指摘事項）

51

月

日

陸水物理学会誌 (Journal of the Japanese Society of Physical Hydrology)

編集後記

陸水物理学会誌第 3 巻をお届けします（表紙は青森県・青池の写真。陸水物理学会誌第 1 巻及び第 2 巻に掲載された論文著者である花石竜治会員より提供頂いた）。会員の皆様の熱心な論文ご投稿により，今巻も 2編の論文を掲載することが出来ました。
2021 年度から前編集委員長から委員長を引継ぎ，まだまだ分からないことばかりですが，皆さまのご協力を頂きながら本誌の編集，発展に力を注げればと思います。また，本学会は 2021 年 3 月に日本学術会議協力
学術研究団体に指定されるなど，本学会及び本誌のますますのアピールが必要になるかと思います。会員の皆様のさらなるご尽力をどうぞよろしくお願い致します。

コロナ禍の影響で，我々の分野で重要なフィールド
ワークに出かけることへの制限や学会もオンラインや
ハイブリッド（対面＋オンライン）などまだまだ不都
合なことばかりであり，新型コロナウィルス感染症の終息を願うばかりです。
この様な状況ではありますが，本学会の発展と世間にインパクト与え，また個性ある学会誌をつくって行
くために会員の皆様のご協力をどうぞよろしくお願い致します。

（編集委員長 藤井智康）

陸水物理学会
事 務 局

会 長
運営委員長

編集委員会

委員長
委 員
委 員

陸水物理学会誌 第 3 巻 第 1 号
発行

鈴木 啓助（信州大学）
知北 和久（北海道大学）
藤井 智康（奈良教育大学）
柴田 智郎（京都大学）
大八木 英夫（南山大学）

2021 年 12 月 1 日発行

陸水物理学会〒001-0021 札幌市北区北 21 条西 11 丁目
北海道大学北極域研究センター知北和久気付
電話・FAX
TEL：011-706-9074(代) FAX：011-706-9623
E-mailchikita@sci.hokudai.ac.jp
URL
http://jsph.sub.jp※本冊子内容の無断の複製・転載を禁ず

52

ISSN 2435-3043

Journal of the Japanese Society of Physical Hydrology
Vol. 3, No. 1

December, 2021

Contents
...

1

...

3

...

21

Announcements........................................................................................

41

Instructions for Authors........................................................................................

43

Submission Guidelines........................................................................................

45

Peer Review Regulations........................................................................................

49

Editor’s Postscript. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..

52

Foreword. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . Hiroshi SUWA

Original Paper
Estimating transpiration from a catchment by sap flow measurements: Comparison with 3D numerical simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . Takeshi SUZUKI and Kazuhisa A. CHIKITA
Original Paper
Intra- and inter-annual thermal variations of a pond in a subarctic mire
- Relations to the vertical movement of a floating islet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Misao OKADA

Published by

The Japanese Society of Physical Hydrology
Hokkaido University Arctic Research Center
Kita-21 Nishi-11 Kita-ku, Sapporo 001-0021, Japan

http://jsph.sub.jp/