人工湖補注地下水之力學分析

人工湖補注地下水之力學分析

蔡佩勳¹ 林建翰 ¹ 邱勢棠¹

1朝陽科技大學營建工程系

摘要

本研究使用FLAC 5.0模擬人工湖湖水下滲地下之過程，改變土壤初始飽和度及討論標準 砂、壤土、粉質壤土等三種土壤對人工湖附近地表高程變化及對地下水補注量進行探討。由 本研究之結果顯示，以較高的初始飽和度使地表隆起量越大，初始飽和度對最終地下水補注 量之影響不大。砂土飽和滲透係數高，其最終地下水補注量越多，但地表最終的隆起量卻較 小。

關鍵字：非飽和土壤、地下水補注、人工湖。

一 、前言

鑿井抽取地下水之使用率逐年提高，然而地下水超抽可能發生地下水位降低、水質惡化、

海水入侵、地下水鹽化及地層下陷等問題。地下水之主要來源來自地表水的補注，當地下水 的抽取量超過地表水的補注量時，則會造成蓄水層中儲水量減少，地層中孔隙水壓的降低，

因而發生地層下陷。

James et al. (1987)[1]提出人工湖地點的選擇會明顯地影響到人工湖之地表水滯留量及補 注能力。倪春發等人(1999)[2]以 2D-FEMFAT 來模擬砂土、粉土、黏土三種均勻土壤之定水頭 補注模擬，其結果顯示在地下水補注之影響區域以及補注量的多寡主要決定於土層飽和的滲 透係數。曾建璋等人(1999)[3]針對嘉義東石鄉鱉鼓地區進行人工湖補注砂槽實驗，其結果顯 示地下水位的升降會影響補注水入滲的進行。當地下水位越高時，則湖水入滲越慢，當地下 水位等於或高於人工湖湖底時，入滲作用即停止。周湘俊等人(2000)[4]以 2D-FEMFAT 模擬人 工湖補注地下水，模擬結果顯示人工湖在同一水深的情況下，寬度越大補注效益越好；而在 同一寬度的情況下，湖水越深補注效益越佳。劉建榮(2001)[5]之研究指出水分特徵曲線參數 之 VG 模式α 與 n 參數對砂箱模型實驗之暫態補注量有極大的影響，且參數α對入滲量之影響 較參數 n 大。翁琮哲(2001)[6]利用 MODFLOW 三維地下水流模式進行不同人工湖補注方案，

由其結果顯示，當地表坡度越陡和地下水位越高時，其補注的水會越早流出地表，而無法完 全補注到下游部分。劉致翔(2011)[7]探討地下水補注對防止海水入侵拘限含水層的影響，當 補注量提高時，能有效阻擋海水往內陸推進，較高的水力傳導係數會使海水入侵現象更為加 劇。

過去有關人工湖補注的文獻中，大多針對人工湖附近的地質條件、人工湖幾何尺寸、淤 泥堵塞湖底的情況、人工湖水位高低、及阻水層位置對地下水補注量之影響分析。對於湖水 下滲至地下時，湖水在滲流過程中間接影響土壤內部孔隙水壓，因此將會影響到地表隆起，

這方面的力學分析以室內試驗居多，參數影響分析之研究仍不多，對於欲以人工湖補注地下 水且兼具提升地表高程的方案中，人工湖設置地點之選擇有賴各影響參數對補注量與地表提 升量之影響的了解，因此本研究以數值方法進行參數影響分析，探討各影響參數對最終地下

水補注量及提升地表高程之影響。

二 、 數值分析方法

2.1 土壤水分特徵曲線

基 質 吸 力 與 土 壤 含 水 量 或 飽 和 度 之 關 係 曲 線 ， 稱 為 土 壤 水 分 特 徵 曲 線 (soil- water characteristic curve)，Harrison and Blight (2000)[8]指出水分特徵曲線為非飽和土壤之基質吸力 (metric suction)與飽和度的關係曲線，受土壤顆粒成分、土壤結構與孔隙大小、型態等因素之 影響，土壤對水分的保水能力會隨著基質吸力之大小而改變。

2.2 二相流(Two-Phase Flow)模式

van Genuchten (1980)[9]根據 Mualem 公式的土壤飽和度與張力水頭(

_Ψ

)之關係式進行積 分，得到水分特徵曲線關係式，此關係式中包含了三個擬合參數，分別為α、n、m，其中α 參數與空氣進氣值有關，n 參數與土壤孔隙分布有關，m 參數與曲線之不對稱相關。van Genuchten 所提出的方程式，如式 1 及式 2 所示。

[ ( )

ⁿ

]

^m

e

1 S 1

Ψ α +

=

(1)

( )

^wr

[ ( )

ⁿ

]

^m

w r w

1 S 1 1 S S

Ψ α +

− +

=

(2)

上式中，S_w=相對於 Se下所對應之飽和度；S =殘餘飽和度；S^w_r e=有效飽和度；

Ψ

=張力水頭。

FLAC(Itasca, 2005)[10] 之二相流模式(two phase flow model)對非飽和土壤中水與氣體之 流量方程式，可以達西定律表示如下：

k)

�x wg ρPw ( xj

∂ w ∂ w r kij

� w =

qi k (3)

k )

� x g g ρ

P g ( x j

∂ g ∂ g r μ w w μ k ij

� i =

q

^g

k

(4)

式中，q 為流量、 w

kij 為飽和滲透係數、 rκ 為相對滲透率、相對飽和度 S^e、µ為動力黏度、

P 為孔隙壓力、ρ 為密度。上標的 g 表示氣體、w 表示水。相對飽和度 Se與相對滲透率 rκ 之 關係公式如下：

]2 )a /a e1

�S

�(1 b[1 Se w =

κr (5)

a ] 2 /a e 1

� S c [1 e )

� S (1 g =

κ r

(6)

式中，a、b 和 c 為經驗常數。毛細壓力 Pw與有效飽和度 Se之關係公式如下：

� a 1) 1 /a �

� 1 (S e P o

= w ) (S

P

_w (7)

上式中，P0為土壤材料之相關參數。飽和度與孔隙壓力對土壤孔隙體積應變之組成律關係公 式如下所示：

t S e t q t

S t P K

n S _w

iw w

w w

w ∂

+ ∂

∂

−∂

∂ = +∂

∂

∂ (8)

t S e t q t

S t P K

n S _g

g g i

g g g

∂ + ∂

∂

−∂

∂ = +∂

∂

∂ (9)

上式中，n 為土壤孔隙率，e 為體積應變，K 與_w K 分別為水與氣體之體積模數，_g S 為氣體_g 之飽和度，它與S 之和等於 1。 _w

2.3 人工湖數值模型

因對稱之關係，人工湖之數值模型採用人工湖左半邊進行模擬分析，且以平面應變進行 分析。元素網格之長度為80m、深度為20m、人工湖底寬為4m、人工湖深度為3m，人工湖邊 坡坡角為27°，地下水位於人工湖湖底往下5m的位置，元素網格共切割成3600個元素，每個 元素長度為0.5m，如圖1所示。為了防止土層加載時，不致於發生整體位移，因此假設模型左 右邊界節點之水平方向位移為零，土層底部邊界各節點為水平及垂直方向位移為零。於本分 析中，滲流分析之邊界條件分為三種，(1)非浸潤邊界、(2)浸潤邊界及(3)入滲邊界。非浸潤邊 界主要假定該邊界的土壤為非飽和；浸潤邊界主要是假定此邊界的土壤為飽和狀態(例如在地 下水位以下或是人工湖表面浸水時)；入滲邊界主要假定當此邊界上的節點其孔隙水壓力為負 值時，土壤之水分能從該邊界流入或流出。由於地下水位於人工湖底下5m的位置，所以在地 下水位左右邊界，設定為浸潤邊界，地下水位以上的左右邊界設為入滲邊界，地表面設為非 浸潤邊界，土層底部設為不透水邊界。由於人工湖蓄水時，人工湖表面係浸在湖水中，所以 人工湖表面節點的邊界條件隨著湖水的上升更改為浸潤邊界，如圖1所示。

(a) 人工湖元素網格

Water Table 水平方向位移為零

(fix x)

水平及垂直方向位移為零 (fix x y)

(b) 力學邊界條件

Water Table 非浸潤邊界(fix nwpp) 入滲邊界(fix seepage)

浸潤邊界(fix pp sat)

不透水邊界

非浸潤邊界(fix nwpp) 浸水後改 浸潤邊界(fix pp sat)

(c) 滲流邊界條件 圖 1.人工湖邊界條件示意圖

2.4 力學-滲流耦合分析

標準砂、壤土與粉質壤土之力學組合律採用莫爾庫倫模式，本研究對標準砂進行三軸排 水試驗得到其力學參數(乾密度、楊氏係數與摩擦角等)，而非飽和土壤滲流分析採用二相流模 式，其相關參數將引用劉建榮(2001)[5]時砂箱滲流試驗反算所得到的水分特徵曲線參數及飽 和滲透係數。除了標準砂外，本研究所用之壤土與粉質壤土之力學參數，以三角座標分類法 來推出壤土與粉質壤土的各砂土、黏土、粉土百分比，再配合Girsang(2001)[11]實驗所得砂土、

黏土、粉土的基本材料參數，以這三種土壤百分比作為基準值，得出壤土與粉質壤土的材料 參數，標準砂、壤土與粉質壤土的力學材料參數，如表1所示。壤土與粉質壤土的二相流參數 將引用van Genuchten(1980)[9]試驗所得的水分特徵曲線參數，如表2所示，標準砂、壤土、粉 質壤土的水分特徵曲線圖，如圖2所示。

表 1 基本力學材料參數 土壤乾

密度 ρ(kg/m³)

體積模數 K(MPa)

剪力模數 G(MPa)

凝聚力 C(Pa)

摩擦 角

φ

土壤種類 1800 16.6 7.67 1 38 標準砂 1826 11.9 3.34 5200 23 壤土 1850 12.7 3.91 3100 26 粉質壤土

表 2 二相流模式之參數

土壤

經驗常數

土壤材料相 關之參數

Po(Pa)

殘餘 飽和 度S _r^w

(%)

孔隙 率

n (%)

飽和滲透係 數 k_s(m/s) a b c

標準砂 0.655 0.5 0.5 4459.1 0.5 0.47 4.8 粉質壤土 0.86 0.5 0.5 9810 0.5 0.47 3.5

壤土 0.51 0.5 0.5 8530 0.5 0.52 3.7

0 -20000 -40000 -60000 -80000 Capillary Pressure (Pa)

20 40 60 80 100

Saturation (%) ^Sand_Loam

Silt Loam

圖2.三種土壤之水分特徵曲線圖

因為每年雨季約為三個月(90 天)，故一年中人工湖假設有 90 天為滿水位，屬於地下水補 注時期，其他月份人工湖處於無水情況，即未進行地下水補注，所以本研究分析之總時間為 360 天，其中 90 天為地下水補注時期(滿水位的情況下)，剩下的天數停止補注時期(水位洩降 至無水位的情況下)。數值模型分為人工湖範圍(0~10m)與鄰近地表範圍(10m~80m)進行討論，

如圖 1 所示，探討人工湖補注及停止補注過程中，地下水位變化對上述兩區域的地表高程之 影響，並得到人工湖最終補注量，討論不同土壤初始飽和度及不同類型土壤對地表隆起量及 補注量之影響。

為了證明人工湖之地表隆起量係受到湖水入滲過程所引起，並非湖水重量所導致，本研 究分別對人工湖進行力學分析以及力學-滲流耦合分析，前者並不進行湖水下滲之滲流分析。

圖3為力學-滲流耦合分析與力學分析在90天與360天時之結果比較圖，由結果顯示，在地下水 補注90天後，如果僅採用力學分析，地表最大隆起量發生的位置在30m的地方(離人工湖左側 邊緣20m處)，隆起量約為0.03cm；如果採用力學-滲流耦合分析下，最大隆起量發生的位置在 14m的地方(離人工湖左側邊緣4m處)，隆起量約為2.8cm。在地下水補注90天後人工湖停止地 下水補注一直到第360天，如果僅採用力學分析，地表最大隆起量發生的位置在30m的地方，

隆起量約為0.002cm；如果採用力學-滲流耦合分析下，最大隆起量發生的位置在8m的地方，

隆起量約為

1 . 3 × 10

⁻²

cm

。

80 60 40 20 0

Distance (m) -0.02

-0.01 0 0.01 0.02 0.03

Y-displacement (m)

Only Mechanical Flow-Mechanical Coupled

80 60 40 20 0

Distance (m) -0.004

-0.003 -0.002 -0.001 0 0.001

Y-displacement (m)

Only Mechanical Flow-Mechanical Coupled

(a) Day 90 (b)Day 360 圖 3.力學-滲流耦合與純力學分析之比較圖

隆起量明顯地比90天小很多，此係因為左側邊界為可滲出之邊界條件，所以當停止補注 地下水時，土層之地下水將因滲流出邊界，使孔隙水壓降低，導致地表隆起量回降，如果有 阻水牆存在於右邊界，則右邊界將為不排水邊界，可以預期此土層之孔隙水壓應不會因停止

補注地下水而降低，地表隆起量也無回降問題。當僅採用力學分析而未加入滲流運算，因只 考慮湖水重量之影響，對地表之隆起量並無明顯的變化。當採用力學-滲流耦合分析，將湖水 入滲效應也考慮的話，其對地表隆起量有顯著的增加，隆起量增加的原因可能是在濕潤面的 推進下，土壤基質吸力的消散，使得土壤飽和度隨之增加，間接導致有效應力的降低。因此，

力學-滲流耦合分析可應用於人工湖補注地下水之數值模擬，而後述之參數影響分析，都會以 力學-滲流耦合的方式來進行。

三、分析結果與討論

3.1 土壤初始飽和度之影響

圖 4 為 90 天與 360 天之不同飽和度下地表各點垂直位移分布圖，由結果顯示，在第 90 天時(圖 4(a))，飽和度 50%的最大隆起量約為 2.8cm(發生位置於 14m 地方)、飽和度 60%的最 大隆起量約為 3.0cm(發生位置於 14m 地方)、飽和度 70%的最大隆起量約為 3.1cm(發生位置 於 14m 地方)，在補注期間，飽和度越大其地表隆起量就越多，但在停止補注(湖水開始洩降 至無水位)後，第 360 天時(圖 4(b))，飽和度 50%的最大隆起量約為

1 . 3 × 10

⁻²

cm

(發生位置於 8m 地方)、飽和度 60%的最大隆起量約為

− 5 . 4 × 10

⁻²

cm

(發生位置於 8m 地方)、飽和度 70%

的最大隆起量約為

− 1 . 2 × 10

⁻¹

cm

(發生位置於 8m 地方)。圖 5 為不同飽和度下補注量之歷時分 布圖，由圖 5 顯示，土壤在飽和度 50%、飽和度 60%跟飽和度 70%時最終補注量幾乎相同，

主要原因可能是補注量多寡跟土壤的飽和滲透係數有關，即使改變了土壤初始飽和度，但因 擁有相同的土壤飽和滲透係數，所以補注速率其差異不大，對應的補注量也差不多。

80 60 40 20 0

Distance (m) 0

0.01 0.02 0.03 0.04

Y-displacement (m) ^S=50%

S=60%

S=70%

80 60 40 20 0

Distance (m) -0.005

-0.004 -0.003 -0.002 -0.001 0 0.001

Y-displacement (m)

S=50%

S=60%

S=70%

(a) Day 90 (b) Day 360 圖 4.不同初始飽和度之地表垂直比較圖

0 100 200 300 400

Time (day) 0

10000000 20000000 30000000

Infiltration Volume(m3)

S=50%

S=60%

S=70%

圖 5.不同初始飽和度之補注量比較圖

3.2 土壤種類之影響

由分析結果顯示，粉質壤土在第 90 天時(圖 6(a))，人工湖之地表隆起量已大於砂土與壤 土的隆起量，砂土的最大隆起量約為 2.8cm(發生位置於 14m 地方)、粉質壤土的最大隆起量 約 3.1cm(發生位置於 14m 地方)、壤土的最大隆起量約 0.48cm(發生位置於 26m 地方)，這可 能是因為壤土飽和滲透係數較其他兩種土壤小之緣故，使孔隙水壓增加的速率較慢，因此隆 起的變量相較於其他兩種土壤少。第 360 天(圖 6(b))時，地表隆起量已趨穩定，最大地表隆起 量由大到小分別為粉質壤土約 3.1cm(發生位置於 6m 地方)、壤土約 2.7cm(發生位置於 4m 地 方)、砂土約

1 . 3 × 10

⁻²

cm

(發生位置於 8m 地方)。因砂土飽和滲透係數較壤土與粉質壤土者 高，使砂土在停止地下水補注時，隆起量回降趨勢明顯。三種土壤對地下水補注量之影響可 從圖 7 得出，砂土(飽和滲透係數

k

_s

= 4 . 8 × 10

⁻⁴

m / s

)最終補注量約為

2 . 5 × 10

⁷

m

³、粉質壤土(飽 和 滲 透 係 數

k

_s

= 3 . 5 × 10

⁻⁵

m / s

) 最 終 補 注 量 約 為

4 . 5 × 10

⁶

m

³ 、 壤 土 ( 飽 和 滲 透 係 數

s / m 10 7 . 3

k

_s

= ×

⁻⁶ )最終補注量約為

1 . 1 × 10

⁶

m

³，上述可知土壤飽和滲透係數越大，最終補注 量越大。

80 60 40 20 0

Distance (m) -0.02

-0.01 0 0.01 0.02 0.03 0.04

Y-displacement (m)

Sand Loam Silt Loam

80 60 40 20 0

Distance (m) -0.01

0 0.01 0.02 0.03 0.04

Y-displacement (m) ^Sand

Loam Silt Loam

(a) Day 90 (b) Day 360 圖 6.不同類型土壤之地表垂直位移比較圖

0 100 200 300 400

Time (day) 0

10000000 20000000 30000000

Infiltration Volume(m3)

Sand Loam Silt Loam

圖 7.不同類型土壤之補注量比較圖

四、 結論

1. 純粹只考慮人工湖湖水重量對鄰近地表變形之影響，而未加入湖水入滲土壤之分析 時，對鄰近地表之隆起量並無明顯的變化。但如果將湖水入滲也考慮的話(即力學- 滲流耦合分析)，對鄰近地表之隆起量有顯著的增加。地表隆起量增加的可能原因 是湖水下滲使土壤基質吸力的消散，土壤飽和度增加，導致有效應力的降低，因此 非飽和土壤之力學-滲流耦合分析可應用於人工湖補注地下水之數值模擬。

2. 不同初始飽和度(70%、60%、50%)中，以初始飽和度 70%之隆起量為最大，飽和度 60%次之，飽和度 50%最小，但初始飽和度對最終地下水補注量之影響不大。

3. 比較砂、壤土、粉質壤土對隆起量與補注量之影響，發現土壤飽和滲透係數越高，

其最終地下水補注量越多，但地表最終隆起量卻較小。所以如果以地下水之補注量 為考量時，人工湖應設置於透水性良好之砂土層，但其所回升的地表高程可能就不 如細粒料含量較多的土壤，如果要維持其地表隆起量，可利用阻水牆，一方面可避 免補注的地下水滲流至其他區域，另一方面也可維持地表隆起量以回升沉陷之土 層。

五、誌謝

本研究承蒙科技部計畫(NSC102-2221-E-324-025-) 之經費支持，使本研究能順利進行，

特此致謝。

參考文獻

1. James, W. P., Bell, J. F. and Leslie, D. L. (1987), “Size and location of detention storage,” J. of Water Resources Planning and Managent, ASCE, 113(1), 15-17.

2. 倪春發、許少華(1999)，人工池補注地下水水庫之效益初步模擬與影響因子，第三屆地下 水資源及水質保護研討會論文集。

3. 曾建璋、簡錦樹、丁澈士、吳銘志(1999)，人工湖與地下水之相互效應研究：實驗調查，

第十屆水利工程研討會。

4. 周湘俊、陳運生、倪春發、許少華(2000)，人工湖補注池之補注量模擬與分析，第十一屆 水利工程研討會。

5. 劉建榮(2001)，van Genuchten 土壤特性曲線參數對濕鋒模擬與暫態補注量之影響，私立 逢甲大學土木及水利工程研究所碩士論文。

6. 翁琮哲(2001)，地下水補注功效之探討，私立中原大學土木工程學系碩士論文。

7. 劉致翔(2011)，地下水補注對沿海含水層海水入侵整治之數值分析，國立成功大學水利及 海洋工程學系碩士論文。

8. Harrison, B. A. and Blight , G. E. (2000), “The use of indicator tests to estimate the drying leg of the soil- water characteristic curve,” Asian Conference on Unsaturated Soil, Balkma, Rotterdam.

9. van Genuchten, M. T. (1980), “A closed form equation for predicting the hydraulic conductivity of unsaturated soil,” Soil Sci. Soc. Am. J., 44, 892-898.

10. Itasca Consulting Group, Inc. (2005), FLAC-fast lagrangian analysis of continua, Ver. 5.0 user’s manual. Minneapolis, Minnesota, USA.

11. Girsang, C. H. (2001), “A numerical investigation of the seismic response of the aggregate pier

foundation system,” Master Thesis of the University of Virginia, Blacksburg.

ABSTRACT

This study was used FLAC 5.0 to simulate underground water recharge of an artificial lake. A parametric study is performed to understand the influence on the final amount of recharge water and ground surface uplift. The influent parameters include soil type and initial degree of saturation. The amount of surface uplift increases with increasing initial degree of saturation of soil during recharging phase. The initial degree of saturation is an insignificant factor on the final amount of recharge water. However, the final amount of recharging water in sand layer is larger than other types of soil due to its high permeability, but the amount of ground surface uplift is less in sand layer.

Keywords: unsaturated soil, groundwater recharge, artificial lake.