国立防災科学技術セソター研究報告第22号

1979年10月

556．14：556，332．6

鉛直降雨浸透の有限要素法による一解法
言十算誤差の発生とその除去法

大

倉

博＊

国立防災科学技術セソター

A皿A皿a1ysis

MethoaofVertica1Rai皿fal1Inmtratio皿by＿Derivationof

FiniteE1ements

Ca1culative趾mr

and

its

E1imination一

By

1Himshi

OhkuraM肋ηα1R・∫・〃I・んαηオθげ・1一〃∫α∫伽P〃τ1ε〃・・，J砂α・AbstractTheequationoffundamenta1equation，based
continuity，of

one

on

Richards

theoryofcapillary；1一￡（・（ψ）；1）
was

solvedweighted
θis
is

after

residuals

volumetriccapi1lary

over

being
for

water

potential；

transformed

space

and

is

total

In

order

to

app1y

the

and

under

the

ground−water

intosimultaneous

the丘nite

content；〔sんpotentia1anddimensionaいn丘1trationdifferencetime；K（ψ）ispotential；

and2equation

in

the

same

surface，the（1）equationsbythemethodmethodfortime．In

pemiabi1ity

forunsaturated

iscoordinate．

verticalsaturatedandequation（1）istheunsaturatedchangedequation（2）istransformed

intopopular

the

space
two

above，〃isposition；thesubscripts

region，theareinterva1of

time

（・）

simultaneous

notationshavingmatrices．Whensolutionsofstep；superscript
one

there

means

equationsasunkmwnん箒1

（・）

time

step；subscript

subscript

areappears

wetti㎎front

c1earlyinerror

causesa1ossequation（3）have

a

a

large

of

into

（鮒・・o穿ノ2）鮒・一c娑㌦1
Inabove，
soil；ψ

reg1ons㍑一￡（凪ψ）；1）・
then

of

vectors；and

which

the

means

notationthehavinganalytic

of

water

balance．
So

as

to

eliminatetheerror，equation（1）isθ，intO

＊第4研究部計測研究室
一145一

directlytransformed

withremaining

国立防災科学技術セソター研究報告第22号

1979年10月

α鮒2附 ム肌＝一（θ訂Lθ1皿）

（4）

〃

But

this

1．緒

equation（4）cannotbeSo

so−vedequation（3）issolutionissmallersolved

than

by

the

reason

being

mixed

hundredth

of

after
one

of

its

with

the

non−convergency−

equations（4），and

error

in

case

of

the

error

of

equations（3）only．

言

雨水の地中への浸透と地表面からの蒸発は，地表表層の土の性質に強く影響を受ける。表層の土は大部分不飽和であるから，土中の不飽和帯の水の挙動の解明は，農地の灌概と排水の問題，河川の流出率と流出係数の決定に有力な情報を与える．

一般に，地下水の滴養は地下水面上の不飽和帯を通して行なわれるから，地下水の酒養機構の解明には不飽和帯の水の挙動の解明が必要になる．

土の力学的な強度は土中の合水率と間隙水圧に密接に関係し，土構造物・白然斜面の崩壌機構の解明も不飽和浸透を考慮しなければ恋ち、ない．

透水係数が毛管ポテソシャルの関数になるため，不飽和帯において流れを支配する運動方
程式が非線形偏微分方程式になる．この非線形偏微分方程式の解析解を得ることは一般に困難であるため，電子計算機を用いた数値解析を行なって近似解を求めざるを得ない．

筆者はかつて降雨浸透の数値解析における有限要素法の有用性を述べ，事前に境界条件の定まらない境界とかなり強い非線形性を有する場合の解析法を定常鉛直二次元の解析例を用いて述べた（大倉，1977）．

非定常の場合は浸潤面の通過に伴う毛管ポテソシャルと含水率の急激な変化のため，数値解析に誤差が生じる．本研究ではこの誤差の発生と除去法について，非定常鉛直一次元の解
析例を用いて述べる．

2．

浸透のモデル

降雨の鉛直浸透の基本式を導くにあたり，次の仮定（1）〜5））を設ける．
土は一様，等方性である．

土は時間・空間的に一定な物理的性質をもつ．

水は時間・空間的に一定な接触角，粘度，密度をもつ．
温度勾配による水の移動は無視する．

土中の空気の移動と空気の圧力変化が土中水に与える影響は考えない．
このとき，Richards（1931）の式
9＝一1（（ψ）9rad（ん）（1）

が成り立つ．ここで，qは比流束，K（ψ）は不飽和透水係数，ψは毛管ポテンシャル（テソシヨソともいう），んは総ポテンシャル（全水頭）である．（1）式と連続の方程式より鉛直

一146一

鉛直降雨浸透の有限要素法による一解法一大倉

一次元の浸透の基本式

繁一ヱ（凪ψ）芸）（・）が得られる．ここで，θは体積含水率（以下含水率という），チは時間，2は鉛直方向座標（上向き正）である．（2）式は3個の未知数θ，ψ，んを持つため解を得るには（2）式の他に
θ，ψ，ん問の2本の式が必要になる．ここに，仮定（6）川8））を追加する．

6）総ポテンシャル乃は位置の座標2と毛管ポテソシャルψの和

乃＝ψ十・

（3）

であらわされる（一般に右辺第2項はgを重力加速度としてg2になるが本報告では
圧力を水柱高で表わすため単にzとした）．ψは不飽和域ではψ〈0，白由水面上で
ψ＝O，飽和域でψ≧0となる．また，白由水面下ではψが水頭，カが全水頭をあらわすものとする．

7）含水率θは毛管ポテ：／シャルψの関数
θ＝θ（ψ）＝θ（ん一2）（4）で，一般にθはψに対してヒステリシスを有する．

8）透水係数Kは飽和浸透に対する透水係数K8と比透水係数K『の積であらわされる．

K（θ）＝KW（θ）（5）

ここで，K『は含水率θの一価関数とする．

境界条件は，境界からの流入量が既定されている境界ではγを単位面積あたりの流入量として∂ん凪θ（ψ））亙＝W（6）となり，総ポテソシャルが既定されている境界では乃。を既定総ポテソシャルとして

々＝机

（7）

となる．ここでWは，境界が上端のとき！V＝1，下端のときW＝一1になる．また，降雨
強度が土の最終浸透能を越えないときは降雨強度をRとすると，地表で，W：1，γ＝Rであるから∂ん

凪θ（ψ））τ一R（8）となる．降雨が最終浸透能を越えるときは地表に湛水が生じることがある．湛水があるときは，加を湛水深とすると，（7）式より地表で

ゐ＝加

（9）

となり，湛水がないときは（8）式になる．しかし，本報告では，降雨強度が最終浸透能を越える場合は考えないものとする．

一147一

国立防災科学技術セソター研究報告第22号

1979年10月

3．重みつき残差法
土中水の鉛直一次元の浸透は（2），（3），（4），（6），（7）式の境界値問題を解くことにより定量的に評価できる．しかし，この境界値問題は，Kがθの関数になり，θが乃の関数になるため非線形になり，特殊な場合を除いて解析解を得ることが困難になる．このため，境界値問題を時間・空間に離散化して，非線形連立方程式に変換した後，この連立方程式を解く

ことにより境界値問題の近似解を得る．この離散化にあたり，（2）式をθを未知数とする連立方程式に離散化する方法（K1ute，1952）と，（2）式をψまたはんを未知数とする連立方程式に離散化する方法がある．θを未矢口数とする方法は仮定2）と仮定3）より飽和域でθの変化が零になり，解析が不可能になる．このため飽和域と不飽和域を一括して扱い難た

く不飽和域のみを解析の対象としている例（Rubin，1967：山村・久楽，1972）が多い．ψ
またはんを未知数とする方法は上記の欠点を持たないため多くの解析例（Per・en・，1977：

赤井・大西・西垣，1977）があるが，解析領域に浸潤面が明確に現われる場合に，水収支の誤差が生じる（後述）．

本報告では，んを未知数として，空間については重みつき残差法，時問については差分法（中央差分）を用いて境界値問題を非線形連立方程式に変換する．

重みつき残差法（Fin1ayson，1972：大島・毛利・中川，1976）は以下のように説明される．

工を微分演算子として，領域「において定義された次の微分方程式を考える．

〃ト0

（10）

この式に任意関数Wを掛けて領域11について積分を行なうと

∫

（11）

W・五［例∂11＝Orとなる．逆に任意関数Wに対して（11）が成立するならぱ関数んは（10）式の解になる．
ここで，関数んを適当な有限関数系炉1，仰，…，伽，…，で展開できるものと仮定し，次式であらわす．

乃＝乃冊ψ冊（12）（以下，表の表現の簡便化を計るため，特にことわらないかぎり，右脚添字を添えた量は，

脚添字が1個のものはベクトル，2個のものはテソソルをあらわし，さらに，式中の同一項に同じ脚添字をもつ量が繰返し現われるときは，その添字について総和をとるものとする）．

同様に，Wもまた適当な有限個の関数ω・，ω・，…，ω冊，…からなる任意の線形結合W＝〃冊ω冊と表わし，（11）式に代入すると

刈、吻・五［W・一・（13）となる．ここで〃刎が任意定数であるから

一148一

鉛直降雨浸透の有限要素法による一解法一大倉

∫

（14）

ω冊・工［ん刎ψ刎］∂r＝0rが得られる．これは未定係数加に関する〃元の代数方程式であり，これを解いて（12）式に代入すると（10）式の近似解んが求められる．関数系伽を試行関数，関数系ω殉を重み関数といい，ω冊を特に伽に等しくおくものをGa1erkin法という．このときんを比較関数，伽を基底関数ともいう．

4．非線形偏微分方程式の離散化（2）式の未知数んを基底関数伽とその未定係数んの有限級数

ん＝ん例（≠畑刎（・）（15）で近似する．甲刎（z）は以下の手続きで求める．

ZZ

図1に示すように2軸上の解析流域（鉛直1

n＋2

n＋2

次元）の全体を11としrを適当な間隔で分割

n＋1

n＋1

する．分割された流域を要素と呼び通し番号をnnつける．要素の境界点を節点と呼び，通し番号n−1をつける．要素および節点の通し番号を下方よ

n＿1

n−2

n−2

り1から順番につけるものとする．また節点〃0の座標を砺とする．このとき，〃番目の要素（以下要素〃という）は，〃番目の節点（以Fi＆1Ei，m、、t

下節点〃という）と節点〃十1との間の区間と

2一砺十1
砺一砺十1
2−2例一1hn．1hnhn，1hdivid、，fu，d，m，nt，1f、。、ti．n

仰・nd・0mP・・i・ion
FEM．

定義する．伽を似z）＝01怜

図1要素分割と基底関数ρκと比較関数ん
fun・ti㎝んfo・

2＞z柵・1
砺十1≧Z＞砺
（16）砺≧2≧2柵十1

砺一砺一1
＝0

2炉1＞2

と定義する．

このとき，ψ柵（2冊）＝1，g㎜（z冊）＝0（舳キ〃）となり，未定係数伽は節点〃上のの値になる．また，（15）式は節点上の乃の値を直線で連結する直線近似になる．θとKも乃と同様に

θ＝θ例（τ）甲刎（2）（17）

K＝K物（ま畑冊（2）（18）

で近似する．

（2）式の右辺を左辺に移項した後，基底関数伽（z）を乗じて全流域「について積分する、

一149一

国立防災科学技術セソター研究報告第22号

1979年10月

1、腸一去（唯）／仰・・一・（・・）

上式に（15），（16），（17），（18）式を代入して整理すると
∂θ㎜（

∂舳

）

∂ま

一Q刎十α舳力肌（オ）＝O

（20）

ここで

伽一一1、ψ繁ψ（・1）

ムー一1、榊ゐ

（・・）

Q冊＝

（23）

ここで，舳＝〃一1，〃，〃十1以外のα舳と∂舳は零になる．ρ物は節点〃が流域内の流入量が既定される境界上にあるとき現われ，篶がその単位面積あたりの流入量になる．

をタイムステップ，右肩添字ゴをタイムステップ回数として，（20）式を時間について中央差分を用いて離散化するとム肌

漱／2鮒1＝一 （θ；，才」θ㍍）十α
〃

（24）

ここで
1K・（K二）｛・1／2＋K8（K；．、）｛十漱ρ＝一／21K5（Kζ。、）｛十

22例一2冊＿1

1K5（K二）也十

十一

ノ2＋K5（κ；）｛十1／2

22冊十1−z刎

ノ2＋Ks（K二一。）｛十1／2α鮎＝α閉＝一一2z冊一2冊＿1
1

1∂冊珊＝一（z冊。r2犯）十一（2r2炉1）2
2
6冊例＿1＝6冊＿1他＝0

（25）
（26）
（27）
（28）

1
（凪）ぜ十1／2＝一（K・（θ三十1）十K「（θ三））2

θ乏＝θ（桝一2例）

（29）

（30）

（22）式を計算すると（27）式と（28）式は
11∂舳＝一（・柵。、一zπ）十一（2ザ2冊一1）
3
3∂仰例．、J弼、、例一⊥（、仰一、物．、）

（31）

（・・）

6
となるが，（27）式と（28）式に変更している．これは以下の理由による．タイムステップ

クからゴ十1における要素〃一1の含水量の変化は，θが節点〃上の値θ柵を直線で結ぶ直線近似であらわされるから，

一150一

国立防災科学技術セソター研究報告
1

1

2

2

第22号

1979年10月

一（砺一2冊一1）（θ三士1一θ三一、）十一（砺一紅1）（θ乏十」θ乏）

1
3

1
6

＝一（砺一Z冊一1）（θ三士1一θ真一。）十一（2ザZ仰一1）（θ三十」θ三）

1
6

1
3

＋一（砺一2冊一1）（θ三士1一θ乏一、）十一（zr2刎一、）（θ三十」θ差）

となる．

（33）

（27），（28）式のf犬りに（31），（32）式を用いると，（24）式の右辺第一項は単位時間における要素〃一1と要素〃の含水量の減少量の離散化による節点〃への配分量
批（砺一仙）（lll・1㌧）・÷（砺一仙）（1■一刈米（舳一砺）（1ガー1l）・÷（舳一砺）（1い1）／

（34）

（35）

の和である．

（24）式の左辺は単位時間の節点〃から節点〃一1と節点〃十1への浸透流量であるから，（24）式は（節点〃から隣接する節点への流出量）

二（節点〃に隣接する要素より配分された合水量の減少量）
十（節点〃が境界にあるときの境界からの領域への流入量）（36）を意味する．

今，浸潤前線が節点〃と節点〃十1の間に到達しθ冊。。が増加中とする．このとき（θ測
一θ三。、）／

は大きな値になる．（24）式の〃番目の方程式に（θ崩一θ美。、）／が含まれるから，（24）式の〃番目の方程式の右辺第1項は急に減少する．このときこの方程式を満たすためには，左辺の未知数雌1，ん島十1，炸1の係数の値がそれぞれ，負，正，負であるから脈｝

が急増または附1が急減または晦1が急増しなければならない．もし，炸1が急増または炸1が急減すると，浸潤面の到達直前に総ポテソシャルが急変することになり矛盾が生じる．このため（31），（32）式を（27），（28）式と変更した．（33）式から，変更された（27），（28）式は各要素の含水量の変化を漏らさず，各要素の境界にある節点に配分していることがわかる．

非線形連立方程式（24）式は右辺に附1の項を持ち，きわめて収束性が悪いため一般に（20）式をムー（影）刎繁一軌・伽一1一一・（37）と変形した後に時間について離散化する．

（枠oチ）鮒一〇芳狐・ρ1

（38）

一151一

鉛直降雨浸透の有限要素法による一解法一大倉炸去（砺一㎞）（者）1＋ψ・去（㎞一砺）（劣ズ（39）c榊；：班鷲＝0（40）（劣ズー（影）、（、榊，（41）

ψ三＝ん三一砺（38）式は右辺の6鮒2が鮒（42）の関数であるにψもかかわらず，（24）式に比べ収束1生が良い．こ漱

れは，左辺のが直接燃

ノ2．

〃燃1により含水率の変化分

ψ苫にフィードバックされるためであ

＿

F

T

る．しかし，々の時間変化率が大きいと（37）式の解は計灘のため水収支関係を満足しなく咋〕

＿＿L＿

，

なる。この誤差の発生を示すため，図2に仮の外1

θ一ψの関係を示す．タイムステップゴのψ差

・θ1

丁

S

θtち1

。

F

θ・旨1θ

図2仮のヒステリシスカーブと数値計算の誤差の発生

時間後にψ三十1になったとする．図2の1∬1曲見積られる節一点〃による含水率の変化は

1ll（ψ三＝〃一z柵）が（38）式を解くことにより，

線が節点〃のθ一ψの関係を示す．（38）式でC

Fig．2 Schematica1 hysteresis
derivation of error．

炉一（景）1㌦一・1）

curve

and

（43）

となる炉は図中の1炸」1lであらわされる直線肝は・α上の点・（1門ψヂ）における接線・㌘に平／寸な点・を通る直線であ／，〃の縦軸に対する傾きが（劣）∴な一ている・ところ／l，1一ψ関係を用いてψ㍑算出される含水率は
θ三十1となり，時間ステップゴからゴ十1への実際の含水率の変化は

猟十ノ2＝θ1＋」θ三（44）となる．1θ美十1／2と1θ二｛十1／2との差が累積して水収支の誤差が生じる．さらに，（38）式は見積られた貯留量の変化に誤差が生じているので，（38）式の解炸においても誤差が生じている．特に，シミュレーショソの領域に浸潤面が現われると，この誤差が大きくなる．浸澗

面の前後でψの差が大きく，浸潤面の通過に伴ない，通過地点のψが大きく変化するからである．この誤差を小さくするにはすると，を十分に短かくすればよい．しかし，を短かくに反比例して同一時間を数値解析するのに要するタイムステップが増大するため，〃を変えずに誤差を小さくする方法を提案する．
収束性の良い一般的な（38）式を解き肘を求め，得られた解が水収支を維持するかどうか

一152一

鉛直降雨浸透の有限要素法による一解法一大倉検査する．すなわち，全ての閉について桝

から附1を算出し，εを水収支の許容誤差とし，（景ズ（鮒一11）一（卯一11）1・1（45）

であるかどうかを検査する．検査に不合格な全ての舳に対して，（38）式の左辺第1項中の（景）肌㌢を（晋）帆にも1した後に・時問について離散化を行な1ことによ／（・・）
式を水収支を維持するように修正する．修正された（38）式のr求解」，r検査」，r方程式の修正」を，未修正の舳に対して新たに（45）式の検査の不合格が生じなくなるまで繰返し
誤差の小さい解を得る．

（38）式の修正は以下のようにすると簡単に行なわれる．収束性に欠ける（24）式と一般的な（38）式を混合した非線形連立方程式

λ鮒2附1＝域十1／2
λ二な／2＝α二な／2＋（1一α例）

（46）

c鮒空

（47）

〃
必丸1／2

1現十1／2＝α㎜．

．

6美丸1／2（θけ1一θ㌦）十（1一α刎）

〃

ん㌦

（48）

〃

を作る．（47），（48）式においては舳についての積和はとらない．α榊はOまたは1の値を取るが，全てのα肌がOならば（46）式は（38）式になり，全てのα㎜が1ならば（46）
式は（24）式になる．

全てのα伽を0とし（46）式を解く，得られた解から附1を算出し（45）式を検査する．

検査に不合格な舳に対するα㎜に1を代入し（46）式を修正する．（46）式の求解，（45）式の検査，（46）式の修正を新らたに1が代入されるα例が生じなくなるまで繰返す．

5．非線形連立方程式の解法∂θ∂ψ

Kがθの関数，θ，一がん一zの関数になるから（46）式は，係数λ鮒2と定数項風十が鮒の関数となり，非線形連立方程式になる．このため，始めに鮒の繰返しにより逐次炊ρを仮定し線形計算の修正を行ない，（47）式の解を得る．さらに得られた解が水収支を十分な精度で維持するかどうかを検査する．検査に合格したときは次のタイムステップに進む．不合格のときは合格するまで，水収支を維持するように（46）式の修正と修正された方程式を解くことを反復する．以下解法のアルゴリズムを述べるが，特に断わらないかぎり

順を追って実行するものとする（図3参照，本文および図3におけるr：＝」は右辺の値を
左辺の変数に代入することを意味する）．

操作1初期値を設定し，タイムステップのサフィックスゴに1を代入する．
操作2全ての〃に対して，1タイムステップ前の鳩を〃三十の仮の値〃十に代入する．

操作3始めに比較的収束性の良い（38）式を解くために，全ての〃に対してα冊＝0と
する．

一153一

国立防災科学技術セソター研究報告第22号

1979年10月

操作4連立方程式の修正反復回数

STAR
三：二1

のカウソターんに1を代入する．

ん宗1二片、1

操作5

んが連立方程式の許容反復

dπ1・0

回数んより小さいかまたは等し

for・allNO

いかどうかを検査する，不合格な

≦ム

YES

らば解が得られないものとして計0．1。。1．t．teメ燃帥仰O・〃．＾i

点1・わ1

SOLVE

を用

三↓1朋

鳩と仮に定めた〃十S舳mNEOUS三十脆一︒伽ん二十1・脇、算を打切る．

操作6

π

点1・1

＾

JUDGl・O
k：・

数項B三十1ρを計算する．

十1

η：：1

NO

≦η

ES

操作7操作6で作成した連立1次

4η1・

YES

方程式（46）式＝風十STOP
SOLUTll
SOLUTlONNO

YES

いて，（46）式の係数λ鮒2と定

＾サ／2燃

NO

土41
后蜆 一κ1・ε1

NO

三：・且 ・1

EOuAT10NS

財1ん

＝

NO

醐一（部棚1
ノ2

（49）

ぴ几：・1

YES

〜1・ハ十1

JUDG：＝1を解く．

操作8操作7で求めたん臭十1と操NO

作6で用いた〃十1との差の絶対

N0

JUDG・

YES
〜≧五

YES

値を求め，差の絶対値が許容誤差

STO

SlMUL＾TlON

ε厄よりも小さいかどうかを検査す

る．全ての〃が検査に合格するならば現時点のα㎜に対する連立

WAS、、、：TOWCOMPLETEO

NODE

図3解法アルゴリズムの流れ図
Fig．3

Flow

chart

of

argorithm

tosolve．

方程式（水収支を満足するとはかぎらない）が解けたものとして操作11に飛び越す．
操作9（46）式の収束解が得られるように〃十を〃十1：＝〃十1＋伽（50）と修正する．伽を求める方法は（大倉，1977）による．その理論は付録に加えてある．

操作10連立方程式の修正反復回数カウソターんを1繰上げる．操作5に飛び越す．
操作11α冊＝Oなる全てのに〃対して，εを許容誤差として1（1ガー1l）一（影ズ（舳）1・1（51）を検査する．不合格ならばα刎に1を代入する．ここで，α腕をOから1に変えることにより（46）式は水収支をより良く満たす（24）式に白動的に修正される．

操作12操作11において新たに1が代入されたα冊がないことを検査する（∫σD0＝0を検査すればよい）．合格ならばタイムステップゴ十1の水収支を満足する解が得られたものとする．検査が不合格ならば連立1次方程式の修正回数のカウソターκを1繰上
げ操作5に飛び越す．
一ユ54一

鉛直降雨浸透の有限要素法による一解法一大倉

操作13

んをシミュレート終了時間に対応するタイムステップとして，シミュレート終了時間に達している（4≧ゴ侶）かどうかを検査する．合格すれば，シミュレート完了として計算を終了する．不合格ならばゴを1っ繰上げ次のタイムステップの解を求めるた
めに操作2に飛び越す．

6．

数値計算例および考察

モデルの検証のために，実験と数値計算との対比を行なった．しかし，検証実験とは別個に土の物理定数，たとえば，Kτとθ，θとψとの関係を算定しておらず，この対比は物理定数の同定ともいえる．一般的な（38）式（以後，一般方式という）を用いた数値計算結果と本報告で新しく提案された（47）式（以後，新方式という）を用いた数値計算結果を検証

実験と比較する．物理定数の同定は一般方式の結果が検証実験に一致するようにし，この物理定数を新方式の数値計算にも流用したため，一般方式の方が新方式より検証実験に良く一致している．しかし，一般方式による数値計算例は水収支関係を満足せず，一般方式にょる
降雨浸透のシミュレーショソは誤差を伴う危険11生を持っ．

1）検証実験
検証実験は科学技術庁国立防災科学技術センターの大型降雨実験施設において行なわれた（富永，1978）．図4に示される4槽に区切られている土槽（各槽の寸法はタテXヨコ×タカ

サ，240cm×240cm×240cm）を作成し，各槽に供試土としてそれぞれ粗砂（粒径0．5−2．0mm），細砂（粒径O・0−0・5mm），関東ローム土，鹿沼土を充填した（本報告のモデルの検証は細砂のみで行なわれている）．土槽には，地下水面が土槽上端から下方170cmに位置す

るようにドレイソパイプが設定されている．なお，各層の底部に層圧40cmの砕石層を置き，ドレイソパイプの土槽側開孔端付近に生じる水平方向の浸透流が供試土中の鉛直浸透に影響しないようにした．測定項目

は，降雨量，毛管ポテソシャル，

240om

体積含水率，土の電気比抵抗，地240㎝

下水の水位，地下水の流出量の6項目である（図5）が・数値計算

FlNE

SAND

K＾NUMAPUM■CE

LO州

240㎝
240、、

240c岬

／

（αo一α5・・〕！

COARSEノ穆

およびその検証に用いるものは・。、α5瑚、〕
次の4項目である．

竃

／穆

／

件の実現を確認する．
ii）

テソシオメーターによる、亨1
トT1PPlNG；昌

あらかじめ設定されるが設定条
図4

BuOXET

TYPE

FLOWMETER

実験用土槽

Fig・4Largescalegroundmodels（1ysimeter）．

一155一

ε
卜

i）転倒枡雨凱！1降雨1＼ノ

量：降雨強度・降雨継続時問は／

国立防災科学技術セソター研究報告第22号

1979年ユ0月
240om

土中各部の毛管ポテソシャル：○

テソシオメーターのポーラスカ

0
○当

40

山0
］」

雨継続中は8点につきおおよそ

5分，その他はテンショソの変

←

く⊃

o

OzLお

u］
ロー
■o一

一

」o

oo〜

←u＿
o一山o﹂崔山oo

o

o

ooo＝一

160

工
れている．測定時間間隔は，降邑

○庄

E

彗
・m，一60・m，一80・m，一120；1・・

cm，一160cmの8点に埋設さ由

L］O＝〇望ωΣ削鴫

ミ；・α＝Φ︺80

L］←

o口⊃暮o訪一ωω山L﹈O＝ocm，一10cm，＿20cm，一40圧

］
←
山
Σ

庄コ匝1鴫

山虻o

ップは供試土の表面より一5

o

山0

Z

oo
く
L］

Σ

」一

O

○つ

TOTOFLOWFLOW

GR州D榊TER

200
E
o
て。、〕寸
o
寸

240

化の緩急により変えている．

LEVEL

O
ト

RUBBLE

一．．一
BAsE

二L⊆o

↑9
r「9

図5計測器配置図
Fig．5

Locationofmeasuring

instruments．

iii）中性子水分計による土中

各部における体積含水率：測定位置は供試土の表面より一10cm，一20cm，一40cm，一60

cm，一80cm，一100cm，一120cm，一160cm，一180cmであり，測定時間間隔はおおよそ30分以上である．

iV）

ドレイソパイプから流出する地下水流出量：転倒枡流量計により自動計測する．

本報告の数値計算は細砂の場合に限られている．この検証実験は，土槽に砂細を充填し，一年間養生した後に行なわれた．図6に実験の結果（降雨量と地下水流出量（A），深さごとのテソションの時間変化（B），時間一深さ空間における総ポテンシャルの等値線（C））を示す．

30分間50mm／hの降雨が継続した後，30分間休止することを3回繰返す間歌的な雨を設定して実験を行なった．しかし，第1波，第2波，第3波の雨は，ピーク値がそれぞれ41．5

mm／h，55．5mm／h，50．5mm／hとなり，ばらついている．総降雨量は72．8mm，降雨継続
期間1．5時間に対する平均降雨強度は48．5mm／hである．

地下水流出は，降雨開始後240分に始まり，350分でピーク値8．8mm／hに達し，780
分までの地下水の流出高は35．7mmになる．地下水流出の時間変化が滑らかにならないのは，測定に用いた転倒枡流量計の転倒枡（11）が大きすぎるためである．

含水率とテソショソの値は正の相関関係がある．よって，図6（B）のテソショソの曲線の■ド降は浸澗過程，上昇は排水過程をあらわす．

深さ，一20cm，一40cm，一60cm，一80cm，一120cmのテソションのそれぞれ12分，90分，140分，170分，252分における急激な増加は浸潤前線の到達をあらわす．

浸潤面が完全に通過した，一5cm，一10cm−20cmでは，問歌的な雨に応じてテソショ
ソの値が変化し，その振幅は深くなるにつれて小さくなる．

第3波の降雨が終了するまでに浸潤面が通過している部分では降雨の終了とほとんど同時に排水過程が始まっている．この排水によるテンショソの回復は表層に近いほど早い．（こ一156一

鉛直降雨浸透の有限要素法による一解法一大倉
8◆RnINF≡≡；＝1LL

INτEN5IτY

0】RuN−0FF

・I・・1

lNτENS−TY

01一

＝lL

Oτ01

日L

R

＝■INFnLL

RuN−0F1＝

01

ギ
0N

（A）

竃。
…N

s

0
8

60

120

1日0

240
240

300
300

368
368

T1HE

RI＝I1NFnLL

＾

・鰭

RND

RUN−0FF

O一
◇一

一5C［：
一6日CH：

420
420

488
488

540
540

688
688

660
660

128
128

760
76

【I−llNUTE】〔EXPERlNENT〕

工F工NE△＝ 一10CH〕
午工 一80CH〕

SRND

50［［／H＊8．5H

十： 一2日C同〕ス：一12日Ch】

3−SOUl＝IRES〕

×： 一49CH〕
Z：一16日Ch〕

・

蓄

（B）一〇
蕎甲

0

68

128

180

2

0

300

368

T1刷…

TENSION〔EXPERユNENT〕

工FINE

SRN1〕

τlnE

○轟

o＝

一

60

・一

128

一一）

三一一＝一

188

I

2ω

一一一

−28

488

540

600

660

720

788

一HlNuτE，

50［［／H＊0．5H

3−SOU∩RES

178，05，22〕

540

1
660

：HlNUτE，

380 一u一
368 ■6一
428 ，一1480●一一一50CH一リmER

600

一，一

一一〇

720

780

一1a6CHmER！

（C）書S■5㏄［一H前Eぺ＝

丁0TI＝lL

POTENTII＝lL〔EXPERI［ENT〕

工FINE

SI＝1ND

50［N／H＊0．5H

3−SQuRRES〕

図6 実験結果（降雨と地下水流出（A）．テソショソ（B）．総ポテソシャルの等値線（C））
Fig．6 ResultS of experiment．
Rainfallandground−waterrun−o舐（A）．

Tension（B）．

Contours

of

total

potential（C）・の一5cmのテソショソの時問変化を後述する数値計算結果と比較すると，実験の一5cmのテソショソの値は約6cm−H・0の測定誤差を持つと推定される．一5cmのテソショソの曲線を縦方向に6cm

H．O上方に平行移動すれば一5cmのテンショソの回復が最も早くなる）．

総ポテソシャルの等値線（C）において，浸潤両の降下は原点から斜め右下方に階段状に伸びる等値線の間隔の狭い部分であらわされる．この浸澗両の降下をあらわす部分が階段状に
一157一

国立防災科学技術セソター研究報告
．；「1

①：

口．二じ

一10〔：！＝■」

・■L−ml1o．30＾〔■：臼ロー一｛）〔・日8⊂一「］Zト1舳／］

0、．二P

→［

第22号

ユ1一≡lo

日

8，10

0，20

タ1

〜．4a

日．50

柵

彗ポヘ

曇い

1義＼

1

＼

毒予。・十
＼

鰯

・＼＼轟∴、；

8．38＾＝
L2竈∩ j芸1丁亨臼こ 1・

一．ヨ；ヨ〔［一Y［．舳［］

ド㍉

1979年10月；一一〇 〇・ニニ ロ・j臼 O■5
VOLu旧Rl［HnTERCONTENτ

3・13

0

図7テソショソと含水率との関係（実験結果）

Fi・・7R・1・ti…hi・b・tw・・・…ti・…dも∴0．19 82日 日．38。藺、・日
VOLu［ETRlC〕RTERCONTENT

O．5日

図8実験によるヒステリシスの主ループの推定

Fig．8E・tim・ti…fm．i．1。。。。fh。。t、、、、i、

・・1・m・t・i・w・t・・…t・・t・

f・・m・。p。。im。。・。。。。1t。．

なるのは，時間のサソプリソグに対して空問（深さ）方向のサソプリングが荒いためである．

図7に検証実験のテソションと含水量との関係を示す．土中各部の同一時刻のθとψを
時刻を変えてθ一ψ平面にプロットしている．θ一ψの関係はヒステリシスを有し，図中の曲線はヒステリシスの主ループ（後述）である．

2）物理定数
図8に数値計算に用いるヒステリシスの主ループと，ほぼ一週問間隔で行なわれた5回の実験の同一時刻同一地点のθとψのプロットを示す．これらの実験の（降雨強度）×（降雨継続時間）はそれぞれl15mm／hx6h，30mm／hx6h，50mm／h，x6h，30mm／h×1hの後70mm／hx1h・50mm／hxo・5hが1時間のサイクルで3回繰返す間歌的な雨である．図中の曲線がプロットされた点から推定したヒステリシス主ルーブである．前述の5回の実験の降雨強度はいずれも細砂の最終浸透能（102〜103mm／h）に比べ僅小なため，排水過程の主ループを推定するθとψがプロットされない．このため排水過程のヒステリシス主ル＿プ
ー158一

鉛直降雨浸透の有限要素法による一解法一大倉はMua1emによる砂質土と砂のヒステリシス主ループ（Mu・lem，1974；Mua1em，1977）を参考に推定した．ヒステリシスの主ループが定められるとMualemモデル（Mua1em，1974）
を用いて，主ループ間の走査線が決定され，θ一ψと（∂θ／∂ψ）一ψの関係が得られる．ただし，飽和域ではMua1emモデルが用いられず，一般に
θ＝θ刎十（〃β十α）ψ

∂θ

万＝

（52）

β十α

（53）

となる．ここで，θ砒は飽和域でψ＝0における体積含水率，〃は土の空隙率，αは圧密などによる土の間隙率の変化，βは水の体積弾性率である．しかし，1章の仮定2）と3）よ
りα＝o／cmH・O，β＝o／cmH・oとなる．θ刎と〃は共に，θ砒＝o．38，〃＝o．38とした．

飽和透水係数Ksは検証実験の供試土をサソプリソグして室内透水試験を行ない，K∫＝1．4
cm／minの値を得た．しかし，数値実験と検証実験との対比によりKs＝6．0・m／minとした．

κとθとの関係はImayの3乗則（Irmay，1954）

一（岩）31・ん
＝1

（54）

θ＞θ仙

（55）

を用いた．ここで，θ。は最小容水量である．降雨開始直前には，前回の実験よりほぼ1週問経過し，地下水の流出強度がほとんど零になる．この状態では懸水帯の水はほとんど移動

せず，この懸水帯のθは最小容小量になると考えた．ただし，地表付近のθは蒸発による乾燥のため最小容水量より小さい値になると考え，深さ一60cmのθの値を用いることにした．各実験の降雨開始直前の深さ一60cmのθはO．06〜0．08になるが，θ田＝O．08と推定した．このときのθとK『の関係を図6に示す．
日

日．1日

a．28

3）要素分割と境界条件とタイムステップ

Q

試土の最下端（以下，最下端という）までを解≧甲析領域とする．検証実験において表面が土槽上…。

端より約5cm沈下していた．この沈下した表昌o

面を鉛直座標2の基準点とすると，供試土の最具r

下端（土槽上端より一200cm）はz＝一195cm…。

になる．また，土槽下底の砕石層における浸透

ゼ

藺．30

日．40

」 9
供試土の表面（以下，地表面という）から供 とFの圧力損失が無視出来るものとすれば，最下端
藺の総ポテソシャルは一165cm

O．1匂H20になる．よ

図9って1境界条件は〃＝一165．mH，O0，20

v0LU旺TRlc岬TER

0，39

8．48

c0NTENTθ

比透水係数K『の体種含水率θによる変化。・。＝＿195cm，Fi・・9V・「i・ti…f1・1・ti・・h・d…li・…一ductivityκwith volumetric water（56）
・・nt・・tθ．
一159一

国立防災科学技術セソター研究報告第22号

1979年10月

降雨強度Rが最終浸透能より小さいから肋（θ）□：R
∂2

Ks

at

2＝0cm

（57）

となる．

節点間隔とタイムステップ間隔は電子計算機の解析プログラムの占有容量が小さく，かつ演算時間が短かくなるように検証実験の精度を損なわない限り大きくすることにした、テソ

シオメーターの最小設置間隔が5cmなので，地表面から最下端まで節点間隔〃＝5cmで節点を設けた．このとき，節点総数は40，要素総数は39になる．

検証実験のテ：／シオメーターの最小読取間隔が5分であることからタイムステップ＝5分とした．計算不安定をおこさないためには，一般に」z／〃が解析する現象の伝播速度より大きいことが必要である．浸潤前線の降下速度は（降雨強度）x（降雨継続時間）が50

mm／hx6hの実験においてo．7cm／分であった．検証実験の浸潤前線の降下速度は，降雨

強度が50mm／hであるが間歌的な降雨のためo．7cm／分を越えない．ル＝5cm，〃＝5分とすると」2／＝1cm／分となり浸潤前線の降下速度より大きくなる．

非線形連立方程式を解くために必要な収束判定許容誤差はε九＝0．0002cm

H．O，水収支許容誤差ε＝0．O0004とした．

4）初期条件ψ，θ，Kと∂θ／∂ψの初期値が必要であるが，K，∂θ／∂ψは乎とθより計算で求められるので，gとθの初期値が必要である．θとgの初期値は検証実験の実験開始直前の測定値を
用いた．

5）

数値計算結果4章において定義したム㎜を変更せずに（31），（32）式を用いた場合の計算結果を図10に示す．浸潤面の通過に伴ない上流（上方）に隣接する節点のθが急激に増大するため，

ゲパ1ポlllザ1為1㌶1
・1

1選

一一r

「

紳．一．

二芸・

ノ麦

念・・

一

4

瀞豊そ帯冷≒一トー｛一一暑一一≡≡一一一…」

○

臼8

12冒

了ENS工0N

図10

r胴

2ω

〔SI［ULRT工0N〕

3日呂

368

42臼

4腱

T］旺川］NUW〕

l

FlNE

SRN口

S8〔N！H＊8山5H

川■■■一ト■」」

548

Eε臼

6臼臼

3−SOURRES

フ20

実験結果に反する数値計算例
ExampleTheofresultss1mulation

withoutwhich

violates

modify1r■g

experimenta1results．

coe冊cientん肌of

一ユ60一

simultaneous

ア冒O

78・8王5・22〕

連立1次方程式の6舳を変更しない場合
Fig．10■equations一

鉛直降雨浸透の有限要素法による一・・解法一・大倉午RnlNF日LL

寵富

十τ0τnL−NτEN511Y

凹RllN−OFF1冊NSIτ、

R日一NPPLLOτ0τ∩LRUN−OFF

一・←一十一→一→i一←一・1一一十一十一■←一←一←一トi＋i一←一←一

一一十一→一→一i←一←→一i←一←一．一一←一十一十

o一
ミー
；
．

（A）

on

o
自NN

8
8

60
60

12日
12日

180180

2＾0
248

3藺0
388

360360τIHE

Rn−NF∩LLΣnND

RUN−0FF

＾フn
4Z0

＾Rn
480

、540藺

目日0
6日0

R60
669

7フ0
720

ラR0
760

−HlNUτE：〔SI［ULl＝IT〕0N〕

〔FユNE

①一

一5C［1

△一

一1日C「〕

十：く！〉〔

一60Ch〕

午一80Ch〕

×一一128C［〕

SRND

一2日Ch〕

58［H／H＊8．5H

×一

3−SOURRES〕

一4日C［〕

Z一一160C［〕

1

（B）畜胃
涙1

÷き岩≒■書喘」姜一呈←｛串≒出粉060TENSION

128

2ω

303

〔Sユ［ULRTION〕

180

〔FINE

・60

・20

S∩N口

58［［1H＊0．5H

．r］nE

τ1nE

o

o

68

12藺

）

166

248

380

−nINuτE】

・80

舳

6日0

3−SOunRES

660

→

フ20

仰

I7b，85，22〕−H！NUτE〕
368

』20

480

548

668

660

720

78日

一50㎝一HmER．舳．岬丁。。ノ（C）蜆

一150CH−HRτE尺

T0TRLPOTENTユPLlS川ULRT10N〕旧NESRN口50［［／H・0．5H3−SOURRES〕

図11 一般方式による数値計算結果
降雨と地下水流出（A）．テソショソ（B），総ポテソシャルの等値線（C）
Fig．11

Simulation

results

Rainfal1andground−waterofpopularmethod．

run−o冊（A）．

Tension（B）一Contours

of

total

poter■t1al（C）．

ψが一時的に減少することがわかる．

ム㎜を変更して，一般的な（38）式を用いた一般方式の計算結果を図11に示し，本報告で提案された（46）式を用いた新方式の計算結果を図12と図13に示す．一般方式の結果が新方式の結果より検証実験に良く一致している．新方式では浸潤前線の降下が一般方式のものより早くなっている．これは以下の理由による．湿潤過程において，ヒステリシスの主ルー
一161一

第22号

国立防災科学技術セソター研究報告
十R

SW＝nNF自LL

十一0TPL−NTENS一τY

凹R州一0FFlNτENS川■■1 □I一 一
■
■一

Oτ0τ

一

1

R＾INFnLL＝lL■1

1979年10月

RUN−0FF■→一一←→一十一十一←一←一←→一十一十一←一トー←一←十一十

1一一トー←一ト→一十一十一←一←→一一．一十一ト→

、
．

（A）嚢

o
ω

0

60

120

RR1NFRLL

RND

180

RUN−0FF

①一
◇一

＾

240

308

360

τ一日E

−nINuτE＝｛SINuLRTION〕

一5Ch〕

△一

460

〔FユNE

一18CH〕

十一

手ト8日㎝〕

一68C［〕

420

540

SRN0

一20CH〕

660

50NNlH＊a．5H
×一Xl−120C［〕．

600

Z

728

780

3−SQURRES〕

一49〔H〕

一160㎝〕

1

：了

（B）看雫

・蹄一喘絆一一一一一書→
0

6〜

一20

180

2』O

300

560

420

τ1NE

TENSION

o

〕

548

600

660

7三0

700

〔SI［ULnTION〕 工FINE SnN口 58NH／H＊0．5H 3−SOUnRES 178．85．22〕

6日

o

480

｛11▲NuT三〕

120

168

））

24

丁ユHE

＝11＝NuTE；

300

360

42⑤

46a

540

688

660

72〜

ンUO

一SOCn−H日1E尺

・100C［一H∩TER

（C）；。
ヨ2

．1。日㎝．。町。r一ノTOTRL

POTENTIRL〔S1［ULRT＝0N〕

｛FINE

S∩N口

5aNN／H＊6．5H

3−SOunRES〕

図12 新方式による数値計算結果
降雨と地下水流出（A），テソショソ（B），総ポテソシャルの等値線（C）
Fig．12

Simulation
Rainfallresu1ts

and

of

new

gromd−water

method．

rm・off（A）．

Tension（B）．

Contours

of

total

potential（C）．

プまたは走査線が図2のSS1曲線であらわされ，節点〃が図2のF点にあり，浸潤前線が要素〃を通過中としよう．新方式と一般方式の炸の変化に対する含水率の変化の見積りは新方式の方が小さくなる．このとき新方式（（46）式）の解炸が一般方式（（38）式）の解より大きくなるから新方式のθ三十1も一般方式のそれよりも大きくなり，新方式の（K；）也十1／2

も一般方式のものより大きくなる．節点〃一1と節点〃の全ポテソシャルの差鳩十」脈1−162一

鉛直降雨浸透の有限要素法による一解法一大倉は新方式によるものが一般方式より大きく（節

竈

。．、口

。．、。

o，30

0．

日

o．50

点〃一1にはまだ浸潤面が到達していない），さらに節点〃一1と節点〃との問の透水係数1／2
×（K岳（K；）｛十1（D〔一10C［1×L−6臼〔［〕

Xl−1ηC［〕＾一一28⊂［〕

φ旦

一8臼c門〕

ヌ正■雄〔；「〕■トー

一且臼C［〕

争＝

一「88⊂［〕Y工一15日〔［〕2＋岬凪．、）｛十1／2）も新方式が一般方

式より大きくなる．ゆえに，新方式のタイムス

テップクでの節点〃から節点〃一1への浸透量が多くなり，浸潤前線の降下が早くなる．

ド

新方式の解は検証実験の以下の現象i），ii），iii）を満す．

i）浸潤前線の通過によりテソシヨ1／が急激に増加する．

ii）第3波の降雨開始以前に浸潤面が完全に通過した部分において，テンシヨソは間歌的な雨に応じて変動し，その振幅は深くなるにつれてノ1・さくなる．

。＼

iii）第3波の降雨が終了するまでに浸潤面が通過している部分では，降雨の終了とほとんど同時に排水過程が始まり，さらに，排水によるテソションの回復は表層に近いほど早い．

、、

臼

ε．l

日．20

VOLuHETR−C

しかし，新方式の解の以下の定量的な5項目

竈．ヨa

−PτER

8．』8

8．58

CONTENT

図13テソショソと含水量との関係が検証実験と良い一致をみない．a）浸潤前線の降下速度．b）浸潤面通過後のテソシヨソの値，

o

Y

工＼

新方式による数値計算結果
Fig．13

Re1ationshipwaterC）地下水流出の開始時刻．d）ピーク流出時刻．

betweensuctioncontent（simulation

of

and
new

m・th・d）．

e）ピーク流出強度．

図14に新方式と一般方式の水収支の誤差の時間変化をあらわす．降雨強度が最終浸透能を越えないとき，蒸発を無視すれば，各時刻において（降雨量）十（土中の初期水分量）＝（土中の水分量）十（地下水流出量）

が成り立たなければならない．1〕を降雨量，1）を土中の水分量の初期水分量との差，Gを地下水流出量とし，水収支の誤差ERRORをERROR：lP＿1）ルτ＿Gl（58）で定義する．このとき，各時刻においてERROR：0になる．数値解析においては丸めの誤差などのためERRORは零にならないが，零に近い値にならなければならない．ERRORの大小により，数値解析のモデルと解法アルゴリズムを評価できる．図14において新方式のERRORはf＝5分（第1ステップ）で0．00087mm，オ＝780分（第156ステップ）で一163一

国立防災科学技術セソター研究報告

凹N［W［ETH0口

第22号

1979年10月①FOPULPR［ETHOD

①
ヒ

ト8

6日

128

188

248

棚

36臼

Tユ［E

428

棚

5畑

600

醐

728

78a

｛［1NUTE〕

岬TERBnLRNcEls川ULRTI0N〕

lERR0R・lF一口［一G1〕

図14水収支の誤差の累積
P：降雨量，D
Fig．14

Cumulation

P＝rainfanof：土中の水分量の初期値との差，G：地下水流出量
error

derived from1oss of water balance．
amount，D〃：decrement of amount water in soil

G：9round・waterfrominitialamount，mn−o冊amount・

0．16mmになる．これに対して，一般方式のERRORはオ＝5分で2．0mm，チ＝780分で23．7mmになるから，新方式の水収支の誤差が一般方式の1／100以下になる．総降雨量が72．8mmであるから，一般方式の水収支の誤差は降雨量の33％に及んでいる．なお，土中の初期含水量は290mmであった．
一般方式による数値計算例は，流出量とテンショソと全ポテソシャルが検証実験に良く一致しているにもかかわらず水収支が満たされていない．このため一般方ヰによる計算例と検証実験とを土層各部の含水量の変化について比較するとあまり良い一致がみられないことが予想される．本報告の数値計算例は，一般方式による計算結果のテソショソと地下水流出量が検証実験に一致するように物理定数を同定しこの物理定数を新方式に流用したため，新方式の計算結果は検証実験と良い一致をみない．しかし，今後，新方式により物理定数を同定すれば，より精度の良い物理定数と数値計算結果が得られると思われる．数値実験による不確定なパラメータの同定には数値実験の多角検討が必要である．また，このパラメータを用いた数値予測は適用範囲に細心の注意を要する．

7．結語

降雨浸透を土中の飽和一不飽和域にわたり解析する計算誤差の小さい数値解析法を提案した．提案された解析法の計算誤差より発生する水収支の誤差は，本報告の解析例では，従来の解析法の誤差の1／100以下であった．

大型水理実験の実験結果と数値解析結果との比較により，数値解析法の検証を行なった。

一164一

鉛直降雨浸透の有限要素法による一解法一大倉提案された数値解析法は検証実験の次の現象（i），ii），iii））を表現している．

i）浸澗面の通過に伴うテソションの急激な増加．

ii）浸潤面が到達している部分において，降雨強度の変動に伴うテソションの変動は表層ほど大きく，深くなるに従って減衰する．

iii）浸潤面が到達している部分では雨が止むとほとんど同時にテンショソが回復し始める．また，テンションの回復する速度は表層に近いほど早い．

ただし，数値解析結果の次の定量的な事項（iv），v），vi））が検証実験と良い一致をみな
レ・．

iV）浸潤面の降下速度，地下水流出開始時刻，地下水流出ピーク時刻．

V）地下水流出のピーク流出強度．
Vi）浸潤面通過後のテソションの値．

これらの定量的事項の不一致は従来の解析法により同定した物理定数の一部（K一θ，ψ一θの関係）を流用したためである．今後，これらの物理定数をより正確に同定すれば，提案された数値解析法で，iV），V），Vi）の定量的な現象もより精度よく解析できると思われる．

従来の方法による解析結果は水収支に大きな誤差（誤差は総降雨量の33％）があるにかかわらず，Vi）を除いて検証実験の結果によく一致している．このように致命的欠陥を持ちな

がら，ある一面で検証実験によく一致する場合があり，不確定なパラメータを有する現象の数f直実験はパラメータの同定と数値実験の適用範囲に細心の注意を払わなければならない．

検証実験は科学技術庁特別研究促進貿・地下水の水収支の解析手法に関する総合研究r水理模型実験による地下水の基本的特性に関する研究」のもとに行なわれた．
謝辞

検証実験において，テ1／ショソの測定に援助していただいた，国立防災科学技術セソター

第3研究部佐藤照子氏，同第4研究部御子柴正氏の両氏に感謝の意を表します．
参

考

文

献

1︶

赤井・大西・西垣（1977）：有限要素法による飽和一不飽和浸透流の解析，土木学会論文栂告集，

2︶第264号．87−96．
Finlayson，B．A．（1972）：The
AcademicPress，New

3︶

Irmay，S。（1954）：

Methodof

Weighted

Residuals

andVariatio施1Principles，Unsaturated

Soils，乃η舳．Agu，Vol．35，

York．
On

the

Hydrau1ic

Conductivityof

No．3，463−467．
4︶K1ute，A．（1952）：A
rated

Numerical

Methodfor

So1vi㎎theF1ow

Equation

for

Water

in

Unsatu−Materials，804Z8d．，73，105−116．

5︶Mua1em，A、（1974）1A

Conceptual

Mode1of

Hysteresis，Wα伽Rωo〃w∫地蜘〃一，Vol．10，

No．3，514＿520．
6︶

Mualem，A．（1977）：Extension

of

the

Similarity

Hypothesis

−165一

Used

for

Mode1ing

the

Soil

Water

国立防災科学技術セソター研究報告第22号

1979年10月

Characteristics，W励ぴル50〃o孤Rθ∫θαグ6ん，VoL13，No．4，773−780．

7）大倉博（1977）：不飽和浸透を考慮した降雨浸透の有限要素法による一解法一定常流一r国立防災科学技術セソター研究報告，第18号，51−72．
8）大島・毛利・中川（1976）：大気汚染制御のための汚染濃度予測の一方式，生産研究，28巻，3号，76−82．

9）Perrers，S．J．，Wats㎝，K．K．（1977）：Numerical
Redistribution，Wακr10）Richards

Analysisof

Two−Dimensional

InmtrationandRθ∫o〃κωR8∫θ〃6ん，Vol．13，No．4，781＿790．

L．A．（1931）：Capinary

ConductionofLiquids

throughporousMediums，Physics，Vol．

1，318−333．
11）Rubin，J。（1967）：Numericaltribution

of

Soil

MethodforAnalyzing

Hysteresis・Af正ected，Post・In丘1tration

Redis−

Moisture，8o〃8cづ．A〃．丹oc．，Vol．31，13−19．

12）富永雅樹（1978）：地下水の酒養推定のための水主里実験（1），国立防災科学技術セソター報告，第20号．123−134．

13）山村・久楽（1972）：堤防への浸透と堤体の安全性，土木研究所報告145号の3，41−71、

付録非線形連立方程式の近似解と真の解との誤差を近似的に求める方法．
係数ん伽が未知数乃j（プ＝1，2，…，刎，・・，〃，…）の単調連続関数となる準線型連立一次方程式
λ閉伽（んj）・乃㎜＝B勿（乃j）（A．1）

を考える（以下，特にことわらないかぎり，式中の同一項に同じ添字をもつ量が繰返し現われるときは，その添字について総和をとるものとする）．近似解を幻とし，舛と真の解幻との誤差をωとする．
〃j＝μ十幻（A．2）（A．1）式の幻に幻を代入し，伽を未知数とする連立一次方程式
A腕㎜（ゐ多）・㎜＝B閉（勿）（A．3）

を解くと，伽は勿の関数となる．
㎜＝

㎜（μ）

（A．4）

上式のμにゐjを代入すると，伽の定義より，〃㎜＝㎜（ゐj）＝肌（μ十ω）（A．5）となる．（A．5）を眺についてTaylar展開して眺の二次以上の項を無視すると∂㎜

ん肌＝柵柵）＝伽（乃5）十∂ん毒腕（A・6）

上式と（A．2）式より，誤差ωは連立方程式○刎ω＝ろ肌

．

∂
α卿＝

㎜

∂ω

6㎜＝〃㍍一

一δ㈹

（A．7）

（A8）

㎜（μ）

（A．9）

を解いて近似的に求められる．ここで，δ刎はクロネッカーのデルタをあらわす．βを適当な定数として，μ十βωを新たに舛として反復を繰返すと，ゐ多の伽への収束が予想される．
（1979年6月22日

一166一

原稿受理）