地下水流

11.1 地下水流の基本式 (Fundamental equations)

11.1.1 Darcy の法則

これは地下水流動に関する法則であり、1856年に Darcy によって図-11.1に示すような均質な砂柱を用いて実験的に得られた。飽和砂柱を通過する比流量 $q = Q/A$ ($Q$: 流量、$A$: 断面積) は、動水勾配 $I = \Delta h / l$ ($\Delta h$: 水頭差、$l$: 砂柱の長さ) のとき

q = K \frac{\Delta h}{l} \tag{11.1}

で与えられる。ここに、$K$ は透水係数と呼ばれ速度の次元をもつ。比流量 $q$ は速度の次元をもつのでダルシー流速あるいは見かけの流速とも呼ばれる。間隙を通過する水の平均流速 $u$ は、有効空隙率 $n_e$ のとき $v = q/n_e$ である。

透水係数 $K$ は水 (一般には通過する流体) の密度 $\rho$ と粘性係数 $\mu$、帯水層 (地下水流動層) の間隙の大きさや構造などにより定まるので $K = k \rho g / \mu = kg / \nu$ ($g$: 重力の加速度、$\nu$: 動粘性係数) ともかかれる。ここに、$k$ は固有浸透係数と呼ばれ面積の次元をもつ。

Darcy の法則を、透水係数が方向により異なる異方性帯水層における3次元流に拡張してテンソル表示すると

q_i = K_{ij} I_j = - K_{ij} \frac{\partial h}{\partial x_j}, \quad i, j=1,2,3 \tag{11.2}

である。ここに、添字の繰返しは和を示し、例えば上式より、$x_i$ 方向の比流量 $q_i$ は $x_j$ 方向の動水勾配 $I_j$ を用いて、$q_i = K_{i1} I_1 + K_{i2} I_2 + K_{i3} I_3$, ($i = 1,2,3$) で与えられ、また $K_{ij}$ は透水係数テンソルであり、実用的には $i \neq j$ のとき $K_{ij} = 0$ である。さらに $I_j = - \partial h/\partial x_j$: 動水勾配、$h = p / \rho g + x_3$: 全水頭 (ピエゾ水頭)、$p$: 間隙水圧、$x_3$: 鉛直座標である。

Darcy 則は理論的には Navier-Stokes の式から誘導される。Darcy 則は地下水流が層流のとき成り立ち、$q$ と $I$ は線形関係にある。レイノルズ数 $Re = qd/\nu$ ($d$: 一般に土砂粒子の代表径) が $1 \sim 10$ 程度をこえると地下水流は徐々に乱流に遷移し、$q$ は $\sqrt{I}$ に比例するようになる[164]。

11.1.2 連続の式

飽和地下水流の連続式は

S_s \frac{\partial h}{\partial t} = - \frac{\partial q_i}{\partial x_i}, \quad (i = 1,2,3) \tag{11.3}

である。ここに、$S_s$ は比貯留量 ($= \rho g (\alpha + n \beta)$)、$\alpha$ は帯水層の圧縮率、$\beta$ は水の圧縮率、$n$ は空隙率である。土砂層の圧縮率は $10^{-6} \sim 10^{-9} \text{m}^2 \cdot \text{N}^{-1}$ であり、水の圧縮率 $4.4 \times 10^{-10} \text{m}^2 \cdot \text{N}^{-1}$ よりかなり大きい。水および帯水層が非圧縮性とみなせるとき、例えば不圧地下水流、すなわち地下水面をもつ流れの場合、式(11.3) は

\frac{\partial q_x}{\partial x} + \frac{\partial q_y}{\partial y} + \frac{\partial q_z}{\partial z} = 0, \quad (i = 1,2) \tag{11.4}

となる。ここに、$q_x$, $q_y$, $q_z$ は $q$ の $x$, $y$, $z$ 成分である。

11.1.3 基礎方程式

基礎方程式は、Darcy の法則と連続の式 (式(11.2) と式(11.3)) より

S_s \frac{\partial h}{\partial t} = \frac{\partial}{\partial x_i} \left( K_{ij} \frac{\partial h}{\partial x_j} \right), \quad i, j=1,2,3 \tag{11.5}

となる。異方性帯水層で $i \neq j$ のとき $K_i \neq K_{jj}$ であれば

S_s \frac{\partial h}{\partial t} = \frac{\partial}{\partial x} \left( K_x \frac{\partial h}{\partial x} \right) + \frac{\partial}{\partial y} \left( K_y \frac{\partial h}{\partial y} \right) + \frac{\partial}{\partial z} \left( K_z \frac{\partial h}{\partial z} \right) \tag{11.6}

となる。ここに $K_x$, $K_y$, $K_z$ は $x$, $y$, $z$ 方向の透水係数である。等方・均質な非圧縮性帯水層で、水も非圧縮性とすれば、式(11.6) は

\nabla^2 h = \frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2} + \frac{\partial^2 h}{\partial z^2} = 0 \tag{11.7}

となり、Laplace の式になる。異方性帯水層の場合でも、(帯水層と水が非圧縮性のとき)、$y_1 = y \sqrt{K_y/K_x}$, $z_1 = z \sqrt{K_z/K_x}$ とおけば、$x$-$y_1$-$z_1$ 座標系では式(11.6)は等方性の場合と同じく Laplace の式になる。

地下水流の境界条件は普通次のようである。

① 給水域との境界: $h = H$ ($H$: 各境界での既知水頭)

② 流量規定の境界: $K_{ij} (\partial h/\partial x_j) n_i = -V$ ($V$: 既知の流量、$n_i$: 境界の単位法線ベクトルの $i$ 成分)

③ 自由地下水面: $h = z_f$ または $K_{ij} (\partial h/\partial x_j) n_i = N n_3$ ($z_f$: 基準面から測った水面の高さ、$N$: 水面を通る鉛直浸透量)

④ 浸出面: $h = z_s$ ($z_s$: 基準面から地下水のしみ出し点までの高さ)

11.2 飽和地下水流 (Saturated flow)

地下水の流れがほぼ水平で鉛直速度成分が無視できるとき、水平方向の速度成分は一つの鉛直線上で一様であるとする仮定を Dupuit-Forchheimer (デュプイ・フォルヒハイマー) の仮定といい、その流れを準一様流という。

11.2.1 堤体内の流れ

(1) 直立堤体内の定常流

図-11.2に示すように直立堤体あるいは地盤の上・下流2地点 (距離 $L$) における水位・水頭が既知のとき、準一様流を仮定すれば、単位幅当りの流量 $Q$ と水位・水頭は、不圧地下水流のとき

Q = K \frac{H_1^2 - H_2^2}{2L}, \quad \frac{H_1 - H_2}{L} = I \tag{11.8}

被圧地下水流のとき

Q = KB \frac{H_o - H_i}{L}, \quad \frac{H_o - H_i}{L} = I \tag{11.9}

である。ここに、$B$ は被圧帯水層の厚さである。

Polubarinova-Kochina [166] は、準一様流を仮定せず、厳密解を求め、浸出面の長さと流量を与えている。

(2) 直立堤体内の非定常流

貯水池の水位が突然 $H_o$ から $H$ だけ上昇する場合を考える (図-11.3)。

準一様流を仮定すると、基礎式は

n_e \frac{\partial h}{\partial t} = K \frac{\partial}{\partial x} \left[ (H_o + h) \frac{\partial h}{\partial x} \right] \tag{11.10}

となる。ここに、$n_e$ は有効空隙率である。いま変数変換

u = \frac{x}{\sqrt{K H_o t}}, \quad \varphi = \frac{h}{H} \tag{11.11}

を行うと、式(11.10)は非線形常微分方程式

\frac{d^2 u}{d \varphi^2} + \frac{2}{\varphi} \frac{du}{d \varphi} \frac{d \varphi}{du} + 4u \varphi = 0 \tag{11.12}

となる。$H/H_o$ のさまざまな値に対して $u$ と $\varphi$ の関係が Polubarinova-Kochina [165], [169] によって与えられている。

11.2.2 構造物の下を回る流れ

この流れの代表的なものは、矢板の下を回る浸透流と不透水性堤体の下部を通る流れであり、このような流れは楕円浸透流に分類される。例えば、矢板で河床などの透水性地盤を仕切るとき、矢板の上・下流の水位差によって生ずる浸透流の流線は矢板下端を共通の焦点とする楕円群で与えられる。また堤体の下部を回る浸透流の流線は堤体の上・下流端を共通の焦点とする楕円群で与えられる。

浸透流量および構造物に沿っての動水勾配 (あるいは浸透圧) は、上・下流の地盤高の差、堤体下部の遮水壁の有無、矢板の数と位置、透水層の構造、不透水性基盤の有無など、さまざまな条件によって異なる。浸透流量の算定に加えて、矢板の場合には浸透水と土砂の局所的な噴出や浸透水圧による地盤の膨れ上りの可能性を検討するために動水勾配および水圧分布を知ることが重要であり、堤体下部の流れでは堤体の安定性を検討するために揚圧力の分布を知ることが必要である。

諸条件下における矢板および堤体を回る浸透流の解は、例えば Muskat [166]、Harr [167]、Polubarinova-Kochina [164b] によって与えられている。

11.2.3 井戸の水理[168]

井戸理論は主として揚水・注水量と水位・水頭 (あるいは水位低下量) との関係を扱うが、ときには帯水層や岩盤の水理定数を原位置試験で決定する際の理論的根拠となる。以下では、帯水層は均質、等方性で無限に広がり、その不透水性基盤は水平で、水位・水頭は基盤から測るものとする。また井戸は完全貫入井戸、すなわち帯水層を貫いて基盤に達するものとする。

井戸水理の基礎式は、流れが準一様流であるとして、円筒座標系 ($r$, $\theta$, $z$) では

\lambda \frac{\partial h}{\partial t} = \frac{1}{r} \frac{\partial}{\partial r} \left( K b r \frac{\partial h}{\partial r} \right) \tag{11.13}

である。上式において、被圧帯水層の場合は $b = B$ (層厚) および $\lambda = S_s$ (貯留係数) であり、不圧帯水層の場合は $b = h$ (水位) および $\lambda = n_e$ (有効空隙率) である。

(1) 被圧帯水層からの定常揚水

式(11.13)を $r$ について2回積分し、境界条件 $r = r_w$ で $h = h_w$ および $r = R$ で $h = H$ を用いると、揚水量 $Q_w$ は

Q_w = \frac{2 \pi K B (H - h_w)}{\log_e (R/r_w)} = \frac{2 \pi T (s_1 - s_w)}{\log_e (r_1/r_w)} \tag{11.14}

で与えられる。ここに、$r_w$ は井戸の半径、$h_w$ は井戸内水位、$T$ は透水量係数 ($= KB$)、$s_1$ と $s_w$ はそれぞれ $r = r_1$ と $r = r_w$ での初期水頭 (水平) からの水頭低下量である。この式は Thiem (ティーム) の平衡式と呼ばれる。

(2) 不圧帯水層からの定常揚水

揚水量 $Q_w$ は、式(11.13) から式(11.14) を導くのとほぼ同様にして、

Q_w = \frac{\pi K (H^2 - h_w^2)}{\log_e (R/r_w)} \tag{11.15}

で与えられる。水位低下量 $s = H_o - h$ が原水位 $H_o$ と比べて十分小さいとき、上式は

Q_w = \frac{2 \pi T' (s_i - s_w)}{\log_e (r_i/r_w)} \tag{11.16}

となる。ここに、$s_i = s_i - s_w^2/2H_o$ および $T' = K H_o$ であり、$s_i$ は $r = r_i$ での水位低下量である。

無限帯水層における定常流は実際には存在しえない。しかし現実には揚水による水位低下圏は有限で、この水位低下圏の半径 $R$ を影響半径と呼ぶ。普通 $R/r_w = 3000 \sim 5000$ あるいは $R = 500 \sim 1000 \text{m}$ の値がとられる。

(3) 被圧帯水層からの非定常揚水

式(11.13)は水頭低下量 $s(r,t) = H_o - h(r,t)$ を用いて表わすと、$Kb = KB = T$ (透水量係数) および $\lambda = S_s$ (貯留係数) として、

S_s \frac{\partial s}{\partial t} = \frac{\partial^2 s}{\partial r^2} + \frac{1}{r} \frac{\partial s}{\partial r} \tag{11.17}

である。いま揚水量 $Q_w$ は一定、井戸半径 $r_w \to 0$ とし、揚水開始後の境界条件として

s(\infty, t) = 0 \quad \text{および} \quad \lim_{r \to 0} \left( 2 \pi r T \frac{\partial s}{\partial r} \right) = -Q_w \tag{11.18}

を設定すれば、式(11.17)の解は

s(r,t) = \frac{Q_w}{4 \pi T} \int_u^{\infty} \frac{e^{-x}}{x} dx = \frac{Q_w}{4 \pi T} W(u) \tag{11.19}

となる。ここに、$W(u)$ は井戸関数と呼ばれ

$W(u) = \int_u^{\infty} \frac{e^{-x}}{x} dx = -0.5772 - \ln u + u - \frac{u^2}{2 \cdot 2!} + \frac{u^3}{3 \cdot 3!} - \frac{u^4}{4 \cdot 4!} + \cdots$

u = \frac{S_s r^2}{4 T t} \tag{11.20}

であり、その数表[168] が与えられている。

(4) 不圧帯水層からの非定常揚水

原水位 $H_o$ に比べて揚水による水位低下量 $s = H_o - h$ が小さいとき、$u = n_e r^2/4 K H_o t$ として、水位低下量と揚水量 $Q_w$ との関係は

s = \frac{Q_w}{4 \pi K H_o} W(u) \tag{11.21}

で与えられる。

(5) 群井戸

接近して存在する複数の井戸で同時に揚水を行うと相互に水位・水頭低下の影響を受けるため、各井戸の揚水量は単独の井戸で揚水するより小さくなる。各井戸での水位・水頭低下量と揚水量との関係を重ね合せると

被圧帯水層のとき

H_o - h_i = \sum_{j=1}^M \frac{Q_j}{2 \pi T} \log_e \frac{R}{r_j} - \sum_{j=1}^{M'} \frac{Q_j}{2 \pi T} \log_e \frac{r_{ij}}{r_i} \tag{11.22}

不圧帯水層のとき

H_o^2 - h_i^2 = \sum_{j=1}^M \frac{Q_j}{\pi K} \log_e \frac{R}{r_j} - \sum_{j=1}^{M'} \frac{Q_j}{\pi K} \log_e \frac{r_{ij}}{r_i} \tag{11.23}

が得られる[169]。ここに、$i, j = 1,2,\ldots,M$ は井戸番号、$r_i$ は $i$ 番目の井戸の半径、$r_{ij} = r_{ji}$ は $i$ 番目と $j$ 番目の井戸の距離、$R$ は群井戸の影響半径、$\sum'$ は $i = j$ の場合を除いた和、である。上式で、各井戸の水位を与えると各々の揚水量が求まり、各揚水量を与えると各井戸の水位が求まる。また任意地点 $(x, y)$ の水位・水頭は

被圧帯水層のとき

h = H_o - \sum_{i=1}^M \frac{Q_i}{2 \pi T} \log_e \frac{R}{L_i} \tag{11.24}

不圧帯水層のとき

h^2 = H_o^2 - \sum_{i=1}^M \frac{Q_i}{\pi K} \log_e \frac{R}{L_i} \tag{11.25}

で算定できる。ここに、$L_i = \sqrt{(x - x_i)^2 + (y - y_i)^2}$、$(x_i, y_i)$ は $i$ 番目の井戸の中心座標である。

井戸理論では上述の条件以外に、部分貫入井戸 (井戸底が帯水層基盤に達しない井戸) の場合、水面からの浸透・蒸発や加圧層からの漏水がある場合、帯水層が傾斜層や累層であったり、その広がりが有限である場合など、さまざまな条件に対応して解や算定式が得られている。それらについては例えば文献5c)、169)、170)を参照のこと。

11.3 不飽和浸透流 (Unsaturated flow)

11.3.1 基礎方程式

地表面と地下水面にはさまれた部分は不飽和帯と呼ばれ、この部分の水は、吸着力によって土粒子表面に薄膜として、また毛管力によって土粒子間の間隙にメニスカスを形成して存在する。不飽和帯の水圧は大気圧よりも低く、圧力水頭 (圧力ポテンシャル) は水分量と密接に関係する。

Richards [171] は、Darcy の法則が不飽和帯の流れにも適用できるものと仮定して、次式を導いた。

\frac{\partial \theta}{\partial t} = \nabla [K(\psi) \nabla (\psi + z)] + \frac{\partial K(\psi)}{\partial z}, \quad (i = 1,2,3) \tag{11.26}

ここに、$\theta$ は体積含水率、$K(\psi)$ は不飽和透水係数、$\psi$ は圧力水頭、$z$ は鉛直座標軸 (上向き正)、である。$\theta$ と $\psi$ が与えられれば、比水分容量 $C(\psi) = d\theta/d\psi$ を導入することで $\partial \theta/\partial t = C(\psi) \partial \psi/\partial t$ となる。そこで $K(\psi)$ と $\psi$ との関係が既知のとき、式(11.26) を解いて不飽和帯の圧力水頭の分布が得られる。この式はまた飽和・不飽和浸透流を同時に取り扱えるので、堤体内の非定常浸透流の解析[172], [173] や雨水浸透施設の浸透解析[174], [175] にも適用される。

不飽和浸透流を気液2相流として扱う解析も行われる。このときの基礎方程式は、水と空気の流れのそれぞれに Darcy の法則が適用できるとして導かれる[176]。高木・森下[177] はこの方法によって不飽和鉛直浸透の問題を解析し、間隙空気が浸透に及ぼす効果について考察している。

11.3.2 水分特性曲線

体積含水率 $\theta$、透水係数 $K$ (あるいは相対透水度 $k_r$) および圧力水頭 $\psi$ の間の関係を示す水分特性曲線は図-11.4のようにヒステリシスをもつが、解析では多くの場合1価関数であると仮定される。例えば Brooks-Corey [178] は排水過程について、

$K = K_s = \left( \frac{\theta}{\theta_s} \right)^{2+\lambda}$

\frac{\psi}{\psi_a} = \left( \frac{\theta - \theta_r}{\theta_s - \theta_r} \right)^{-\frac{1}{\lambda}} \tag{11.27, 11.28}

を与えている。ここに、$\psi_a$ は空気侵入値と呼ばれ排水過程で空気が最初に土中に入るときの圧力水頭、$K_s$ は飽和時の透水係数、$n$ は空隙率、$\theta_r$ は残留水分量、$\lambda$ は定数、である。Campbell [179] もほぼ同様の関係を与えている。

11.4 地下水流における物質と熱の移動 (Mass transport and heat conduction in groundwater flow)

11.4.1 地下密度流と移流分散現象

ここでは臨海帯水層への塩水侵入を例として、地下密度流および地下水流による溶質の移動について述べる。

海水の比重が $1.02 \sim 1.03$ であるため、臨海帯水層では、海に向かって流れる淡水の下に海水がくさび状に侵入する。塩水くさびの挙動を解析する際に、淡水と塩水が混合しない (あるいは混合域がきわめて狭い) という前提による方法と、両者が混合するという前提に立つ方法がある。前者では淡塩水境界面の形状 (界面形状) および移動が、後者では塩分濃度の分布が主たる解析対象である。

(1) 地下2層流としての塩水くさび

淡水と海水が混合しないとき、図-11.5のような流れ場を考える。普通水平流速は鉛直流速より十分大きいので、流れは準一様流として扱える。いま淡塩水境界面の静的平衡および淡水内の静水圧分布を仮定すれば

h_s = h_f / \varepsilon \quad (\text{ただし } \varepsilon = \Delta \rho / \rho) \tag{11.29}

が得られる。ここに、$h_f$ は海面から測った地下水面の高さ、$h_s$ は海面から淡塩水境界面までの深さ、$\rho$ は淡水の密度、$\rho_s$ は海水の密度、$\Delta \rho$ は海水と淡水の密度差 ($= \rho_s - \rho$)、である。式(11.29) は Ghyben-Herzberg (ガイベン・ヘルツベルグ) の近似と呼ばれる[168b]。

淡水流と海水流の各々に Darcy の法則を適用し、式(11.29) と淡塩水境界面での水圧のつり合いを考慮すると、淡水と海水の比流量 ($q_f$ と $q_s$) は

淡水域:

q_f = - K \frac{\partial h_f}{\partial x} \tag{11.30}

海水域:

q_s = - K \sqrt{\varepsilon} \frac{\partial h_s}{\partial x} + K \varepsilon \frac{\partial h_f}{\partial x} \tag{11.31}

となる。流れの基礎方程式は、連続の式と上記の2式から

淡水域:

\frac{\partial h_f}{\partial x} + \frac{\partial}{\partial x} \left[ K (h_f + h_s) \frac{\partial h_f}{\partial x} \right] = 0 \tag{11.32}

海水域:

\frac{\partial h_s}{\partial t} + \frac{\partial}{\partial x} \left[ K (B_s - h_s) \left( \frac{\partial h_s}{\partial x} - \sqrt{\varepsilon} \frac{\partial h_f}{\partial x} \right) \right] = 0 \tag{11.33}

である。

a. 定常淡塩水境界面流れが定常であれば塩水くさびは静止しており $q_s = 0$ である。不圧帯水層内の界面形状は、淡水の単位幅流量 $Q_L (= \text{const.}) = - K (h_f + h_s) (\partial h_f/\partial x)$ なる関係に式(11.29) を代入し、$\varepsilon << 1$ を考慮すれば

\frac{\varepsilon K}{Q_L} h_s^2 = \left( \frac{2 \varepsilon K x}{Q_L} \right)^{1/2} \tag{11.34}

で表わされる。この式は、$x = 0$ で $h_s = 0$ を仮定して導かれており、したがって淡水の海への流出速度が無限大であるという矛盾をもっている。海岸線の近くでは淡水流速はもちろん有限であり、また鉛直流速成分も無視できない。玉井・嶋[180] はこの点を考慮して、界面形状が

\frac{h_s}{h_0} = \left[ \left( \frac{2 \varepsilon K x}{Q_L h_0} + 0.5 \right)^2 - 0.25 \right]^{1/2} \tag{11.35}

で与えられることを示した。

被圧帯水層内の界面形状についても式(11.34)や式(11.35)に相当する関係式が導かれている[169]。佐藤・渡辺[181] は透水係数の異なる水平な2層間内の界面形状を求めている。

浸透領域・境界条件が比較的単純な場合、淡塩水境界面は Zhukovsky (ジューコフスキー) 関数

$\Omega = \varphi + i \psi = - i e^{- \zeta}, \quad \zeta = \xi + i \eta, \quad z = x + i y$

\varphi = \phi - \varepsilon K y, \quad \psi = \phi - \varepsilon K x, \quad i = \sqrt{-1} \tag{11.36}

を用いて解析できる。ここに、$\zeta = \phi + i \psi$、$z = x + i y$ であり、$x$ と $y$ は物理面内の水平および鉛直座標、$\phi$ は速度ポテンシャル ($= K (y + p/\rho g)$)、$\psi$ は流れ関数である。上田・崎山[182] は注水堰による河口貯水池の塩水遮断、崎山[183] は河口貯水池地盤からの淡水取水時の塩水混入限界、上田・杉尾[184] は円形暗渠からの淡水注入による塩水侵入防止について解析している。

境界条件が複雑でも浸出面がなければ、Cauchy-Riemann (コーシー・リーマン) の条件と逆 Laplace の式を適用し、数値解析を援用して淡塩水境界面を定めることができる。この手法により崎山ら[185] は河口貯水池における塩水排除について解析している。

塩水侵入の注水井群による制御問題では、上田ら[186] は定常地下密度流の流量ポテンシャルを導き、有限要素法と線形計画法を用いて最適注揚水状態を求めている。

b. 非定常淡塩水境界面

塩水くさびの侵入過程で、はじめ上に凸であった淡塩水境界面は進行とともに下に凸になる。嶋[53c], [187] は海岸線部での淡水層厚さの時間的変化および流量分布の2次近似をとり入れ、上記の侵入過程をほぼ再現できる解析結果を得ている。

複雑な境界条件のもとで注水や取水に伴う塩水くさびの挙動を追跡するには数値解析が有力な手段となる。川谷[188] は注水による塩水くさびの変形を差分法で、吉野・上田[189] は有限要素法によって海面変動に伴うくさびの挙動も含めて解析し、佐藤ら[190] は、境界積分方程式法を用いて、地下空洞掘削によって誘発される塩水侵入の問題を扱っている。

(2) 移流分散現象としての塩水侵入

淡水と塩水が混合するという前提では、海水侵入を地下水流における塩分の移流分散現象として解析する。すなわち塩水くさびの伸縮・変形を塩分濃度分布の変化としてとらえ、その解析にあたっては「浸透流」と「移流分散」の基礎方程式を連立させて解く。

飽和浸透流における移流分散の基礎方程式は、吸着や減衰が無視できるとき

n \frac{\partial C}{\partial t} = \frac{\partial}{\partial x_i} \left( n D_{ij} \frac{\partial C}{\partial x_j} \right) - u_i \frac{\partial C}{\partial x_i}, \quad i, j=1,2,3 \tag{11.37}

である。ここに、$C$ は溶質の濃度、$D_{ij}$ は分散係数、$n$ は空隙率、$q_i$ はダルシー流速である。

この立場の解析は、主として塩水くさびの非定常な挙動や地質構造の複雑な帯水層における形状を把握するために利用され、数値解析によって実行される。Segol-Pinder [191] は塩水と淡水の密度勾配による流れも考慮して有限要素解析を行い、Frind [192] や川谷[193] は解法の簡易化と適用範囲の拡張を行っている。また河野ら[194] は不飽和浸透流における塩水化現象を取り扱っている。

11.4.2 地下水流における伝熱

飽和地盤・岩盤における伝熱現象の解明は工学的な見地からも、地中埋設物からの熱損失、燃料の地下備蓄や放射性廃棄物の地下処分による高温伝熱、液体ガスの貯蔵や地盤凍結工法に伴う凍結・融解などに関して重要である。

地盤は固体、液体、気体から構成された複雑な系であり、一般に熱移動も伝導、対流、輻射がすべて起っている。しかし実際には伝導が最も大きな比重を占め、対流がこれに次ぎ、輻射は無視できるほど小さい。

伝熱現象の解析は「地下水流」と「伝熱」の基礎方程式を連立させて行う。飽和浸透流の場における伝熱の基礎方程式は

(pc) \frac{\partial T}{\partial t} + (pc)_w q_i \frac{\partial T}{\partial x_i} - \frac{\partial}{\partial x_i} \left( K_{ij} \frac{\partial T}{\partial x_j} \right) = 0 \tag{11.38}

である。ここに、$T$ は温度、$(pc)$ は飽和土の等価体積熱容量、$(pc)_w$ は水の体積熱容量、$K_{ij}$ は熱伝導率 (テンソル)、$q_i$ はダルシー流速、である。

飽和土における熱伝導は、① 固体 (土粒子) と液体 (水) が一体となって熱を伝える、② 固体のみが熱を伝える、③ 液体のみが熱を伝える、という3ルートで起る。それゆえ、飽和土の熱伝導率は土と水の構成比を考慮して算定される[195]。また土の等価体積熱容量 $(pc)$ は単位体積中の固体、液体、気体の熱容量の和で与えられる。

また水は相変化するとき潜熱を発生・吸収するので、凍結・融解にかかわる問題を取り扱うには相変化による物性値の変化とともに潜熱も考慮しなければならない。

流れの場および境界条件が複雑な場合、地下水流における伝熱現象の解析には数値解析が用いられる。佐藤・伊藤[196] は地下空洞周辺の岩盤および風化帯における伝熱現象を解析しており、このとき対流の卓越する場の解析にも適用性の高い変形 FLIC 法を用いている。地下水流における凍結・融解の問題を、川谷・渡部[197] は潜熱を含む見かけの体積熱容量を導入して解析している。