日本ヒートアイランド学会論文集 Vol.5 (2010)
Journal of Heat Island Institute International Vol.5 (2010)

学術論文

都市散水効果の定量評価に向けた
大気・地下連成‐水・熱環境モデルの開発
Development of an Atmosphere-Subsurface-Coupled, Hydrothermal Environment Modelfor Quantitative Evaluation of the Effect of Urban Water Sprinkling
吉岡 真弓*1
Mayumi YOSHIOKA

登坂 博行*2
Hiroyuki TOSAKA*1 産業技術総合研究所 Geological Survey of Japan, AIST
* 東京大学大学院工学系研究科 Graduate School of Engineering, The University of Tokyo
2

Corresponding author: Mayumi YOSHIOKA, yoshioka-mayumi@aist.go.jp
ABSTRACTFor quantitative evaluation of the effect of water sprinkling as a mitigating method of urban “Heat Island”, acomprehensive numerical model is proposed. The developed model describes both air-water-vapor movement andheat transport in atmosphere and subsurface, by using the generalized concept of multi-phase, multi-component
solid-fluid system. It includes hydrodynamic phenomena, such as pseudo static air flow in atmospheric boundary
layer, the heat balance of solar radiation and long-wave radiation, evaporation, transient heat transfer among thesolid, liquid and gas phases. In order to confirm the performance of the model, several case studies have been maderegarding the basic modeled processes. The calculated results are reasonable and consistent with the analyticalsolutions and/or laboratory experiment with the reasonable parameters
キーワード: 数値モデリング，散水，気温低下，定量評価
Key Words : Numerical model, Water sprinkling, Temperature decreasing, Quantitative evaluationて屋外散水実験において多点温度観測を行い，0.8 ºC 程

１．はじめに度の気温低下を報告した(1)．
都市域での舗装面の増加や人間活動・エネルギー消費

散水の特徴として，現存の舗装に行うことができ即効の集中化を主因とする夏季の高温化現象（ヒートアイラ

性があること，保水性舗装や地下水利用との組み合わせ

ンド現象）は現代的な環境問題の一つである．その緩和による環境改善効果も見込まれること，などが挙げられ対策のため，これまで様々な方法が研究・提案され，一よう．一方，例えば大都市中心部の夏の気象条件（日射，部実施されてきた．例えば，植物を利用した対策（屋上

気温，風速，湿度など）を考えた時，どの程度の範囲に緑化や緑地化），風を利用した対策（風の道の利用，川風

どの程度の散水を行えば，どの場所でどの程度涼感効果の利用），水を利用した対策（路面散水，保水性舗装や高が生じるのか生じないのか，或いは湿度の増大などの逆反射性舗装敷設）などが挙げられる．

効果（不快感）はないのか，などといった疑問にはまだ著者らは，特に路面散水（以下，散水）によるヒート

答えられない段階であり，実用化には至っていない．今アイランド緩和効果に注目し，屋外実験や解析に関する後の議論では，太陽・大気・地圏・人工物の相互作用の(1)研究を行ってきた ．ヒートアイランド緩和対策として中にある都市環境中の水・熱過程を詳しく捉える数理技の散水に関しては，これまで多くの屋外実験が行われて術の開発が必要であり，それを通して散水の費用対効果(2)，(3)いる(4)．例えば，土屋ら は 2004 年の「大江戸打ち水を客観的に評価してゆく必要があろう．

大作戦」において散水場所付近の多点気温・湿度観測を

以上を踏まえ，本研究では，散水時に生じる都市熱環

行い，平均 0.48ºC の気温低下を報告している．また，小境変化，即ち，地表面近傍における大気・地下環境の水・

作ら(5）は東京丸の内や神田の保水性舗装への散水を行い，水蒸気・熱の動きを忠実に反映する数理モデルの開発を散水日と非散水日の気温・湿度変化を統計的に比較し，行う．本論文では，モデルの基礎方程式，基礎パラメー

夕方から夜間の散水において平均で 0.2～0.4 ºC の気温低タを詳細に論じると共に，モデルの基本的動作に関する下を観測している．著者らは 2004 年から 2006 年にかけ

検証結果について報告を行う．なお，続報にて屋外実験- 11 -

(1)の再現性の検討結果を報告する．

果も物理的に表現可能とする．

２．水・熱輸送モデルの既存研究と本研究の特徴

３. 連成モデルの基本仮定

大気の流れや熱，水蒸気輸送を扱う数値モデルでは，

3.1 解析対象とする空間および時間スケール(6)

地球規模を対象としたモデル（例えば MIROC ）や，メ

解析対象とする空間は，例えば都市空間内の一つの街

ソスケールモデル（例えば MM5(7)など），また都市気候

路や広場のスケールであり，地表からの高さ 10m オーダ

や街区といった比較的小さなスケールを対象としたモデーの接地境界層下部，地下 10m オーダー，水平方向にル (8)， (9) などがある．これらのモデルでは大気の運動や100m 程度のオーダーを考えることとする．現象を評価

熱・水蒸気輸送の追跡を主とし，地表面過程（潜熱・顕する時間スケールは特に制限はないが，例えば 1 時間平熱輸送）は境界条件・ソース項などにより表現されるこ均の観測記録の入力を基に，秒～分単位の時間分割で日とが多く，地表面下の細かな熱や流体流動は扱われない．

～週オーダーの解析を行うことを考える．

一方，地下の流体（水・空気など）や熱輸送を対象とした数値モデルとしては，MODFLOW(10)，TOUGH2(11)，3.2 モデル化の仮定

FEFLOW(12)などがあり，地下水温の追跡，地中熱利用，

上記の時間・空間スケールにおける大気・地表面・地深部地熱資源問題などの評価に用いられている．これら

下における水（液体水および水蒸気）・熱輸送を物理モデのモデルに地表面過程や大気圏の状態を反映させることル化するために以下のような基本的仮定を置く．モデル

はソース項や外部ルーチンにより可能な場合もあるが，の概念図を図 1 に示す．

一般的には考慮しないことが多い．

(1) 対象スケールにおける地表面付近の大気運動では，

近年では，大気や地表と地下を結びつけた数値モデル

コリオリ力や慣性の効果は小さく無視できるものとも開発されつつある．HYDRUS(13)は，地表面から飽和・する．また，乱流渦構造は直接扱わず，平均化した

不飽和土中の水分や溶質，熱移動を詳細にモデル化した乱流拡散場と考える．

モデルであり，可能蒸発量を与えることで蒸発過程を表(2) 接地境界層中において水平風速が卓越する場を仮定現することができる．また，萩島らの開発したし（平行流），温度・比湿変化による空気の動きはス(14 ）AUSSSMでは鉛直１次元的に大気モデル，土壌サブ

トークス型の密度差置換（圧力の拡散）として取りモデル，都市気候サブモデル，建築サブモデルを連成させ，都市高温化要因の評価などを行っている．登坂ら

扱う．

(15)，

(3) 大気中の熱・水蒸気の動きは，空気の移動に伴う乱(16)は，地表水と地下水の動きと共に，大気気温と液相水

流拡散および移流によりモデル化する．

温，固体温度を分離して追跡する多相流動数値モデルを(4) 対象とする大気・地下は，２相（気相，液相）
・2 成

開発し，実流域の河川水温・地下水温の観測値の再現に

分（空気，水）の流体（空気・液体水・水蒸気）と(17)

は広域の地下水，地表水

固体（岩石やアスファルト粒子）により構成されるの熱・流動解析に加え， 地表面上に 2 層モデル（地表層系と考える．大気中では固体および液体水の体積をと植皮層）を配置することで大気の熱収支も含めた水・非常に小さな値と考えることで，大気から地下まで

熱循環系の統合モデル化を試みている．

を同様の方程式系を用いて表現する．

成功している．また，稲葉ら

散水に注目すると，それによる局地地表面の変化が風(5) 気相，液相，固相は独立した温度を持ち，相互に非上・風下側大気中に及ぼす変化や地下で起こる昼間の蓄

定常熱交換をするものとする．これにより，地表面

物理的に扱うことが望ましいが，前記の統合型モデルに

大気放射

風
熱

おいてもそれは実現されていない．

水蒸気

熱・物質輸送（乱流状態）

・大気および地下における 3 次元流れを扱うと共に，固・蒸発
舗装

現（相変化や相間の熱交換）を詳細化する．

熱考慮する．これにより，散水などの 3 次元的効果（風

水蒸気

A 2層

φ ≈1

S w≪1

浸透
（液体水）

a

水層
(w)

G1層
G2層

w

φ=x

0 < Sw < 1・・・接地境界層の流れによる水平・鉛直方向の移流分散を

a

Sw = 0赤外放射 A1層
散水

液・気 3 相系として一般化することで，現象の物理表

φ ≈1

短波長放射

以上を踏まえ，本研究では既存研究を一歩進め，次のような特徴を持つ連成モデルを開発する．

Ana層

・・・

熱・夜間の放熱過程を，多次元的水・熱輸送過程として

Gng層

s

a
w

上・風下側の違い，高度による違い）を考慮する．
図 1 散水時の地表面付近での熱・物質移動の概念図と各領域における格子番号および相割合の模式
図．ここで，φ ：間隙率[-]，Sw：水飽和度[-]，a，・地下に関しては，多相流動として一般化し，舗装面（不透水面や保水性・透水性）だけでなく，土壌面などへの降水や散水の浸透，地盤中での蓄熱・夜間の放熱効

w，s：それぞれ気相，水相，固相を表す．

- 12 -

格子にまかれた水の温度と大気の温度の違いを表現Γa = Π a +すると共に，地下への水の浸透時の固体粒子と液体

Pa

(8)

ρa

ここで， Π p ：内部エネルギー[Jkg-1]（ただし，添え字 p

水の温度差を表現する．
(6) 地下での液体水の浸透，空気の流れは多相ダルシー則に従うものとする．また，地下での熱移動は，液は相を表し，s：固相，a：気相，w：水相のいずれか），-1-1 -1
Γa ：気体のエンタルピー[Jkg ]， cP, p ：比熱[Jkg K ]で体水，空気，水蒸気の移動に伴う移流分散と固体間ある．

の熱伝導，異相間の熱交換により表現する．各相間(5)空気および水の水理ポテンシャル
Ψ a = Pa + ρ a gzの熱交換では，熱平衡に達するまでの過渡過程も表
現する．

(9)

Ψ w = Pa − PC ( S w ) + ρ w gz

(7) 地表面付近では，短波・長波放射，顕熱・潜熱輸送(10)

ここで，Ψ p ：p 相の水理ポテンシャル[Pa]，z ：高度[m]，に伴う熱収支現象をモデル化する．

S w：水飽和率，PC (Sw)：水飽和度 Sw のときの毛管圧力[Pa]，-2g ：重力加速度[ms ]である．なお，大気中では PC は 0

４．大気・地下連成系水・熱輸送解析の基礎方程式

であり，地下間隙中では飽和率の関数として表わす．
4.1 基本物性の計算式本モデルで使用する基本的な物理量は以下の式で算出
する．
(1)水密度(18)および湿潤空気密度(19)(1)ρ w = ρ w0 {1 + c f ( P − P0 )}{1 − cT (T − T0 )}⎞m p⎛1
ρa = d* ⎜⎟R T ⎝ 1 + 0.608q ⎠4.2 全領域の共通基礎方程式
本モデルでは，固・液・気 3 相，2 相流体 2 成分系を取り扱う．ただし，空気中の水蒸気についても個別に追跡するため，3 成分系と同様の方程式系となり，未知数として気相圧力 Pa ，水飽和度 Sw，比湿（水蒸気）q，固(2)

ここで，ρw：水の密度[kgm-3]，ρa：湿潤空気密度[kgm-3]，

ρ w 0 ：基準流体密度[kgm-3]，P0：基準圧力[Pa]，T0：基準(

温度[ºC]， cf：水の圧縮率[Pa-1]，cT ：水の熱膨張率[K-1]，
md：乾燥空気の分子量（=0.028964）[kgmol-1]，R*：普遍
気体定数（=8.314）[Jmol-1K-1]，p，P：圧力[hPa またはPa]，T：温度[ºC]，q：比湿[kgkg-1]である．
(18)

液体水の粘性係数は以下(2) 水・空気の粘性係数の Andrade の式を用いる．

B' ⎞
⎛
−3
(3)
⎟ × 10T −T'⎠⎝
ただし， A' = −1.5668, B' = 230.298, T ' = 146.797

μ w = exp⎜ A'+

)

(

)

(

)

∂t

固相のエネルギー保存式：(∂ ( ρs (1 − φ ) Π s ))∇ DQ, S ∇Ts + ΘQ, s =また，空気の粘性係数は，以下の式を用いる．

∂t

(14)

気相のエネルギー保存式：
JJG∂ ( ρaφ Sa Π a )−∇ ⋅ ρa ua Γ a + ∇ ρ a cP ,a DQ ,a ∇Ta + ΘQ ,a =∂t

32

⎛ T + C ⎞⎛ T ⎞(4)
μa = μ0 ⎜ 0⎟⎜ ⎟⎝ T + C ⎠⎝ T0 ⎠ここで， μ w , μa ：水，空気の粘性係数[Pa・s]， μ 0 ：基

(

) (

)

(15)
水相のエネルギー保存式：

準状態での空気粘性係数[Pa・s]，C：サザランド定数（空気の場合，C = 117）である．

JJG∂ ( ρ wφ Sw Π w )−∇ ⋅ ρ w uw Γ w + ∇ ρ w cP , w DQ , w∇Tw + ΘQ , w =
∂t

(

) (

)

(16)

ただし，

(3)比湿および飽和水蒸気圧(19)⎛∂ ∂ ∂⎞∇ = ⎜⎜ , , ⎟⎟⎝ ∂x ∂y ∂z ⎠JJG
Tu p = ( u p , v p , wp )ρ vapor
0.622e(5)=ρap − 0.378e(6)
eSAT = 6.1078 × 10 aT (b +T )-1
-3
ここで，q：比湿[kgkg ]，ρ vapor ：水蒸気密度[kgm ]，eSAT：q≡⎛ Dr , p , xx⎜
Dr , p = ⎜ Dr , p , yx
⎜ Dr , p , zx⎝
Sa + S w = 1

飽和水蒸気圧[hPa]，e：水蒸気圧[hPa]，a,b：Tetens の係数（水面上の場合，a = 7.5，b = 237.3）である．
(4)内部エネルギーおよび気体のエンタルピー

Π p = c P, pT p

相温度 Ts，水相温度 Tw，気相温度 Ta の 6 つを扱う．解析対象となる大気から地下までのこれら各状態量について以下の 6 つの基礎方程式が成り立つ．
空気の質量保存式：
JJG∂ ( ρ aφ S a (1 − q ) )(11)−∇ ⋅ ρ a (1 − q ) ua + Θ M , a =∂t
水（液相水）の質量保存式：
JJG∂ ( ρ wφ S w )(12)−∇ ⋅ ρ w uw + Θ M , w =∂t
水蒸気の質量保存式：
JJG∂ ( ρ aφ S a q )(13)−∇ ⋅ ρ a qua + ∇ ( ρ a DM , a ∇q ) + Θ M ,v =(17)

(ただし，p = a, w)(18)

Dr , p , xz ⎞⎟Dr , p , yz ⎟Dr , p , zz ⎠⎟(19)

Dr , p , xy
Dr , p , yy
Dr , p , zy(20)

ここで，添え字 r は質量 M，または熱量 Q を表す．φ ：(7)間隙率（大気中では１）[-]，S a ：気相飽和率[-]， Θ r , p ：- 13 -

p 相，r のソース項[kgm-3s-2, または Jm-3s-1]，Dr,p：p 相，5.2 風速の算出方法

r の拡散係数テンソル[m2s-1]，t：時間[s]である．
本モデルでは，大気，地表面近傍，地下の各領域や層

接地境界層において水平風速が卓越する条件（平行流）に異なる拡散係数 Dr,p およびソース項 Θ r , p を与えること下において，温度・比湿変化に伴い生じるわずかな密度で全領域の連成解析を行う．

差による空気塊の置換をストークス型の置換と考えれば，
風速は以下の式で表現される．

4.3 ソース項
式(11)～(16)中のソース項 Θ r , p は，各状態量や層によJJG
K a ,ij∇Ψ aua = −μ

り異なり，以下の式で表す．

Θ M ,a = M a ,bnd
Θ M , w = − M w→v + M w, sp + M w,bnd

(21)
(22)

Θ M ,v = M w→v + M v ,bnd

(23)

ΘQ , s = −Qs → a − Qs → w + Qrad + Qs ,bnd

(24)

ΘQ , a = Qs → a + Qw→ a + QH + Qa ,bnd

(25)

ここで，Ka,ij：大気の置換係数[m2]である．
運動量保存則を考慮した大気の流動を本モデルに組み入れることも可能であるが，計算負荷が大きくなるという欠点を伴う．本研究では対象としている都市の熱環境を時間平均や日平均や時間平均の風速を用いて解析する
ことから，式(30)に示した定式化により風速を表現する．

ΘQ , w = Qs → w − Qw→ a + Qw, rad − QE − QH + Qw, sp + Qw,bnd (26)ただし，
M w→v = M E + M int
Q p1→ p 2 = Fp1→ p 2 + f p1→ p 2(30)

5.3 大気中の境界入出項（大気境界条件の設定）
境界において気相（空気・水蒸気）および気相の熱量の注入を表すために最も風上側（x = 1）および風下側格
子列に与える項であり，以下の式で表す．
JG G
M a ,bnd = ρ a (1 − q ) A ⋅ u a ,bnd
JG G
(31)
M
= ρ q A ⋅ u a ,bnd

(27)
(28)

である．ここで， M p ,bnd ：境界格子での風速に伴う質量流入量， Q p ,bnd ：境界格子での熱量流入量， M w→v ：蒸発量（ME：地表面からの蒸発量，Mint：地下間隙中の内部蒸発量），Mw,sp：散水による水の質量注入量，Q p1↔ p 2 ：v , bnda

JG G
Qa ,bnd = ρ a c p Π a A ⋅ u a ,bnd

相p1 と相p2 間の熱交換量（ Fp1→ p 2 ：格子内部の相間熱交

ただし，，Qrad：固相換量， f p1→ p 2 ：異なる格子間での熱交換量）

JG
A = (1 / Δx,1 / Δy ,1 / Δz )G
ここで， u a ,bnd ：流入風速[ms-1]である．定常場を作るたへの短波長放射量，QH：顕熱輸送量，QE：潜熱輸送量，Qw,rad：水相での長波長放射収支，である．式(21)～(28)の各項については次章以降で詳しく説明する．

め風下側境界でも同量が流出するように条件を与える．
境界以外の内部格子については M a ,bnd = 0 ， M v ,bnd = 0 ，Qa ,bnd = 0 である．

4.4 数値計算における領域の離散化
基本方程式(11)～(16)を数値的に解くために，図 1（右5.4 大気中の熱・水蒸気輸送の定式化

図）に示すような空間離散化を行う．大気の最下層（地表面と接触している層）をA1 層（厚さ 0.1m程度）とし，大気中における気相のエネルギー・水蒸気輸送を移流

上方へA2 層，A3 層…Ana層と次第に大きくする．また，と拡散効果の和としてモデル化する．拡散は，さらに分

地下についても，地表面層をG1 層（0.01m程度）とし，子拡散と乱流拡散の和としてモデル化する．また，水蒸下方へG2 層，G3 層，…Gng層と分割する．
次章以降は，各離散格子の大きさを Δx, Δy, Δz ，格子気の乱流拡散と温度の乱流拡散の大きさは同一と考える．

の体積を Δv = ΔxΔy Δz として，支配方程式の数値解析的一般化した乱流拡散テンソルD（ここでは式(19)の添え字

な取り扱いを説明する．

p，rについては省略する）を次のように表現する．

Dxx = Dm , p +

+
+
U
U
U
ϖ T u 2 ϖ L v 2 ϖ T w2
(32)
Dyy = Dm , p +
+
+
U
U
U
ϖ u 2 ϖ v 2 ϖ w2
Dzz = Dm , p + T + T + L
UUU((ϖ − ϖ T )uw
ϖ L − ϖ T )uv
D xy = D yx =, D xz = D zx = LUU(
ϖ L − ϖ T )vw
2
2
D yz = D zy =, U = u + v + w2
U

５．大気中（A2～Ana 層）における水・熱輸送の取り扱い
5.1 大気中における固相・水相の取り扱い
接地境界層中では，液体水の流動，固相，および水相のエネルギー輸送は生じないと考えることができる．そのため，大気中の水飽和度 Sw，空気飽和度 Sa，間隙率 φ ，
毛管圧 Pc，水相および気相の相対浸透率 krw，kra については次のような基本設定をする．

Sw = ε s ,S a = 1 − ε s (ただし，ε s = 10−4 )φ = 1 − ε ϕ (ただし，ε ϕ = 10−4 )ϖ L u 2 ϖ T v 2 ϖ T w2

ここで， Dm, p ：風速がない場合の分子拡散係数[m2s-1]，(29)

ϖ L ,ϖ T ：縦分散長，横分散長[m]である．縦分散長は流Pc = 0, krw = 0, kra = 1- 14 -

する．短波長放射の水中での吸収は無視する．

れに沿う乱流拡散効果，横分散長は流れに垂直な乱流拡
散効果を表す係数である．概念図を図 2 に示す．接地境

・長波は A1 層の湛水層に吸収される（長波長放射は実界層付近で水平風速の卓越した場での鉛直方向の熱・蒸質無限小の厚さの水に吸収される(19)．水深が下限値
Zw,Lb の場合でも同様とする）．

気輸送は，水平風速の横分散効果によって表現される．

・顕熱交換は湛水層表面と大気の間で行われる．湛水

また，一般的に乱流拡散効果は乱流のスケールに依存層が下限値の場合は，A1 層の水相は G1 層と同一温度（高度に依存）する係数であり，ここでは以下のようにとなる．

高度の関数として考える．

ϖ L = γ L z, ϖ T = γ T z・蒸発は，A1 層湛水深が下限値より大きい場合には湛(33)

水面から生じる．下限値に達した場合は，G1 層水相

ここで，γ L , γ T ：乱流スケールの高度効果を示す係数[-]，から地下多孔質体を通して蒸発が行われる．

z：高度[m]である．これらの値の算出については第９章

・A1 層の水と G1 層の固体および水は常に熱交換を行う．

のモデルの検証の際に述べる．
表 1 に，式(21)～(28)の各(2) A1 層・G1 層の熱収支項に関して，A1 層，G1 層のみに生じる項を整理して示す．
表中に示した項は他の層では 0 として取り扱う．表中の各係数は以下の式で表わされる．
(36)
QH = DH Tw, A − Ta , A

６．地表面付近の熱・物質移動の定式化
6.1 散水・降雨入力（A1 層）
地表面への散水や降雨を A1 層（大気最下層）への水の注入と考え，式(19)，(24)中の Mw,sp，Qw,sp を以下のよう
に表わす．
1
(34)
M w, sp = ρw Z w
Δz A1
(35)
Qw, sp = M w, sp c p , wTw, sp

(

1

1

)

QE = ιρβ DE (qSAT (Tw, A1 ) − qA1 )(37)

Qs ,rad = (1 − ref ) S

(38)

↓

Qw,rad = L − L
↓

↑

(39)

M E = ρβ DE (qSAT (Tw, A1 ) − qA1 )
-2

(40)

ここで，DH：熱伝達率[Wm K ]，L ：大気放射[Wm-2]，

ここで， Z w ：湛水深[m]，Tw,sp：散水（降雨）水温[ºC]-1↓L ↑ ：地表面からの赤外放射[Wm-2]， S ↓ ：全天日射量[Wm-2]，QH：顕熱輸送量[Wm-2]，QE：潜熱輸送量[Wm-2]，
ref：アルベド[-]，ι ：水の気化潜熱[Jkg-1]，DE：地表面からの蒸発に関する物質伝達係数[ms-1]， β ：蒸発効率[-]，である．水平に近い舗装面への散水では，湛水深は微小で殆どのものが流出するため，湛水深（Zw）には上限（最大水層厚さ Zw,Ub）を設ける．散水量が上限を超えた場合は，現存の水層と熱交換後，系外に排水されるものとした．また数値解析上の必要性から，湛水深に微小な下限（Zw,Lb）を設定する．なお，注水による A1 層の水飽和度・
水相温度の変化は支配方程式から計算される．これらの項は A1 層以外では 0 として与える．

q,SAT(Tw,A1)：A1 層の水相温度 Tw における飽和比湿[kgkg-1]，
qA1：A1 層の比湿[kgkg-1]である．なお，地表面からの長波放射は以下の式で表現される．
(41)
L↑ = σ (Tw, A1 + 273.15)4
大気からの長波長放射 L↓ には近藤 20)の推定式を用いた．
表 1 中の fp1→p2 は，A1 層の水相と G1 層の固気相との熱交換を示しており，非定常熱伝導理論解を用いて以下のよ
うに与える．

6.2 地表面付近の熱収支（A1，G1 層における熱収支）
地表面付近で起こる短波・長波輻射，顕熱・潜熱移動，相間の熱伝達を以下のように取り扱う．図 3 に地表面付
近における各相の熱収支を示す．図中の各変数については後述する．

L↓ L↑
QH QE

S↓

(1)基本仮定
地表面付近の熱収支に関しては以下のような仮定を置く．
・湛水層は薄いため，日射（短波）は G1 層の固相に直S ↓ ⋅ ref

f
固相

時間 t

Gs

Dxz

t2
熱源

t1

QH
A1層

Zw

F

F
Gw

水相

f

f

Ga
気相

固相

Gs

Z w > Z w, Lb

風速 v

A2層

L↑
QE

S ↓ ⋅ ref

水層

接入射し，G1 層のアルベドに従い反射されるものと

L↓

S↓

F

F
Gw

水相

f

G1層

Ga G 層2

気相Z w ≤ Z w , Lb

図 3 散水時（左），非散水時（右）それぞれの地表面付近における各相熱収支・放射収支の概念図．

横分散

表 1 各層に生じるソース項．
分子拡散
地表面

Dm , a

層

縦分散

Dxx

固相

A1

図 2 大気の移流拡散効果における分散係数
の概念図．

G1

- 15 -

気相(水蒸気)Q H ，(M E )，Q s,rad ，-f s-w-f a-w

液相-Q H ，-M E
Q w,rad ，f s-w ，-Q H, ，-M E (ただしZ w <Z w,lb )f s→w = 2

λ ρ cP(Ts,G1 − Tw, A1 ) (1 − φG1 ) / ΔzG1
πΔt

(42)

7.2 地下多孔質体内部での蒸発プロセス
ある程度湿り気のある地下浅層の不飽和帯では，液体λ ρ cP
fa→w = 2(Ta,G1 − Tw, A1 ) φ (1 − Sw.G1 ) / ΔzG1 (43)πΔtここで， λ , ρ , cP ：それぞれ平均熱伝導率[Wm-1K-1]，平

水と空気が共存していることから，地下内部であっても蒸発（あるいは凝縮）
による液体水の質量変化が生じる．
この質量変化は式(27)中のMintに相当し，以下の式で表さ

均密度[kgm-3]，平均定圧比熱[Jkg-1K-1]である．
なお，A1 層と G1 層では水相同士でも熱交換が行われるが，それは支配方程式に陰に含まれている．

れる．

M int = ρ a Dm ,v ( qSAT (Tw ) − q ) Caw / Δv

(47)

-1

ここで， Dφ ,v ：間隙中の分子拡散係数[ms ]，qSAT(Tw)：水(3)蒸発効率
蒸発効率は A1 層の湛水深が下限値より大きい場合には湛水深に依存した値をとり，下限値となった場合は G1 層からの蒸発が起こるものとして，次式のように変化させる．

相温度における飽和比湿[kgkg-1]，Caw：当該格子中の水相と気相の総接触面積[m2]である．Caw の詳細については後で詳しく述べる．

Zw≧Zw,Lb
Zw<Zw,Lb

7.3 地下の水蒸気輸送の定式化

DE = Dw , β = f ( Z w )⎛ S − Sr ⎞DE = Dg , β = ⎜ w⎟⎝ 1.0 − Sr ⎠（A1 層）(44)

地下（多孔質体）での水蒸気輸送は，大気中と同様に移流拡散によって表現する．ただし，分子拡散係数およn（G1 層） (45)び分散長は多孔質体に対するものを使う．多孔質体中の分散係数については，そのスケール（実験室，フィール

ここで Dw，Dg：それぞれ水層面および土壌面（多孔質体

ドなど）に依存することが知られており，例えば

表面）での物質伝達係数[ms-1]，Sw：G1 層の水飽和率[-]，Neuman(22)により具体的な値が提案されている．詳しくは

Sr：蒸発が生じなくなる最小飽和率，n：係数である．壌

文献(23)を参照されたい．

土や砂土などの水分量と蒸発効率の関係については，近藤(19)など，いくつかの実験により提案されている．

7.4 地下の多相熱輸送
本モデルでは，稲葉ら(17)を参考に，地下の固相（岩石(4)その他の取り扱い，水相，気相それぞれの熱輸送過程を，同相間の熱伝等）
導，異相間の非定常熱伝達，水・空気の移流による熱輸

散水（降雨）時には地表面のアルベドが変化する．

送を考慮して追跡する．

本モデルでは乾燥時（Zw<Zw,Lb）と湛水時（Zw≧Zw,Lb）に

一般的な土壌や舗装面中の水や空気の流速は小さく，異なるアルベドを与えることで，その変化を表現する．
ただし，蒸発の過程で生じるアルベドの経時変化までは

固体と流体は同一温度状態にあると考えられる．しかし，追跡しない．

散水のような急激な温度変化を伴う場合や地下が高透水の場合には，固相と流体相の温度は瞬時に同温とはなら

地中への固相伝導熱および水・空気による移流拡散量

ず，過渡過程を伴うと考えられる．
相 p1 と相 p2 が接しており，各相の温度が異なる場合，については，式(14)～(16)に陰に含まれている．

DH ，DE の具体的な算定には，風速に応じた熱伝達率

接触面積付近における相間の熱伝達 Fp1,-p2 を非定常熱伝を屋内外の実験(20)から得る方法，レイノルズ数，グラス導理論解で次のように表現する．

ホフ数，ヌセルト数などから求める方法，バルク係数を用いる方法などがある．

Fp1→ p 2 = 2

７．地下における熱・物質移動の定式化

λ ρ cP(Tp1 − Tp 2 ) C p1, p 2 / Δv
πΔt

(48)

ここで，λ , ρ , cP は p1 相と p2 相の平均熱伝導率[Wm-1K-1]，
平均密度[kgm-3]，平均比熱[Jkg-1K-1]，Cp1,p2 は相 p1 と相7.1 空気・液体水の流動本モデルでは，地下多孔質体中の 2 成分（空気・水）

p2 の総接触面積[m2]である．一般に，不飽和状態であるの流体流動について，以下の一般化された 2 相ダルシー程度の水分が残っている場合，多孔質体中の液体水は固則を用いる．

相粒子を取り囲むように存在し，水相の飽和度に関わら
JJG
Kij kr , p∇Ψ p
up = −μp

ず，固相と水相の接触面積は最大となり，粒子の総表面

(46)

SAT
）で与える．気相と水相の接触面積に相当する値（ Cws

ここで，u p ：p 相の流速[ms-1]，Kij：浸透率テンソル[m2]，と水飽和度に依存すると仮定し，積については， CwSAT→skr,ｐ：相対浸透率[-]， μ ：粘性係数[Pa・s]を表す．相対

本モデルでは以下のように推定するものとする．
SAT(49)
Caw = Cws
(1 − Sw )浸透率および毛管圧力は，理論や実験結果に基づき，水飽和度との関係式が提案されている(21)．

- 16 -

計算のフローチャートを図 4 に示す．安定性を確保す
7.5 地表面での水の浸透るため，自動タイムステップや非線形反復時の収斂改善

散水や降雨により湛水した透水性の地表面からは液体機能を付け加えた．

水の浸透が生じる．このときの水の浸透はダルシー則により表現できる．地下との連続性を考え，以下の式(18)を９．モデルの検証導入し大気最下層に擬似的な毛管圧力 Pc*を与える．
Pc* = ρ w g Δz A1 ( 0.5 − S w, A1 )ここでは，開発したモデルを実観測結果(1)の再現に用(50)いる前の段階として，いくつかの動作確認を行った結果

こ こ で ， Δz A1 ： 大 気 最 下 層 の 格 子 厚 さ [m] で あ る ．
(51)
S w, A1 = Z w / Δz A1を報告する．検証内容は以下の項目である．
・大気中の風速場の表現の妥当性検討
・大気中の熱・水蒸気輸送過程の検討

大気最下層の毛管圧力を式(50)で表すことによって，地表と地下のポテンシャル差が表現され，水の浸透（あ・土壌面からの蒸発速度に関する実験結果との比較，

るいは湧出）が地下の多孔質体流れと全く同様に扱うこ・水の 2 相浸透過程，熱の移流拡散過程の検討(18)とができる．詳細については文献

・大気・地下全体の解のプロファイルの妥当性の検討を参照されたい．

9.1 大気中の風速場の表現の妥当性の検討

８．数値解法

静止した大気条件にある領域に側方境界から一定風速
基本数値解法として，同時完全陰的差分解法を用いる．

で大気を流入させ，内部に風速場が出来上がる様子を計

本モデルは，前述した 6 つの支配方程式の残差形（左辺算する．

－右辺 = 0）を，6 個の未知数に関して保存形で３次元計算領域を図 5 に示す．ここでは，一定の粗度を持つ

陰的差分展開を行い，Newton 法により反復的に解く．た

なめらかな地表面を有する 2 次元平面を仮定し，側方かだし，乱流拡散係数および風速分布の算出に関しては，らは図 6 に示す対数分布に従う風速で温度 30ºC，比湿一部陽的に計算を行う．なお，行列解法として前処理付

0.013kgkg-1 の大気を流入させる．各状態量（圧力，水飽き共役残差法（PCR 法：Pre-Conditioned Conjugate Residual

和度，比湿，3 相温度）の初期値を表-2 に，用いたパラ

Method）を用いた．

メータを表-3 に示す．計算時間は側方から風を入力し始
上部：閉境界
400m
風

42m
解析領域

大気

開始

100m

格子セットの読込

地表面

格子物性設定（間隙率・浸透率・固相物性・残留飽和率）

下部：地表面境界（圧力・飽和度・比湿・各相温度一定）

図 5 計算領域．

初期値入力（圧力，水飽和度，比湿，固相・液相・気相温度）

表 2 各状態量の初期状態．

気象条件の読込み(全天日射量・風速・比湿・大気温度)

大気地表面[Pa] 静水圧分布(地表面101325)
圧力0.4
水飽和度 [-]
1.0×10-5
-1
0.013
0.013
比湿
[kgkg ]
3030
固相温度 [ºC]3030気相温度 [ºC]3030液相温度 [ºC]∆t の短縮収斂回数？

40

NewtonLoop

残差およびヤコビアンマトリックスの生成

30

YES

マトリックスソルバーによる求解
NO

収束判定？

高度[m]

散水条件の読込み(散水量・散水温度)

NO

10m

20
表 3 計算に用いたパラメータ

YES

10

計算値出力
∆t の延長全ステップ終了？

NO

YES
終了

図 4 計算フローチャート．

0

2

4

風速[ms-1]

6

図 6 流入風速分布．

- 17 -

係数
単位
-2 -1
17～25
熱伝達率
[Wm K ]
-2 -1
-2
蒸発係数
[kgm s ] 2.0～7.0×10[-]
0.4
地表面間隙率
[-]
0.4
残留飽和度
-4
[m]
地表面粗度
1.0×10
2
-5
置換係数
[m ]
1.0×10 at 1m

めてから 10 時間とし，計算タイムステップは最大 600

は 1.0×10-4[m]とする）， κ ：カルマン定数（= 0.4）で秒まで収束条件に応じて変化させた．図 7 に領域内部中ある． 計算では，地表面温度を 32 ºC，比湿を飽和比湿

央部の断面の風速の経時変化を示した．流入開始から速状態（0.026kgkg-1）で固定し，初期状態として大気温度

やかに領域内部に入力風速に応じた対数的な風速場が形

30 ºC，比湿 0. 013kgkg-1（相対湿度 50%）として，一定成されることから，定常に近い風速場が本モデルで表現

境界流入条件の下で領域内の鉛直方向の温度，比湿分布できることがわかる．

が定常状態に近づくまで約 10 時間の計算を行った．その他の計算設定や条件は(1)と同じである．

9.2 大気中の熱・水蒸気輸送過程の検討

図 8(a)(b)は，地表面での熱伝達率および物質伝達係

鉛直方向の熱・水蒸気フラックスが一定に近づくと，数を変化させた場合の計算最終時点の温度・比湿分布を

温度・比湿変化に伴う空気の密度変化を考慮しなければ，示したものである．図 8(a)(b)より地表での熱伝達率や物温度および比湿の鉛直分布は理論的には高度に対して対

質伝達係数の変化に伴い傾きが変化し，片対数でほぼ直数的関係となる．しかし，温度・比湿の変化に伴う密度

線的な鉛直分布が得られている．これにより，モデル計

変化を考慮した場合には地表付近は対数的関係から外れ算が温度・比湿場を適切に表現していることがわかる．

ると考えられ，その状態が表現可能か否かについて検討

なお，地表面から高さ 7m 付近以上で片対数直線からわする．

ずかにはずれているがこれは上部境界設定の影響であると考えられる．今回の計算条件では，熱伝達率を(1)密度一定条件および乱流拡散効果による計算

20Wm-2K-1，物質伝達係数を 5.0×10-2ms-1 として与えた場理論解と比較するため，計算上，式(32)の分散長を全合定常解の温度・比湿分布（実線）をよく再現しており，て 0 とし， Dm ,a に乱流拡散係数（熱および水蒸気共に同

これらの値は近藤(24）による交換速度 CHU に換算した値

じ値）を与えて計算を行う．一般に安定度が中立に近いに近い値である．

環境では，乱流拡散係数は次式から算出される．
Dm , a = κ u * z(52)

u* = u κ ln ( z z0 )(2)密度変化を考慮した熱・水蒸気輸送の計算(53)

ここで，u*：摩擦速度[ms-1]，z0：地表面粗度（ここで
x=30m

0.02s
0.05s
0.3s

x=60m

じ条件下において，式(32)を使い，縦分散・横分散モデ
ルを用いて乱流輸送過程の計算を試みる．対数則が成り

x=90m立つような条件下では式(52) ，(53)よりそのオーダーを推測することが可能であり，

8

γ T = γ L = κ 2 ln(z z 0 )(54)と表すことができる．このとき，分散長は，

0
0

1

2

0

1

2

0

1

0.002 < ϖ T = ϖ L < 0.2

2

10

DH = 15.0
DH = 17.0
DH = 20.0
DH = 25.0
対数則

0.1

29.8

30.6
温度[ºC]

31.4

10

1

0.1
0.01

10

γ = 0 .1

Dw = 2.0
Dw = 3.0
Dw = 5.0
Dw = 7.0

γ = 0.01
γ = 0.02
γ = 0.05

γ = 0.01
γ = 0 .1
γ = 0.005
対数則

対数則

対数則
高度[m]

1

場合の温度および比湿分布の計算を行った．計算では縦

高度[m]

10

(55)

とおくことができる．この範囲で分散係数を変化させた

図 7 領域内部断面の風速の経時変化．x は風上側からの水平方向の位置を示す．

1

0.1
29.5

0.02

比湿 [kgkg-1 ]

図 8 (a) 熱伝達率[Wm-1K-1]を変化させた場合の鉛直温度分布，(b) 物質伝達係数[×10-2ms-1]を変化させた場合の鉛直比湿分布の計算結果．

図 9- 18 -

高度[m]1
2
風速[m/s]

高度[m]

高度[m]

16

次に，

温度や比湿変化に伴う空気の密度変化も考慮し，(1)と同30.5

温度[ºC]

31.5

1

0.1
0.01

0.02

比湿 [kgkg-1 ]分散長[-]を変化させた場合の(a) 鉛直温度分布の計算結果，(b) 鉛直比湿分布の計算
結果．

0.03

分散長および横分散長は同じ値を与えた．なお，熱伝達

えず，ビーカー内の砂層と表面からの蒸発のみを対象とするため，大気層については 1 層のみ設定する．計算に率および物質伝達係数は(1)での最適値を用いる．

用いた物性値およびパラメータを表 4 に，各状態量の初図 9(a)(b)は，分散長を変化させた場合の温度および比湿の鉛直分布の計算結果である．いずれの場合も片対数期状態を表 5 に示す．用いた砂の物性のうち，間隙率・上で直線に近く，分散効果が大きい程その傾きが大きく

固相密度・温度・比湿に関しては，神島の実験により測定されている値を用い，その他の値に関しては代表的ななり（即ち，鉛直方向に一様になり），分散効果を小さく値を用いた．格子設定では，格子の断面積および高さがすると，特に密度変化の大きな地表付近では直線的分布実験に使用した 500mlビーカーとほぼ等しくなるようにから外れる様子が表現されている．これらは，理論的に設定した．計算に用いた砂の毛管圧力および相対浸透率予想される現象と整合的である．

にはVan Genuchtenの式(21)を用い，式に必要な係数については相原(26)の実験結果を用いた．地表面における蒸発効9.3 多孔質体表面からの蒸発に関する実験結果との比較

率は式(45)より，水飽和度の関数として図 11 に示す曲線

土壌面や保水性舗装などの多孔質体からの蒸発・輸送で与えた．多孔質体表面からの水蒸気の物質伝達係数Dg

モデルの妥当性を検証するため，
ここでは，多孔質体（砂）については，0.015～03ms-1 まで変化させて計算を行った．

からの蒸発速度の測定を行った神島 (25) の室内実験を対計算タイムステップは最大 600 秒まで収束条件に応じて象に，本モデルによる再現性を検討する．

変化させた．

神 島 (25) は， 砂を ビ ー カー

(1) 神 島 (25) の実 験 概 要

図 12 は多孔質体から(3)計算結果と実験結果の比較（500ml）に充填し，水で飽和させたあと，恒温恒湿槽の物質伝達係数 Dg を 0.3，0.15，0.015kgm-3 とした場合のに入れ，数時間毎に全質量の計測を行っている．この作計算結果と実験結果それぞれの蒸発量－経過時間の関係業を数日間行うことで，蒸発量－経過時間の関係を求めを示したグラフである．図 12 より，物質伝達係数を 0.15ている．

ms-1 とした場合，本モデルによる計算結果と神島による実験結果はよく一致した．

本モデルによる計算に用いた 1 次元格(2)計算設定

実験結果を最もよく再現する物質伝達率 0.15 ms-1 の妥子設定を図 10 に示す．実験は恒温恒湿槽で行われている

当性を確認するため，バルク式との比較を試みる．

ため，本モデル計算においても大気中での乱流輸送は考地表面からの蒸発は，バルク式を用いると以下の式で
表される．

1.2

大気層（1層）最大：1×10 -2 m
最小：2×10 -3 m

0.12m
（17層）

Dg =0.3200

0.6

蒸発量[g]

蒸発効率 β [-]

1
0.8

0.4
0.2

S r = 0.2, n = 0.3

0
0

0.2

0.4

7.6×10 -2 m
（1層）

0.6

0.8

図 11

Dg =0.0150蒸発効率と水飽和度の関係曲線．

[Wm-2K-1]

物質伝達係数

0.3，0.15，0.015[ms ][-]／浸透率
-1

4.0

-1

2

-2

-1

密度[kgm ]／比熱[Jkg K ]／熱伝導率[Wm K ]残留飽和率 [-]
Van Genuchtenの係数

0

-2

-1

比湿 [×10 kgkg ]固相温度 [ºC]
気相温度 [ºC]
水相温度 [ºC]40

60
80 100
経過時間 [hour]

120 1400.01

比湿[kgkg-1 ]
0.02

0.03

0

水飽和度[-]
0.5

1.0

0.41／1.0×10
2650／840／2.1
0.2
a : 14.4，n : 1.43

表 5 各状態量の初期状態．
空気圧力 [Pa]水飽和度 [-]

20

-11

[m ]

-1

0

図 12 蒸発量の時間変化．実線は計算結果を，白丸は実験結果(26)を示す．

深度[m]

熱伝達率

Dg =0.15

1

表 4 計算に用いた物性およびパラメータ．

間隙率

100

水飽和度 S w [-]図 10 1 次元格子．

-3

実験結果

大気（1層）
土壌
静水圧分布（地表面：1013.25hPa）
1.0×10-5
1.98

0.99

30

30

30

30

30

30

0.04

0.08

2.4
0.12

1時間後
10時間後
1日後

図 13 (a) 比湿の鉛直分布，(b)水飽和度の鉛直分布の計算結果．

- 19 -

M E = ρ a β CH ua ( qSAT − q )カー内部の比湿は上方の大気との接触部分が最も減少し，(56)

水飽和度に関しても時間と共に上方から全体的に減少している様子が計算により表わされている．

ここでCH：バルク輸送係数[-]， u a ：平均風速[m/s]である．式(37)と式(56)を比較すると，本モデル中における物

9.4 地下の熱伝導・水の浸透・熱輸送に関する検討

質伝達係数は上式中のバルク輸送係数 C H u a に相当す

(1)地下の熱伝導計算と理論解との比較

る．今回計算対象としたような面積の狭い範囲（ビーカーの断面積，直径約 0.1m）からの蒸発では，バルク輸送

まず地下多孔質体を単一固体と単純化した時の熱伝導に係数 C H u a は通常の裸地などにおける値よりも大きくついて理論解との比較を試みる．

なることが知られており，次式で表すことができる(24)．
1.82U
C H ua =( Re = 10 ~ 2500)
Re0.6
0.83U
C H ua =( Re = 2500 ~ 15000)
Re0.5

本モデルで

は 3 相非平衡状態の温度分離型計算を基本としているが，一方が一定温度に保たれた一次元固体中の熱伝導過程の理論解は，以下の式で表わされる．

(57)⎛ x ρ cP ⎞T − TL= 1 − erf ⎜⎜ 2 λ t ⎟⎟TH − TL⎝⎠

また，X：風速方向の物体の寸法[m]として，(60)

ここで，TL：媒体の初期温度[K]，TH：熱源温度[K]，λ ：Re = X u a ν(58)

媒体の熱伝導率[Wm-1K-1]である．
計算格子は 1 次元格子とし，地表面付近の層厚を 0.01m，である．ここで， Re ：レイノルズ数，ν ：空気の動粘性係数[m2s-1]である．上式(57)，(58)に実験条件を当ては地下深部ほど層厚が大きくなるように設定した．最大深めるとRe = 6000～7000 程度（恒温恒湿室の風速を 1ms-1

度は 2m である．表 6 に使用したパラメータを示す．初程度として換算）と見積もられ， CH ua = 10−2 ~ 10−3 程期状態として，格子内部の温度を 30 ºC とした．なお，度となる．本モデルでは，土壌面からの蒸発は，土壌最固相物性として花崗岩の値を使用している．熱源は G1

上層の固相以外の面積（全面積A×間隙率 φ ×飽和度Sw）

層に与えるものとし，熱源温度は 50ºC とした．計算時間
は 3 時間である．なお，ここでは媒体中の内部蒸発およとして扱っているため，

び地表層からの蒸発は生じないように設定した．本モデ

JJG
Dg ≈ CH ua S wφ

(59)

ルによる計算結果および理論解による鉛直温度分布を図と表すことができ，Dg は 10 ～10 ms 程度となる．今回，14(a)(b)に示す．ここでは，3600 秒，10000 秒の結果を示

-1

-2

-1

(26)

神島

-1

の実験結果と最もよく一致した 0.15ms という値

表 7 9.4 (2)の計算に用いたパラメータ．

はバルク式と比較しても妥当な範囲であると考えられる．

間隙率 [-]

図 13(a)(b)は．計算されたビーカー内の水飽和度およ

浸透率 [m ]-3
密度 [kgm ]-1 -1
比熱 [Jkg K ]-2 -1
熱伝導率 [Wm K ]残留飽和率 [-]
Van Genuchtenの係数

び比湿分布の時間変化のグラフである．この結果より，ビーカー内の砂上面から蒸発が生じることにより，ビー
表 6 9.4 (1)の計算に用いたパラメータ．
0.001
0.1
2000
840
1.7
30
50

-3

固相密度 [kgm ]-1 -1
固相比熱 [Jkg K ]-2 -1
固相熱伝導率 [Wm K ]初期温度 [ºC]
熱源温度 [ºC]

深度 [m]

200.0

温度[ºC]
40

60 20

温度[ºC]
40

60

00

60分後
30分後

350分後
180分後
350分後

本モデル
GETFLOWS

1.0
1

Case1
Ca se2

(a)
04
0.4

06
0.6

0 8 1.010.8

飽和度[-]

1.5
2.0

180分後

60分後

5
0.5

180分後

計算値
理論解

(a) 3600s

5分後

5分後

5分後

0.5
1.0

1.0×10-12
2000
840
1.7
0.4
a : 4.37，n : 0.37

60分後

深度 [m]

間隙率 [-]
残留飽和度 [-]

0.4

2

(b)
0.4

0.6

0.8 1.0 20

飽和度[-]

4060

温度[ºC](c)80

図 15 (a)本モデルによる計算結果と GETFLOWS による計算結果の比較．(b)浸透水の温度差（地表層温度
Case1：20 ºC，Case2：80 ºC）による浸透速度の比較． (c)Case2 における温度分布（気相温度）の時間変化．

(b) 10000s

図 14 (a)3600 秒後，(b)10000 秒後の温度分布．

- 20 -

している．どちらの時間帯においても，本モデルによる（Case2，G2～Gng 層の初期温度は 20ºC）の水飽和度分布

計算結果は理論解とよく一致していることが分かる．

の比較結果である．図 15(b)より，高温の水を浸透させた
場合，温度上昇によって水の粘性が低下することにより浸透速度が速くなっていることが認められる．図 15(c)(2)地下への水の浸透に関する既存モデルとの比較
多孔質体中の 2 相流動モデルの妥当性を確認するためは水の浸透に伴う温度分布（3 相同一温度）の変化を示

1 次元での浸透計算を行い，既存の数理モデル

した図である．図 15(b)(c)より，液体水の浸透速度よりも温度の移流拡散（伝導）速度の方が遅いことが分かる．

(15), (16））との比較を行った．

（GETFLOWS

格子設定では最大深度を 2m とし，格子厚さは最小

これは，温度の高い液体水が，周りの固相や液相と熱交0.01m，最大 0.1m とした．ここでは地下への水の浸透を換をしながら浸透するためである．これらの計算（熱の計算対象としているため，大気層（A1 層）は 1 層のみと移流拡散を伴う 2 相浸透過程）は比較できる理論解や実し，この層の温度・圧力を一定とした．使用したパラメ

験値がないため定量的な比較は今後の課題であるが，定ータを表 7 に示す．初期状態として，間隙圧力は大気圧

性的には妥当な結果と考えられる．

-1

状態とし，水飽和度は 0.4，比湿は 0.013kgkg （相対湿9.5 大気・地下－熱・水輸送プロセス全体の動作確認

度 60％），媒体温度 20 ºC を与えた．以上の条件を設定し，G1 層（地表面層）の飽和度を常に 0.8 として水の浸透過

現実的な気象条件変化（気温，風速，日射など）を入程の計算を行った．計算時間刻みは 10 秒一定とし，全計力した場合，本モデルが対象とする全領域（大気・地表

算時間は 3 時間である．
便宜上，
上記設定での計算を Case面・地下）における計算状態量のプロファイルと水文学

1 とする．図 15(a)に本モデルによる計算結果と参照モデ的知見などの整合性を確認するため全体的な動作確認を

ル（GETFLOWS）による水飽和度の鉛直分布を示す．両行った．
計算領域は水平方向約 100m，鉛直方向に大気約 16.7m，モデルの計算結果はほぼ一致し，本モデルの地下流動計

地下 10m の範囲とする．入力気象条件（図 16）は東京

算は基本的に妥当なものと考えられる．

都における 2008 年 8 月の快晴日（かつ日照時間 7 時間以
上を対象）の平均気象データである．なお，風速に関し(3)水の浸透に伴う熱輸送に関する検討ては抽出データ全体の時間平均値（2.4 ms-1）を一定としとして同様の浸透計算を行い，浸透水の温度の影響や相て与えた．
図 16 の各気象要素の昼夜の時間変化を繰り返

間の熱交換が適切に表現されているかの検討を行った．

し与え，約 6 日間の計算を行った．さらに，地表面が土図 15(b)は，全体温度 20 ºC とした場合（Case1）と地

壌面のケース（CaseA）と舗装面のケース（CaseB）を設

表第 1 層（G1 層）の各相温度を 80 ºC 一定とした場合定した．両ケースに共通で用いた計算パラメータを表-8

気温
日射量

400

20

比湿

0
0

1

時間[day]10

表 8

両ケースに共通の計
算パラメータ．

熱伝達率
蒸発係数(土壌)
置換係数
分散係数γ
地表面粗度
アルベド
蒸発効率の係数-2 -1
10[Wm K ]
-2 -1
-2[kgm s ]2.9×102-5[m ]1.0×10 at 1m-2[-]
1.0×10-4[m]
1.0×10[-]0.3 (0.2)‐4
Nv ： 0.4
Sr ： 10

図 16 入力気象条件．

40

1.0×10-16
0.1
0.1
1000
840
1.0

1.0×10-12
.0.4
0.1
2300
840
1.0

1.0×10-12
0.4
0.1
2300
840
1.0

表 10 大気・地下における各状態量の初期状態．

20
0.6

600

Rn

400

QE

-1
比湿[kgkg ]温度（3相）[ºC]

10cm
50cm

0.4
0.2

G

200
0

土壌 (Case A)
舗装面 (Case B)
地表面
地下
地表面
地下
静水圧分布（地表面：1013.25hPa）
0.9
0.6
0.4
0.4
(1.0-5)0.02450.01440.01440.01440.0144
30
20
20
20
20
大気

圧力[hPa]
水飽和度[-]

0cm
20cm

水飽和度[-]

密度[kgm-3]
比熱[Jkg-1K-1]
熱伝導率[Wm-1K-1]

1.0×10-12
0.6
0.1
2300
840
1.0

30

フラックス[Wm -2 ]浸透率[m2]間隙率[-]
残留飽和率[-]

舗装面 (Case B)
地表面層
地下

30
1.5m

表 9 計算に使用した地下の熱・水理物性値．
土壌 (Case A)
地表面層
地下

0.05m

35
気温[ºC]

30

温度[ºC]

800

気温[ºC]，湿度[×10-3kgkg-1]

全天日射量[Wm-2]

Case 2 として，Case 1 の設定下で浸透水の水温を 80 ºC

QH

2

経過時間[day]46

図 17 Case A における計算結果．(上)気温変化，（上中）
地表面地中温度変化，(下中)地下の水飽和度変化，(下)熱量フラックス変化．地表面地中温度変化の凡例は水飽和度変化のグラフに同じ．

- 21 -

に，各ケースの地下水理・熱物性を表 9 に，また，各状ら下方への熱流を正とする）を見ても，CaseB では日中

態量の初期状態を表 10 に示す．
計算時間刻みは最大 3600に舗装面から最大 80Wm-2 程度の熱量を地下へ輸送して秒とし，収束状況に応じて自動的に変化させた．

おり，地下温度の上昇が CaseA より大きい．図 19 は，

図 17 は CaseA（土壌面）の気温，地表面・地中温度，

地下 1m から地上 3m までの温度の鉛直プロファイルを

地下の水飽和度，各種熱量フラックスの時間変化を示し示している．日中(12:00)，CaseA では地表面の温度およたものである．図 17 より，地中温度は地表面温度に遅れ

び地下の温度は CaseB よりかなり低い状態が保たれる．

を伴いながら日変化し，その振幅は地下深部になるほど

また，夜間(0:00)のプロファイルを見ると，気温分布の差

小さく，深さ 0.5m 程度で日変化はほとんど表れていなはほとんど見られないが，CaseB では CaseA と比較してい様子が分かる．また，水飽和度も日変化を伴いながら地中温度が 5～7 ºC 程度上昇する（蓄熱する）結果とな

徐々に蒸発により減少している．各種熱量フラックスでっている．これらの結果より，舗装面では，土壌面と比は，潜熱輸送量が最大で 400Wm-2 程度生じており，顕熱

較すると，潜熱輸送量がないため，
地表面温度が上昇し，-2

輸送量が 100Wm 程度である．また，気温では高度によ

地下温度や日中の大気温度が上昇する様子が計算された．

り多少の温度差が生じている様子が見られる．なお，高これらの結果は一般的な水文学的知見や既存研究 (27)

さ 1.5m に見られる夜間の微小な温度のゆれは入力気温と定性的に整合しており，本モデルは全体として水・熱のゆれによるものである．

輸送の物理過程の定量的追跡に適用性があるものと考え

図 18 では，CaseA と CaseB での地表面温度および熱られる．

量フラックスの差（顕熱輸送量，地中伝導熱）の時間変化（計算開始後 3 日～6 日目）を示している（なお CaseB

１０．結論の潜熱輸送は 0 であるため図示しない）．図 18 より，舗
装面（CaseB）の場合，地表面温度が最大 55 ºC 程度にま

本研究では，都市域の熱環境改善策としての散水実用化で上昇しているのに対し，土壌面（CaseA）の場合は，に向けた散水効果の定量評価を目的として，大気・地下を最大でも 40 ºC 程度であり，潜熱の効果が表現されてい連成した水・熱輸送数値モデルを開発した．その要点は以る．また，顕熱輸送量を見ると，CaseB の方が日中最大下のようにまとめられる．

で 50Wm-2 程度大きい．また，地中伝導熱（地表面層か(1) 全領域を 3 相（気・液・固相），2 相 3 成分流体系（水，空気，水蒸気）として一般化し，大気中の物質・熱Case B

40
30
Case A3
4

5

6

経過時間[day]

乱流輸送，地表面付近の放射熱収支・顕熱・潜熱移ΔQH

50

動，地下の気液 2 相流動・熱輸送を含めたモデルの基礎方程式および数値解析手法について詳細を述べ

0

た．
(2)

-50 ΔG345

経過時間[day]6

開発したモデルの動作確認，解の妥当性を検証する
ため，大気中での熱・水（水蒸気）輸送，地下多孔質
体表面（砂）からの蒸発，地下多孔質体での水・空気

図 18 (左)CaseA と CaseB における地表面温度の変化，(右) CaseA および CaseB の顕熱輸送量(QH) と 地 中 伝 導 熱 (G) の 差 (Δ ={CaseB-CaseA})を示す．

の 2 相流れと熱輸送過程，岩石を模擬した熱伝導過程，
および全体的な解のプロファイルの整合性，に関して計算し理論解や実験結果と比較を行った．その結果，既存の文献値に整合する物理パラメータを与えるこ

0:00とで，各物理プロセスが適切に表現できることが示さ

12:00

れた．
本モデルは，ここで報告した以外にも多数のケース計算

2を行っており，2 次元領域において 6 つの状態量を比較的高度(深度)[m]温度[ºC]

50

熱量フラックス差[Wm-2 ]100

安定に計算できることが認められてきた．今後は，各種物理プロセスについて，さらなる実験との比較などの定量的な検証を行うと共に，3 次元問題への適用性，メソスケー
ル大気モデルとの結合の可能性についても検討していきたい．続報では，著者らがこれまでに行ってきた屋外実験

0

結果(1)に本モデルを適用し，観測結果の再現性について論
-1

CaseA
CaseB
20 30

20

じる予定である．

40

温度[ºC]

図 19 CaseA，CaseB の温度分布．

- 22 -

参考文献(1) 吉岡真弓，中川康一，登坂博行，ヒートアイランド
現象緩和のための都市散水効果に関する実験的検討，土木学会論文集 G，64-1(2008)，pp.68-77．
(2) 例えば，木内豪，神田学，栗城稔，小林裕明：都市内散水による気温緩和効果の現地観測，水工学論文
集，38(1994)，pp.381-386．
(3) 例えば，狩野学，手計太一，木内豪，榊茂之，山田正，打ち水の効果に関する社会実験と数値計算を用いた検証，水工学論文集，48(2004)，pp.193-198．
(4) 土屋修一，加藤琢磨，手計太一，山田正，打ち水による市街地の熱環境緩和効果，水工学論文集，49-1(2005)，pp.367-372．
(5) 小作好明，春日井哲夫，橋本一郎，古賀睦，高根澤
浩二，保水性舗装に散水した場合の気温・湿度・黒球温度への効果，平成 19 年度都土木技術センター
年報(2007)，pp.153-164．
(6) K-1 model developers，k-1 coupled model (MIROC)
descriotion, K-1 technical report,1(2004), H. Hasumi
and S. Emori (eds.)，Center of Climate System Reseach，
University of Tokyo．
(7) Grell. G., J. Dudhia and D. Stauffer，A Description of
the Fifth-Generation Penn Statue/NCAR, Mesoscale
Model (MM5), NCAR Tech. Note(1994) ，
NCAR/TN-398+STR，National Center for AtmosphericResearch.
(8) 例えば，ヴタンカ，足永靖信，浅枝隆，都市大気環
境層の乱流モデリング 都市建築計画における都市気候予測システムの開発 その 1，日本建築学会
計画系論文集，536(2000)，pp.95-99．
(9) 例えば，近藤裕明，亀掛川幸浩，玄地裕，井原智彦，
大橋唯太，東海林孝幸，都市気候・ビルエネルギー
連成モデルによるヒートアイランド対策の評価，日本ヒートアイランド学会論文集，1(2006)，pp.53-57．
(10) Harbaugh, A.W., Banta, E.R., Hill, M.C., andMcDonald, M.G. ， MODFLOW-2000, the U.S.
Geological Survey modular ground-water model- User
guide to modularization concepts and the Ground-Water
Flow Process(2000)，U.S. Geological Survey Open-File
Report 00-92, 121p．
(11) Pruess, Karsten ， TOUGH2 –A general purpose
numerical model simulator for multiphases fluid andheat flow(1991) ， Lawrence Berkeley Preport
LBL-29400．
(12) Diersh, H.J.G.，WASY Software FEFLOW 5.4 User’s
Mannual(2008)．
(13) Šiｍůnek, J., M. Th. Van Genuchten and M. Šejna，The
HYDRUS Software Package for Simulating Two-abd
Three-Dimensional Movement of Water, Heat and- 23 -

Solutes in Variably-Saturated Media, Technial Manual,
Version 1.0(200&), PC Progress, Prague ， Czech
Republic，pp.241．
(14) 萩島理，谷本潤，片山忠久，大原賢志，改良・建築－都市－土壌連成系モデル（AUSSSM）による都市
高温化の構造解析 第 1 報 モデルの理論，日本建築学会計画系論文集，550(2001)，pp.79-86．
(15) 登坂博行，伊藤彰，田中将希，岩井卓，自然水理系の流体・熱移動統合モデリングの試み その 1．流体・熱移動の新しい定式化と実験的検討，地下水学
会誌，41-3(1999)，pp.147-158．
(16) 登坂博行，伊藤彰，岩井卓，自然水理系流体・熱移
動統合モデリングの試み その 2． フィールドシ
ミュレーションによる適用性の検討，地下水学会誌，
41-3(1999)，pp.159-176．
(17) 稲葉 薫，登坂博行，平井 卓，吉岡真弓，地圏水・熱循環系の統合モデルの構築－熱・流体移動の定式
化 と 解 析 的 検 討 － ， 水 工 学 論 文 集 ， 51(2007) ，pp.469-474．
(18) 登坂博行，地圏水循環の数理 流域水環境の解析法(2006)，東京大学出版会．
(19) 近藤純正，水環境の気象学－地表面の水収支・熱収支－(1994)，朝倉出版．
(20) 萩島理，谷本潤，成田健一，都市表面の対流熱伝達率に関する既往研究のレビュー，水文・水資源学会
誌，17-5(2004)，pp.536-554．
(21) 例えば， M. Th. Van Genuchten，A Closed-form
Equation for Predicting the Hydraulic Conductivity of
Unsaturated Soils，Soil. Sci. Soc. Am. J，44(1980)，
pp.892-898．
(22) Neuman, S.P ， Universal scaling of hydraulicconductivities and dispersivities in geological media,，
Water Resources Resarch,，26(1990)，pp．1749-1758．
(23) 日本地下水学会 原位置トレーサー試験に関するワーキンググループ編，地下水のトレーサー試験
地下水の動きを知る(2009)，技法堂出版．
(24) 近藤純正，地表面に近い大気の科学(2000)，東京大学出版会．
(25) 神島要，陸域水循環系における蒸発散現象のモデル化に関する研究，東京大学大学院工学系研究科地球システム工学科修士論文(2004)．
(26) 相原貴，不均質・不飽和地層中における溶質移流分
散現象のモデル化についての研究，東京大学工学部
システム創成学科 E＆E コース卒業論文(2006)．
(27) 例えば，近藤裕昭，応用気象シリーズ 8 人間空間の気象学(2001)，朝倉書店．
(Received August 31, 2009, Accepted November 20, 2009)