國立交通大學
土木工程學系
博士論文
河岸退縮數值計算模式之發展與應用
The development and application of numerical model of
riverbank retreat
研 究 生
：
姜世偉
指導教授
：
楊錦釧
博士
蔡東霖
博士
河岸退縮數值計算模式之發展與應用
The development and application of numerical model of
riverbank retreat
研 究 生 ： 姜世偉 Student ： Shih-Wei Chiang
指導教 授 ： 楊錦釧 Advisor ： Jinn-Chuang Yang
蔡東霖 Tung-Lin Tsai
國 立 交 通 大 學 土木工程學系
博 士 論 文
A Dissertation
Submitted to Department of Civil Engineering College of Engineering
National Chiao Tung University in Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy
in
Civil Engineering June 2011
Hsinchu, Taiwan, Republic of China
誌謝
承蒙楊教授錦釧與蔡教授東霖的悉心指導使論文得以完成，在此致上由衷的感 謝，並且感謝口詴委員黃教授良雄、詹教授錢登以及潘教授以文悉心給予指正與建議， 使本論文更臻完善。此外，亦感謝系上老師在課業上的啟發與教誨，求學期間使學生 受益匪淺。 感謝所有支持我的朋友以及多位學長、同學與學弟妹們給予的幫助與關懷。 最後對於家人的支持與鼓勵，感激之情，永銘於心。摘要
本研究之目的旨在發展並應用河岸退縮預測數值計算模式，模式中採用未飽和地 下水流理論，利用數值計算獲致河岸土體內部孔隙水壓隨時間變化之分布，藉以建構 河岸邊坡塊體破壞機制。此方法不但可以改進以往過去相關研究須仰賴地下水面線且 須以靜水壓分布假設始能求得沿破壞面上之孔隙水壓，並可利用邊界條件的處理技巧 將河川水位變動以及降雨因素之影響一併納入計算。此外，更進一步整合懸臂型破壞 與水流沖蝕機制，完整地考量河岸退縮的主要關鍵因素，據以探究其互制行為與演變 過程。 本研究首先以所建立的模式分別針對塊體破壞、懸臂型破壞與水流沖蝕三種破壞 機制進行探討，利用案例測詴結果歸納出河川水位升降決定其所提供之靜水壓力，以 及土壤滲透性與降雨特性交互作用下造成土體孔隙水壓變化之差異，對塊體破壞以及 懸臂型破壞二種機制均具決定性之影響，而水流沖蝕則主要與水力條件(底床坡降、河 川水位、水位歷線型態)以及土壤抗沖蝕能力(臨界剪應力、沖蝕係數)相關。接續整合 三種破壞機制，並由模擬結果得知河岸退縮是塊體破壞、懸臂型破壞與水流沖蝕三種 破壞機制反覆循環發生的過程，而在水流沖蝕較大的條件下，河岸退縮主要受水流沖 蝕的程度與範圍所影響。最後於濁水溪河段之案例應用結果顯示，此整合模式可提出 不同機制下的泥砂產量，並且對於河岸退縮長度可獲一定程度之預測。 關鍵詞：河岸退縮、未飽和地下水流、極限平衡法、邊坡穩定、塊體破壞、懸臂型破 壞、水流沖蝕Abstract
In this study, a numerical model of riverbank retreat has been developed and applied to a practical case. The governing equation of unsaturated groundwater flow is solved by implementing numerical method to obtain the transient distribution of the pore water pressure to evaluate riverbank stability with respect to mass failure. However, previous studies to compute the pore water pressure were usually based on groundwater table with hydrostatic pressure distribution hypothesis. The approach proposed in this study not only improves this shortcoming but also takes the effects of river stage variations and rainfall into account by defining boundary conditions. In addition, cantilever failure and fluvial erosion are incorporated into the model in order to further understand the interaction and the process of riverbank retreat.
First, mass failure, cantilever failure and fluvial erosion are respectively investigated by a series of hypothetical scenarios. The simulated results indicate that the occurrence of mass failure and cantilever failure mainly depend on the fluctuations in pore water pressure determined by river stage variations, soil permeability and rainfall condition. Fluvial erosion is determined by hydraulic conditions (i.e. channel slope, river stage and stage hydrograph) and soil erodibility (i.e. critical shear stress and erodibility coefficient). Subsequently, mass failure, cantilever failure and fluvial erosion are combined to estimate riverbank retreat. According to the conclusions of analyses, riverbank retreat is the process of repeated failure events and is primarily influenced by the magnitude and range of fluvial erosion with remarkable fluvial erosion. Finally, the results of the study reach of Jhuoshuei River reveal that the proposed model is capable of quantifying sediment yield and well predicting riverbank retreat length.
Keywords: riverbank retreat, unsaturated groundwater flow, limited equilibrium analysis, slope stability, mass failure, cantilever failure, fluvial erosion
目錄
誌謝 ... I 摘要 ... II Abstract ... III 目錄 ... V 表目錄 ... VIII 圖目錄 ... IX 符號說明 ... XIII 第一章 緒論 ... 1 1.1 研究動機與目的 ... 1 1.2 文獻回顧 ... 2 1.3 研究方法與步驟 ... 7 1.4 章節介紹 ... 8 第二章 分析模式之建立 ... 9 2.1 地下水流計算 ... 9 2.1.1 水流控制方程式 ... 9 2.1.2 土壤水力特性 ... 12 2.1.3 離散方程式與數值求解步驟 ... 14 2.1.4 模式驗證 ... 20 2.2 河岸退縮機制 ... 21 2.2.1 塊體破壞 ... 21 2.2.2 懸臂型破壞 ... 26 2.2.3 水流沖蝕 ... 26 2.3 模式演算流程 ... 30第三章 分析模式之測詴 ... 32 3.1 案例設定 ... 32 3.2 塊體破壞 ... 32 3.2.1 模式概述 ... 32 3.2.2 河川水位升降之影響 ... 33 3.2.3 降雨特性之影響 ... 37 3.2.4 河川水位與降雨特性之綜合影響 ... 40 3.3 懸臂型破壞 ... 41 3.3.1 河川水位升降之影響 ... 42 3.3.2 降雨特性之影響 ... 42 3.3.3 河川水位與降雨特性之綜合影響 ... 43 3.4 水流沖蝕 ... 43 3.4.1 河川臨界水位之探討 ... 43 3.4.2 水位歷線型態之影響 ... 44 3.5 河岸退縮 ... 44 3.5.1 案例設定 ... 45 3.5.2 水力條件之影響 ... 45 3.5.3 土壤條件之影響 ... 47 3.5.4 降雨特性之影響 ... 48 第四章 應用案例 ... 49 4.1 案例說明 ... 49 4.2 數值模式設定 ... 49 4.3 結果與討論 ... 51 第五章 結論與建議 ... 53 5.1 結論 ... 53 5.2 建議 ... 54
參考文獻 ... 56 附錄一 Hydrostatic Model 簡介 ... 63 附錄二 Boussinesq Model 簡介 ... 65
表目錄
表 1.1 相關研究彙整表 ... 74 表 2.1 等效摩擦角彙整表 ... 75 表 3.1 砂、壤土砂與黏壤土之未飽和土壤水力參數彙整表 ... 75 表 4.1 濁水溪主流河道斷面曼寧糙度係數率定結果彙整表 ... 75 表 4.2 濁水溪應用河段模擬結果彙整表 ... 76 表 4.3 濁水溪應用河段河岸退縮長度比較表 ... 76圖目錄
圖 1.1 河岸邊坡破壞型態示意圖 ... 77 圖 1.2 邊坡臨界高度與坡度關係示意圖 ... 77 圖 2.1 控制體積示意圖 ... 78 圖 2.2 土壤顆粒間氣液界面變化歷程示意圖 ... 78 圖 2.3 典型 S 型土壤保水曲線圖 ... 79 圖 2.4 二維 Richards 方程式計算概念圖 ... 79 圖 2.5 二維 Richards 方程式數值計算驗證結果圖... 80 圖 2.6 二維 Richards 方程式數值計算驗證結果圖... 80 圖 2.7 廣義莫爾-庫倫破壞包絡線示意圖 ... 81 圖 2.8 壓力水頭計算示意圖 ... 81 圖 2.9 FS 與 FSc計算概念圖 ... 82 圖 2.10 懸臂型破壞型態示意圖 (a)剪力型破壞 (b)樑型破壞 (c)張力型破壞 ... 83 圖 2.11 梯形渠道剪應力分布圖 ... 83 圖 2.12 邊界剪應力計算概念圖 ... 84 圖 2.13 沖蝕係數與臨界剪應力關係圖 ... 84 圖 2.14 沖蝕係數與臨界剪應力關係圖 ... 85 圖 2.15 分析模式演算流程圖 ... 86 圖 3.1 塊體破壞模擬參數及河岸幾何形狀概念圖 ... 88 圖 3.2 砂、壤土砂與黏壤土之 (a)保水曲線 (b)滲透函數 ... 88 圖 3.3 不同土質與河川水位上升速度之 FS 計算比較圖 ... 89 圖 3.4 砂質河岸於河川水位上升速度為 0.15m/hr 之地下水面線分布圖 ... 90 圖 3.5 砂質河岸於河川水位上升速度為 0.15 m/hr 沿破壞面壓力水頭分布圖 ... 91 圖 3.6 砂質河岸於河川水位上升速度為 0.15 m/hr 之壓力水頭與地下水流場分布圖 92圖 3.7 不同土質與水位下降速度之 FS 計算比較圖 ... 93 圖 3.8 不同河川水位下降速度之最小安全係數比較圖 ... 94 圖 3.9 不同降雨強度之平均入滲率比較圖 ... 95 圖 3.10 不同貣始水位之降雨強度與平均入滲率關係圖 ... 96 圖 3.11 不同貣始水位與降雨強度之 FS 計算比較圖 ... 97 圖 3.12 不同貣始水位之破壞包絡線分布圖 ... 98 圖 3.13 砂質河岸之壓力水頭與地下水流場分布圖-降雨強度為 (a)5 mm/hr (b)10 mm/hr ... 99 圖 3.14 砂質河岸之壓力水頭與地下水流場分布圖-降雨強度為 (a)20 mm/hr (b)40 mm/hr ... 100 圖 3.15 砂質河岸之壓力水頭與地下水流場分布圖-降雨強度為 (a)60 mm/hr (b)80 mm/hr ... 101 圖 3.16 黏壤土質河岸之壓力水頭與地下水流場分布圖-降雨強度為 (a)5 mm/hr (b)10 mm/hr ... 102 圖 3.17 黏壤土質河岸之壓力水頭與地下水流場分布圖-降雨強度為 (a)20 mm/hr (b)40 mm/hr ... 103 圖 3.18 黏壤土質河岸之壓力水頭與地下水流場分布圖-降雨強度為 (a)60 mm/hr (b)80 mm/hr ... 104 圖 3.19 不同累積降雨量與破壞包絡線分布比較圖 ... 105 圖 3.20 不同降雨強度與河川水位上升速度之 FS 計算比較圖 ... 106 圖 3.21 不同河川水位上升速度之破壞包絡線分布圖 ... 107 圖 3.22 不同降雨強度與河川水位下降速度之 FS 計算比較圖 ... 108 圖 3.23 不同河川水位下降速度之破壞包絡線分布圖 ... 109 圖 3.24 懸臂型破壞模擬參數及河岸幾何形狀概念圖 ... 110 圖 3.25 不同土質之 FSc計算比較圖 (a)河川水位上升 (b)河川水位下降 ... 111 圖 3.26 河川水位上升速度 0.15 m/hr 之壓力水頭與地下水流場分布圖 (a)砂質 (b)黏
壤土質 ... 112 圖 3.27 河川水位上升速度 0.9 m/hr 之壓力水頭與地下水流場分布圖 (a)砂質 (b)黏壤 土質 ... 113 圖 3.28 不同降雨強度之 FSc計算比較圖... 114 圖 3.29 不同降雨強度與河川水位上升速度之 FSc計算比較圖 ... 115 圖 3.30 水流沖蝕臨界水位分布圖 ... 116 圖 3.31 不同河川水位歷線型態沖蝕率與沖蝕長度比較圖 ... 116 圖 3.32 不同底床坡降之 FS 與 FSc計算比較圖(1/2) ... 117 圖 3.33 不同底床坡降之 FS 與 FSc計算比較圖(2/2) ... 118 圖 3.34 河岸退縮隨時間變化圖(1/2) ... 119 圖 3.35 河岸退縮隨時間變化圖(2/2) ... 120 圖 3.36 三種破壞機制隨時間變化圖 (a)泥砂產量 (b)泥砂累積產量與河岸退縮長度 ... 120 圖 3.37 不同流量之泥砂累積產量與河岸退縮長度比較圖 ... 121 圖 3.38 不同水位歷線型態之泥砂累積產量與河岸退縮長度比較圖 ... 121 圖 3.39 臨界剪應力為 10 Pa 之泥砂累積產量與河岸退縮長度比較圖 ... 122 圖 3.40 臨界剪應力為 100 Pa 之泥砂累積產量與河岸退縮長度比較圖 ... 122 圖 3.41 砂質河岸之泥砂累積產量與河岸退縮長度比較圖 ... 123 圖 3.42 黏壤土質河岸之泥砂累積產量與河岸退縮長度比較圖 ... 123 圖 3.43 不同降雨強度之泥砂累積產量與河岸退縮長度比較圖 ... 124 圖 3.44 臨界剪應力為 10 Pa 不同降雨強度之泥砂累積產量與河岸退縮長度比較圖 ... 124 圖 4.1 濁水溪應用案例河段相關位置圖 ... 125 圖 4.2 土壤剪力強度圖 (a)應力路徑圖 (b)莫爾圓圖 ... 125 圖 4.3 jest-test 詴驗結果圖 (a)詴驗一 (b)詴驗二 ... 126 圖 4.4 jest-test 詴驗結果套疊比較圖 ... 126
圖 4.5 海棠颱風水位率定結果圖 (a)溪洲大橋 (b)自強大橋 ... 127 圖 4.6 碧利斯颱風水位驗證結果圖 (a)溪洲大橋 (b)自強大橋 ... 127 圖 4.7 柯羅莎颱風水位驗證結果圖 (a)溪洲大橋 (b)自強大橋 ... 128 圖 4.8 民國 94 至 96 年彰雲大橋流量歷線圖 ... 128 圖 4.9 斷面編號 43 河岸退縮歷程圖 ... 129 圖 4.10 斷面編號 44 河岸退縮歷程圖 ... 129 圖 4.11 斷面編號 45 河岸退縮歷程圖 ... 130 圖 4.12 八場颱洪事件之泥砂累積產量與河岸退縮長度比較圖 ... 130 圖 4.13 民國 93 年灘線位置比較圖 ... 131 圖 4.14 民國 96 年灘線位置比較圖 ... 131
圖 A1.1 Hydrostatic Model 孔隙水壓計算示意圖 ... 64
圖 A1.2 Hydrostatic Model 計算流程圖 ... 64
圖 A2.1 一維控制體積示意圖 ... 71
圖 A2.2 一維 Dupuit assumption 示意圖 ... 71
圖 A2.3 一維 Boussinesq 方程式計算概念圖 ... 72
圖 A2.4 一維 Boussinesq 方程式數值計算驗證結果圖 (a)第 24 hr (b)第 48 hr ... 72
圖 A2.5 Boussinesq Model 孔隙水壓計算示意圖 ... 73
符號說明
A：水流通過截面積 AB、BC、CD、DF 、FG、GH、HA：河岸土體邊界範圍 a：沖蝕率計算經驗參數 ai、bi、ci、di：數值計算格點之權重係數 aij、bij、cij、dij、eij、fij：數值計算格點之權重係數 C：比含水量 c'：土壤有效凝聚力 D：泥沙顆粒直徑 50 D ：中值粒徑 FD：驅動力 FR：抗剪力 FS：塊體破壞安全係數 FSc：懸臂型破壞安全係數 FSmin：最小安全係數 H：河岸高度 Hb：土層底部厚度 h：水力水頭
t h0 、hl
t ：邊界上隨時間變動之水頭 hgwb：土層內部邊界上之地下水面高 i h：河岸邊坡土層之貣始水頭 hp：破壞面上任一點距地下水面線的垂向距離 w h ：河川水位高 hwc：河川臨界水位hwi：貣始河川水位 hwp：洪峰水位 I：總土體單元數 Ir：降雨強度 i、j：空間網格座標 K、K
、K
：水力傳導係數 d k ：沖蝕係數 s K ：飽和水力傳導係數 L：沿破壞面之土層總長度、延破壞面方向 Le：沖蝕長度 Lme：最大沖蝕長度 Lr：河岸退縮長度 Lsat：沿破壞面之飽和土層長度 Luns：沿破壞面之未飽和土層長度 l：計算域長度 M、N：數值計算格點數 m：Picard 迭代次數 v m 、nv：曲線密合參數 N：總正向力 Ni：源流項 n：時間網格座標 P：河川水位靜水壓力 Pt：總累積降雨量 Q：流量 x Q ：x 方向流量 q、q 、q 、q ：水流通量R：水力半徑 Ra：地表平均入滲率 Ri：地表入滲率 S：基質吸力 S0：底床坡降 Sf：能量坡降 y S ：比出水量 s：流線方向 T：總剪力 Tsat：沿破壞面之飽和土層剪力 Tuns：沿破壞面之未飽和土層剪力 t：時間、降雨延時 U：上舉力 ua：孔隙氣壓 uw：孔隙水壓
ua uw
：基質吸力 V：流速 Vas：泥砂累積產量 Vs：泥砂產量 Vw：河川水位升降速度 W：破壞土體重 x：空間座標 y：沿河岸邊坡方向的距離 z：空間座標、位勢水頭 j i z, ：z 方向高程 Δt：時間間距Δx、Δz：為空間間距 ：有效飽和度 ：河岸坡面角度 v
：進氣潛能因子 β：破壞面角度 b
：等分線與底床之夾角 d
：土壤乾單位重 m
：土壤濕單位重 s
：土壤單位重 w
：水單位重 ：數值計算收斂條件 ε：沖蝕率 ：含水量 r
：殘餘含水量 s
：飽和含水量 w ：地下水面線坡度 v ： 2 h w ：水的密度 σ：正向應力
ua
：淨正向應力
uw
：有效正向應力 b ：邊界剪應力 c
：臨界剪應力 * c ：無因次臨界剪應力 cb
：底床臨界剪應力 τ ：飽和土壤剪應力τuns：未飽和土壤剪應力 ：摩擦角 ：土壤有效摩擦角 b ：基質吸力造成剪應力增加所對應之角度 * ：等效摩擦角 ：壓力水頭
t 0
、
l
t ：邊界上隨時間變動之壓力水頭 h ：水力梯度第一章 緒論
1.1 研究動機與目的
台灣地區地狹人稠，土地高度開發利用，以致多採築堤束水治理工法以降低水患 之威脅。然而，台灣特殊的氣候與地質條件，河川大多流短而急，以致暴雨來臨時， 洪水歷程急促，水流沖刷劇烈，除了垂向侵蝕導致河床嚴重刷深外，側向侵蝕造成的 河岸退縮，不但影響兩岸土地之安定，亦可能危及堤防和水工構造物之安全，更甚者， 造成人民生命財產的嚴重損失。另外，隨著數值工具的發達，目前國內河川治理規劃 大多採用數值計算模型來代替水工模型詴驗進行水理輸砂的分析工作，迄今所累積的 研究成果與工程經驗雖有相當程度的可靠性與準確性，然綜觀目前常用的河道動床數 值模式則多著重於底床沖淤之預測，對於側向變遷之計算，仍未臻完備。因此，對河 岸退縮的可能型態、形成機制及演變過程等實須進一步研析與瞭解。 一般而言，河岸退縮之趨勢，可利用現場調查量測、經驗法則或數值計算模式等 方法進行評估。雖然透過回歸相關地文、水文條件與量測資料可達到預測河岸退縮之 目的，但此種方法需要大量長期的現場調查資料且缺乏嚴謹之理論基礎。另外，經驗 法係利用流域面積、河寬、流量、底床坡降、河床質等資料，藉由回歸分析建立其關 係式，然而，此法可能因時因地而存在適用性的問題，且時間尺度侷限在長期河寬之 預測，並無法得知其隨時間變化的情形。而數值計算模式則是根據物理控制方程式並 以數值方法進行求解，因其具泛用性且能定量分析處理複雜的問題，因此為一準確且 有效率之方法。故本研究之目的，從較具物理意義的力學理論切入，發展一套預測河 岸退縮行為之數值計算模式，以深入探究洪流來臨時，河岸遭受破壞及其退縮過程， 以做為高灘地管理或河防構造物安全評估之工具，更期能應用於水道治理計畫線、河 川區域線或堤防預定線劃設之參考依據。1.2 文獻回顧
天然河川總是處在不斷地變化與發展的過程中，在空間上主要以兩種型態呈現： 縱向變形和側向變形。縱向變形為泥砂與水流交互作用下的底床沖淤行為，側向變形 為近岸流速的改變導致河岸淤長(riverbank advance)或河岸退縮(riverbank retreat)。造成 河岸淤長主因如：主流偏離、近岸為緩流或局部渦流致使泥砂沉滓落淤，而造成河岸 退縮主因通常是河岸邊坡土體與近岸水流相互作用下的結果，河岸邊坡表土與坡趾遭 受沖刷，再者使其坡度變陡或岸高增加，最終因重力作用而產生崩塌或滑動破壞，河 岸退縮即為如此反覆循環作用下之結果。河岸邊坡受河川水流作用的受蝕程度與其組 成的材料性質相關，Johnson and Stypula (1993)依據河岸所構成的材料將其分為岩床 (bedrock)、非凝聚性(cohesionless)、凝聚性(cohesive)與層狀(stratified or interbedded)等 四類。岩性河岸大多將之視為河床整體的一部份進行探討，而岩床的侵蝕是眾多因素 交互作用下的結果，如：磨蝕(abrasion)、塊體抽離(plucking)、穴蝕(cavitation)、顆粒 彈跳沖蝕(saltation)、風化侵蝕(weathering)等(Whipple et al., 2000；Sklar and Diereich, 2004；Stock, et al., 2005)，其過程與機制非常複雜，欲進行完全定量之預測確有其困難。 因此，目前主要分析方法依其破壞機制大致分為水力沖蝕(hydraulic scour)以及泥砂磨 蝕(sediment and rock abrasion)等兩種類型進行岩床侵蝕速率的評估(Annandale, 1995； Sklar and Diereich, 2004, 2006；Lamb et al., 2008)。相較於其他河岸材料，岩性河岸一 般相對穩定而不易發生大規模的河岸退縮。
相對於岩性河岸，Julien(2002)將沖積型河川(alluvial river)的破壞依其邊坡材料性 質分為非凝聚性、凝聚性與層狀等三種，並指出非凝聚性河岸其坡趾遭受沖刷而使其 坡度大於其安息角(angle of repose)時將會產生滑動破壞。Johnson and Stypula (1993)則 描述非凝聚性河岸多由粉砂(silt)、砂(sand)或礫石(gravel)所摻混組成，由於該類型的材 料顆粒間缺乏鍵結能力，因此，其破壞型態可視為多個單一顆粒受水流沖刷作用而脫 離的行為，其破壞程度主要是受材料特性(顆粒的大小、形狀、級配、含水量與相對密
度)以及水力條件(近岸流速、紊流強度與剪應力)而定。Thorne and Tovey (1981)發現此 類型的河岸邊坡多為粗顆粒所組成而透水性較佳，其穩定與否受土體內的孔隙水壓 (pore water pressure)變化影響較小，河川水流的作用才是造成破壞的主因。Thorne (1991)則提出非凝聚性材料的破壞類型其一為平面或圓弧形的淺層剪力破壞，其二為 河川水流挾帶泥砂顆粒的表層剝離。Terzaghi and Peck (1996)則指出隨著深度增加，因 剪力增量較剪應力大，使得此類型的破壞面大多發生於淺層。上述文獻僅對非凝聚性 河岸破壞的特性與行為進行論述，而 ASCE Task Committee (1998)則綜合相關研究後指 出此類型的河岸穩定分析方法，主要分為兩種：Pizzuto (1990)與 Li and Wang (1993, 1994)假設當河岸邊坡坡度大於其安息角時，則土體將沿破壞面產生滑動；Wiele (1992) 與 Kovacs and Parker (1994)則建構泥砂連續方程式，以側向通量處理河川水流侵蝕的 過程。然而，Thorne (1991)則認為天然河岸邊坡通常摻雜具有凝聚性的土壤，並且須 考量由植物根系或未飽和土層所提供額外的凝聚性，因此，非凝聚性河岸邊坡的應用 範圍因實用性而有所限制。
對於凝聚性河岸邊坡而言，Johnson and Stypula (1993)說明此類型的河岸因富含黏 土(clay)，相較於非凝聚性河岸，其滲透性較差，因此其破壞機制與土體內的地下水息 息相關，如：滲流(seepage)、管湧(piping)、凍脹(frost heaving)作用等，並提出此類型 的河岸邊坡破壞常發生於河川水位洩降的過程。Hagerty (1991)則根據河岸邊坡土質特 性與坡度差異提出不同的破壞類型，如圖 1.1 所示。圖 1.1(a)與圖 1.1(b)分別為發生於 均質土壤且陡峭邊坡的平面型(planar)以及低緩邊坡的圓弧型(rotational)塊體破壞(mass failure)，並且指出破壞常發生於降雨或河川水位的漲退變化過程中，河岸土體濕潤而 增加其自重並減少土壤顆粒間的鍵結強度則為其肇因。而圖 1.1(c)所示則為河岸底層 土壤受水流侵蝕作用下，致使上方土層形成類似懸臂樑形式之土塊並因重力作用而發 生懸臂型破壞(cantilever failure)，Thorne and Tovey (1981)則依據不同的破壞機制，提 出三種形式及其分析方法，後續則有 Amiri-Tokaldany and Samadi (2007)、Langendoen and Simon et al. (2008)、Simon et al. (2009)與 Simon et al. (2010a)等研究結合此類型的破
壞機制建構其分析模式，並指出欲探討完整的河岸退縮過程，此類型的破壞機制必須 加以考量方能符合實際情形。
前述針對凝聚性河岸之塊體破壞，依據破壞面的類型大致可分為平面與弧面的破 壞(見圖 1.1(a)與圖 1.1(b))，然而根據許多現場調查資料顯示凝聚性河岸的破壞型態大 多為沿平面破壞面滑動的塊體形式(Thorne, 1982；Osman and Thorne, 1988；Darby and Thorne, 1996)，因此相關研究多採此種破壞型態進行解析。分析方法則多參照大地工 程邊坡穩定的理論架構，亦即在均質土體且不考慮其變形下，以莫爾-庫倫破壞準則 (Mohr-Coulomb failure criterion)建立河岸土體沿破壞面之正向應力與剪力強度關係，並 依據土體的抗剪力與驅動力，採用極限平衡分析法(limited equilibrium analysis)評估其 穩定性。基於上述分析理論，Lohnes and Handy (1968)以坡度與剪力參數建立臨界坡高 的關係式。Chen (1975)同樣以臨界坡高來評估邊坡的穩定性，並依不同河岸坡度彙整 結果如圖 1.2 所示。Osman and Thorne (1988)則結合 Arulanandan et al. (1980)所提出的 水流剪應力與臨界剪應力關係式，計算河岸邊坡受河川水流沖蝕作用下的退縮長度， 並據以評估河岸受水流沖蝕後對河岸穩定之影響。然而，上述研究均未考量河川水位 變化以及河岸土體孔隙水壓之影響，其研究結論顯示河岸破壞發生與否僅與邊坡的幾 何條件以及土壤強度有關。因此，後續學者則致力改善前述分析理論不足之處，Huang (1983)根據破壞土體內地下水位的分布比例提出孔隙水壓比(pore pressure ratio)，據以 評估其對邊坡穩定的影響。Springer et al. (1985)則分析一砌形土體受河川水位與地下水 位對其穩定之影響。Simon et al. (1991)在簡單的河岸幾何形狀條件下，將河川水位與 地下水的影響機制同時納入其分析模式。Darby and Thorne (1996)則基於前述的研究， 並進一步考量河川水流沖蝕作用下對河岸穩定之影響，且與前人分析方法比較其差異 性。然而，上述分析方法對於河岸土體內之孔隙水壓，僅考量地下水面線以下飽和區 域所產生的上舉力(uplift force)，而忽略未飽和區域形成基質吸力(matric suction)的影 響。為釐清土壤未飽和狀態之機制，Casagli et al. (1999)與 Rinaldi and Casagli (1999)依 據 Fredlund et al. (1978)所提出之廣義莫爾-庫倫破壞準則(extended of Mohr-Coulomb
failure criterion)，進一步結合基質吸力分析義大利 Sieve River 河川水位與河岸土體內 孔隙水壓的變化情形，研究結果顯示基質吸力對河岸邊坡穩定的貢獻甚大，尤其在地 下水位偏低的狀況下。然而，天然沖積型的河岸，其土壤性質多呈現非均質的層狀分 布，Simon et al. (2000)遂將不同的土壤強度參數建構至所發展的模式以反映層狀土壤 的效應，另依據美國 Goodwin Creek 長期監測的資料，評估河岸穩定程度隨時間的變 化情形，並探討基質吸力在河川水位漲退過程中的變動及其對河岸穩定的影響，最後 在 結 論 中 強 調 基 質 吸 力 對 河 岸 穩 定 的 重 要 性 。 而 由 美 國 農 業 部 農 業 研 究 局 (USDA-ARS)所發展並被廣為應用之河岸邊坡穩定分析模式(Bank Stability and Toe Erosion Model, BSTEM)亦以上述分析理論為其主要架構，爾後相關研究大多以該模式 或依循上述分析理論進行現地案例之分析與應用，如：Simon et al. (2002)以美國 Missouri River 為例，探討不同河川水位條件下河岸的臨界穩定狀況；Simon and Thomas (2002)在美國 Yalobusha River 以不同地下水位高度的條件下，建立河岸高度與臨界坡 度的關係；Amiri-Tokaldany et al. (2003)則以美國 Hotophia Creek 為研究對象，探討河 岸的長期破壞趨勢；Amiri-Tokaldany and Darby (2006)則分別分析義大利 Sieve River 與美國 Goodwin Creek 的河岸穩定情形。然而，上述研究僅針對塊體破壞在河川水位 變化單一機制下進行分析與探討。
後續學者則進一步結合 Partheniades (1965)所提出的沖蝕率計算方法探討河岸破 壞與退縮行為，如：Langendoen and Simon (2008)模擬美國 Goodwin Creek 河岸幾何形 狀受破壞而變化的過程；Simon et al. (2009)則評估美國 Blackwood Creek、Upper Truckee River、Ward Creek 等多條河川受蝕的情形；Shields el at. (2009)則分析美國 Little Topashaw Creek 河岸在不同的地下水位條件下的穩定狀況；Simon et al. (2010a)在美國 Big Sioux River 進行不同重現期距流量下泥砂產量之計算。然而，前述研究的分析理 論架構，皆將地下水面線簡化為一水平狀態，並使地下水位維持在河川水位上升歷程 的最高水位以評估河岸塊體破壞的臨界條件，此假設雖可快速分析河岸穩定與否，但 卻無法真實反映河川水位升降、土壤特性或降雨因素造成孔隙水壓變化對破壞的影
響。有鑑於此，Chiang et al. (2011)藉由求解一維非穩態 Boussinesq 地下水流方程式， 將河川水位視為邊界條件以考量河川水位與地下水位的互制關係，利用求得之地下水 面線計算孔隙水壓分布情形，據以評估河岸之穩定性。另有部分學者為能更精確掌握 孔隙水壓的變化情形，則利用商用軟體 SEEP/W(Geo-Slope International Ltd)模擬河岸 土體內孔隙水壓的分布情形進行現地案例的分析，如：Dapporto et al. (2001)與 Dapporto et al. (2003)評估義大利 Arno River 因河川水位變化的河岸穩定狀況，並歸納最高河川 水位與河岸穩定之關係；Rinaldi et al. (2004)模擬一場完整的洪水事件探討義大利 Sieve River 河岸穩定隨時間變化的情形，並提出河岸穩定條件與洪水型態、河岸土質及其幾 何形狀有關。然而，上述研究皆僅限於塊體破壞之分析與預測，且降雨的影響亦未加 以考量與探討。
影響河岸退縮的因素眾多，Thorne (1982)、Lawler (1995, 1997)認為河岸退縮主要 包含三個破壞機制之集合：表面侵蝕(subaerial processes)、塊體破壞(mass failure)與水 流沖蝕(fluvial erosion)。其中表面侵蝕係指因各種氣候因素造成土壤表面弱化、風化等 過程，相較於其他兩個機制，表面侵蝕對於河岸退縮的影響甚微，如：Lawler (1993) 指出英國 River Ilston 年侵蝕率約 27 mm；Prosser et al. (2000)調查澳洲 Ripple Creek 年 侵蝕率為 132 mm；Couper and Madock (2001)於英國 River Arrow 數個監測點發現年 侵蝕率多介於 10 至 40 mm，最大達 181 mm；Veihe et al. (2011)於丹麥 Harrested Stream 測得年侵蝕率為 17.3 至 30.1 mm。因此，Rinaldi and Darby (2008)認為塊體破壞與水 流沖蝕才是造成河岸退縮的主要因素。Thorne (1982)以及 Darby et al. (2010)則指出欲 探討長期的河岸退縮行為，河川水流造成坡趾的掏刷為其關鍵因素。另外，Hagerty (1991)、Simon and Curini (1998)、Casagli et al. (1999)、Rinaldi and Casagli (1999)、Simon et al. (2000)等則提出降雨亦為造成河岸破壞的主要原因之一，然而，過去相關研究甚 少針對其進行分析與討論。
(1975)計算概念為基礎並與 NETSTARS 一維動床模式結合，分析河道斷面在颱洪事件 前後的變化情形；林恩添(2005)採用 Osman and Thorne (1988)評估方法，在不考慮地下 水的作用下與二維動床模式結合，以探討河道垂向沖淤與側向變遷的問題；陳晉琪 (2007)控制水位在緩慢上升的條件下進行室內渠槽詴驗，並歸納河岸臨界破壞條件， 最後與 Chen (1975)分析方法之計算結果相互驗證；駱建利(2009)以室內詴驗進行降雨 機制對邊坡破壞的研究，並與 SEEP/W 數值模擬結果進行比較；黃群玲(2010)則結合 一維非穩態 Boussinesq 地下水流方程式以及 Green-Ampt 入滲理論，藉以評估降雨入 滲對河岸塊體破壞的影響。 茲將本節所述相關研究彙整於表 1.1，由表 1.1 可綜覽過去研究之發展過程與未盡 完備之處。
1.3 研究方法與步驟
如前節所述，過往相關研究尚有考量未盡周詳之處，且多數研究以現場案例為其 分析主軸，僅能就單一案例予以分析與探討，而並未通泛地進行不同誘發條件對河岸 破壞以及退縮機制之研析與歸納。因此，本研究擬以較為嚴謹之未飽和地下水流理論 計算河岸土體內部隨時間變化的孔隙水壓分布，並利用邊界條件的處理技巧，將河川 水位變動以及降雨因素之影響納入考量，並整合塊體破壞、崩懸臂型破壞與水流沖蝕 三種主要破壞機制，建立河岸退縮行為數值計算模式，期能藉由模式的測詴與應用， 對於河岸退縮過程的現象能有進一步的瞭解。 本研究方法與步驟首先闡述本研究所應用之相關理論基礎，包含地下水流計算以 及塊體破壞、懸臂型破壞、水流沖蝕等河岸退縮機制，並說明模式整合方法與演算流 程。接續應用所建立之模式於不同滲透性土壤之河岸，以河川水位變化和降雨強度個 別探討對上述三種破壞機制之影響，再以假設案例，評估在不同水力、土壤與降雨條 件下，河岸受破壞所造成的泥砂產量以及河岸退縮長度，並將模式應用於濁水溪河段，藉由現場案例之模擬展示模式之實用性並對模式不盡完善之處進行討論。最後，對本 研究分析結果統整結論，並提出建議。
1.4 章節介紹
本研究擬以五個章節探討河岸退縮的過程與行為，茲將本文各章節內容扼要說明 如下： 第一章為緒論，主要闡述本研究的研究動機與目的，同時回顧過去相關研究，並 提出本研究所採用之方法與步驟，最後簡要說明本文主要內容。 第二章為分析模式之建立，針對本研究所採用之分析理論以及建立方法進行介 紹，主要包含地下水流計算以及河岸破壞機制兩大部分，其相關的控制方程式、數值 求解步驟、邊界條件處理、假設條件與模式演算流程等均詳述於本章。 第三章為分析模式之測詴，主要分為三部分，首先以本研究所建立之模式與前人 相關研究進行比較，接續針對塊體破壞、懸臂型破壞以及水流沖蝕三種破壞機制，個 別探討不同條件對其之影響，最後，整合上述三種機制，以假設的案例進行泥砂產量 與河岸退縮長度之分析與探討。 第四章為應用案例，以所建立之模式應用於濁水溪河段，藉由實際案例之模擬展 示模式之實用性。 第五章為結論與建議，除對本研究成果做綜合性之歸納說明外，並對不盡完備或 供未來可改進方向提出建議。第二章 分析模式之建立
本章首先介紹地下水流連續方程式以及描述地下水分別在飽和與未飽和土壤孔隙 中運動所依循 Darcy 定律與 Darcy-Buckingham 方程式之相關理論，據以推導二維 Richards 方程式，接續說明土壤水力特性、離散方程式與數值方法求解步驟，並以簡 單的案例驗證所建立的地下水流數值計算模式。其次概述所考量的河岸退縮機制，包 含塊體破壞(mass failure)、懸臂型破壞(cantilever failure)與水流沖蝕(fluvial erosion)等三 種類型之理論與分析方法。最後，說明模式整合與建立的步驟以及演算流程。
2.1 地下水流計算
2.1.1 水流控制方程式 圖 2.1 為控制體積示意圖，假設土壤中水份為連續分布，其移動遵循質量守恆(mass conservation)定理，則可推得地下水流連續方程式為：
CV CS w wd V ndA t ˆ (2.1)式中，CV 表示控制體積(control volume)；CS 表示控制表面(control surface)；t 表示時 間；w為水的密度；V為流速；等號右側負號表示流入控制表面。上式的物理意義係
描述地下水在控制體積內的質量增加率等於進入控制體積內的質量流率。
Darcy 定律
地下水最基本的運動形式為滲透，亦即受重力或壓力差的作用下通過孔隙介質的 流動行為。在重力作用下，水由位勢水頭(z, elevation head)大向小的方向流動；在壓力 差作用下，水由壓力水頭( , pressure head)大向小的方向流動，而水力水頭(h, hydraulic head)為上述位勢水頭與壓力水頭之和，因此，地下水的流動方向取決於水力水頭的水
頭差，即水力梯度(hydraulic gradient)。法國水利工程師 Henry Darcy 透過均質飽和砂 濾層的滲透詴驗，於 1856 年提出水流通過飽和孔隙介質的流量與水力梯度關係，即 Darcy 定律(Darcy’s Law)：
h KA
Q (2.2)
z
h (2.3) 式中，Q 為流量(flow rate)；K 為水力傳導係數(hydraulic conductivity)；A 為水流通過 的截面積；h 為水力水頭； h 為水力梯度； 為壓力水頭；z 為位勢水頭。 式(2.2)亦可表示如下： h K q (2.4) 式中，q 為水流通量(flux)，亦稱為達西流速(Darcy velocity)。式(2.2)與式(2.4)中，等號 右側負號表示地下水朝水力梯度遞減的方向流動。 Darcy-Buckingham 方程式 水在飽和土壤中的運動情形可依循 Darcy 定律加以描述，而探討未飽和土壤中的 水流行為，土壤的孔隙比(void ratio)以及飽和度(degree of saturation)或含水量(water content)將變成影響地下水流動的重要因子。如圖 2.2 所示，圖中編號由 1 至 7 氣液介 面的變化，代表土壤由飽和轉為未飽和的過程，當含水量逐漸減少時，土壤顆粒間的 孔隙將逐漸由空氣所取代，能容許水份運動的孔隙則逐漸減少且流動通道亦逐漸曲 折，而不利地下水的流動。換言之，土壤滲透性隨含水量減少而降低。Buckingham (1907) 修正 Darcy 定律，提出水力傳導係數並非為常數，而是壓力水頭或土壤含水量之函數， 則未飽和水流通量與水力梯度的關係可表示如下：
h K q (2.5)
h K q (2.6) 式中， 為含水量(water content)；K
與K
表示未飽和土壤之水力傳導係數分別為 壓力水頭以及含水量之函數。 二維 Richards 方程式 若僅考慮 x 與 z 方向水流，則(2.1)式可分別表示如下：
t dxdz dxdz t d tCV w w
(2.7)
z q dxdz x q dxdz q dx dz z q q q dz dx x q q dA n V z w x w CS z w z z w x w x x w w ˆ (2.8) 式中，x 與 z 為空間座標；qx與qz分別表示 x 與 z 方向之通量。將式(2.7)與式(2.8)分別 代入式(2.1)後可得： z q dxdz x q dxdz t dxdz z w x w w (2.9) 本研究以壓力水頭 為單一變數，且考量未飽和土壤中水流的運動行為，依據 Darcy-Buckingham 方程式，式(2.9)中qx與qz可分別表示為：
x K x z K h K qx (2.10)
1 z K z z K h K qz (2.11)將式(2.10)與式(2.11)代入式(2.9)後整理可得：
z K z K z x K x t (2.12) 上式即為二維 Richards 方程式(Richards,1931)，亦是被廣泛地用以探討未飽和土壤中水 流運動的控制方程式。欲求解式(2.12)須給予適當之貣始條件(initial condition)及邊界條 件(boundary condition)，上述條件將在 2.1.3 節與數值求解過程中一併說明。 2.1.2 土壤水力特性 式(2.12)中，土壤含水量與水力傳導係數均為壓力水頭之函數，因此，上述的關係 式須先建立方能進行求解二維 Richards 方程式。土壤保水曲線(soil water retention curve; soil water characteristic curve)係描述土壤含水量與壓力水頭的關係，圖 2.3 為一 典型的土壤保水曲線，圖中
s為飽和含水量(saturated water content)，當孔隙壓力下降 時，含水量隨之減少，最終當壓力持續下降但含水量未有明顯變化，此時含水量即為 殘餘含水量(
r, residual water content)。上述未飽和土壤含水量的變化關係常以有效飽和度(, effective saturation)來描述，其定義如下： r s r (2.13) 式中，為有效飽和度；
r為殘餘含水量；
s為飽和含水量。 當土壤含水量達完全飽和時，則此時有效飽和度為 1，若土壤呈未飽和狀態，則 有效飽和度將視土壤含水量、殘餘含水量與飽和含水量三者而定。探討未飽和土壤保 水曲線的相關研究甚多(Leong and Rahardjo, 1997b)。本研究採用較為廣泛應用由 van Genuchten (1980)所提出的 S 形分布的土壤保水曲線關係式，其具有曲線平滑、相關參數取得容易，且曲線於土壤接近飽和時能滿足 0 之優點，該關係式如下：
0 / 1 1 , 1 1 0 1 v v m n v n m v v (2.14)式 中 ，
v為 進 氣 潛 能 因 子 (air entry value);nv 與mv 為 曲 線 密 合 參 數 (curve fitting parameters)。基於上述的定義，則土壤含水量可表示如下：
0 0 r s r s (2.15) 水力傳導係數為描述水在土壤內流動的重要參數，一般而言，與土壤種類、粒徑 分佈、孔隙率、總體密度、飽和度等土壤物理性質有關。當土壤含水量達飽和時，水 份佔據土壤孔隙而可在孔隙間自由流動，此時水力傳導係數為最大值且維持一定值， 稱為飽和水力傳導係數(Ks, saturated hydraulic conductivity)。但當土壤呈未飽和狀態 時，部分的土壤孔隙被空氣所佔據，可供流動的孔隙相對減少，且孔隙間由水與空氣 界面形成的毛細壓力(capillarity)亦限制了水份於孔隙間的移動，且隨著飽和程度降 低，水份流動受到毛細壓力的作用越強，導致水力傳導係數下降。滲透函數(permeability function)即為描述土壤含水量與水力傳導係數的關係，該函數利用理論或由實驗歸納 所得之半經驗公式研究頗多(Leong and Rahardjo, 1997a)，本研究採 van Genuchten (1980)改良自 Mualem (1976)所提出之關係式：
0 / 1 1 , 1 1 0 2 1 2 1 v v m m s s n m K K K v v (2.16)式中，Ks為飽和水力傳導係數; nv與mv為曲線密合參數。 2.1.3 離散方程式與數值求解步驟 考量二維 Richards 方程式具有高度非線性的特性，本研究依據 Clement(1994)所提 出之數值方法架構進行模式之建置，其具有易於建構、減少質量守恆誤差、數值穩定 等之優點，以下說明其步驟與流程。 以 Euler 後項差分法結合 Picard 迭代法對式(2.12)時間項進行離散化： t t n j i n m j i , 1 , 1 ,
(2.17) 式中，i、j 為空間網格座標；m 表示 Picard 迭代次數；n 為時間網格座標； t 為時間 間距。為解決應變數過多的問題，將上式中 1, 1 , n m j i 對 以一階泰勒級數展開可得：
m n
j i n m j i n m j i n m j i n m j i d d 1, , 1 , 1 , , 1 , , 1 , 1 , 1 ,
(2.18)另外，定義 C 為比含水量(specific water capacity)如下所示：
m n j i d d C , 1 ,
(2.19) 結合式(2.17)、式(2.18)與式(2.19)並整理後可得： t C t t n m j i n m j i m n j i n j i n m j i , 1 , 1 , 1 , , 1 , , , 1 , (2.20)式(2.12)空間項則以全隱式有限差分法(fully implicit finite difference method)結合 Picard 迭代法進行離散化：
2 2 1 2 2 1 2 2 1 , 1 1 , , 1 , , 1 1 , , 1 , 1 , 1 1 , 1 , 1 , , 1 1 , , 1 , 1 , 1 , 1 , 1 1 , , 1 1 , , 1 , 1 , 1 , 1 1 , 1 , , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i K K K K z z K K z K K z x K K x K K x z K z K z x K x (2.21) 式中， x 與z為空間間距。結合式(2.20)與式(2.21)並重新排列整理後可得參考座標 點( 1, 1 , n m j i )與其鄰近格點的關係式如下： ij n m j i ij n m j i ij n m j i ij n m j i ij n m j i ij b c d e f a ,11, 1 1,1, 1 ,1, 1 1,1, 1 ,11, 1 (2.22) 式中，aij、bij、cij、dij、eij、fij 分別表示鄰近格點之權重係數，係數 aij、bij、cij、dij、 eij、fij可表示如下：
, 1
1 , 1 , , 2 2 1 mn j i n m j i ij K K z a (2.23)
, 1
, 1 1 , , 2 2 1 mn j i n m j i ij K K x b (2.24)
, 1
, 1 1 , , 2 2 1 mn j i n m j i ij K K x d (2.25)
, 1
1 , 1 , , 2 1 mn j i n m j i ij K K e (2.26)t C e d b a c n m j i ij ij ij ij ij 1 , , (2.27)
, 1
1 , 1 , 1 , 1 , , 1 , , , 1 , , 2 1 mn j i n m j i n m j i n m j i n j i n m j i ij K K t t C t f
(2.28) 若將土層在 x 與 z 方向分別以 M 與 N 個格點進行細分成
M 1
N1
個元素， 則非邊界格點與其鄰近格點的關係式可表示如下： 1 1 1 , 1 , 1 1 1 1 , 1 1 , 1 1 1 , 1 1 , 1 1 1 1 , 1 1 , 2 1 1 1 , 1 2 , 1 1 1 23 1 , 1 4 , 2 23 1 , 1 3 , 3 23 1 , 1 3 , 2 23 1 , 1 3 , 1 23 1 , 1 2 , 2 23 32 1 , 1 3 , 3 32 1 , 1 2 , 4 32 1 , 1 2 , 3 32 1 , 1 2 , 2 32 1 , 1 1 , 3 32 22 1 , 1 3 , 2 22 1 , 1 2 , 3 22 1 , 1 2 , 2 22 1 , 1 2 , 1 22 1 , 1 1 , 2 22 N M n m N M N M n m N M N M n m N M N M n m N M N M n m N M N M n m n m n m n m n m n m n m n m n m n m n m n m n m n m n m f e d c b a f e d c b a f e d c b a f e d c b a (2.29) 貣始條件如圖 2.4 所示，貣始條件之設定採用靜水壓分布(hydrostatic pressure distribution) 之假設如下表示： j i i b m j i H h z, 0 , 1 , (2.30) 式中，Hb為土層底部厚度；hi為貣始地下水位；zi,j為 z 方向高程，下標 i、j 為空間 網格座標。 邊界條件 如圖 2.4 所示，本研究應用於 HA、AB、BC 與 FG 為固定水頭邊界條件(Dirichlet
boundary condition)，可分別表示如下： HA、AB、BC： im,j1,n1 Hb hw
t zi,j (2.31) FG： Mm,1j,n1 Hb hgwb
t zi,j (2.32) 式中，Hb為土層底部厚度；hw
t 為河川水位高；hgwb
t 為土層內部邊界之地下水位高。 若 FG 為一不透水邊界，則 x 方向流量邊界條件係以 Darcy 定律描述水在土壤中的流 動行為，故依據式(2.5)，則 FG 零流量邊界條件可表示為： FG：
x K x z K h K 0 (2.33) 上式以差分式可表示為： FG： 0 1 , 1 , 1 , 1 , 1 x n m j M n m j M (2.34) 上式重新排列後可得： FG： 1, 1 , 1 1 , 1 , m n j M n m j M (2.35) 土層表面(CD、DF)在積水效應(ponding effect)尚未發生時，降雨可完全入滲至土 體，故此時地表入滲率(Ri, infiltration rate)將等於降雨強度(Ir, rainfall intensity)，而土層 底部(GH)假設為不透水之土層，則在 z 方向流量邊界條件同樣地以 Darcy 定律加以描 述如下： CD、DF：
z z K z h K h K I Ri r (2.36)GH：
z z K z h K h K 0 (2.37) 式(2.36)以差分式可分別表示為： CD： cos 1 2 sin 2 cos 1 , 1 1 , 1 , 1 , 1 , 1 1 , 1 , 1 , 1 , 1 , 1 1 , 1 , 1 , 1 , 1 1 , 1 , z K K x K K I n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i r (2.38) DF： 1 2 1 , 1 1 , 1 , 1 , 1 , 1 1 , 1 , 1 , z K K I n m N i n m N i n m N i n m N i r (2.39)式中，Ir為降雨強度，為河岸坡面角度(angle of riverbank)。在描述降雨入滲行為時， 先假設降雨可完全入滲通過地表進行求解後，再判斷地表各計算點壓力水頭是否大於 零，若否，則表示土壤未達積水(ponding)，降雨可完全入滲，可進行下一時階(time step) 之計算。若降雨持續發生使地表計算格點壓力水頭大於 0，此時地表將產生積水，在 不考慮積水深度的假設下，亦即令im,j1,n1 0為邊界條件並於同一時階內重新計算。 考量上述的降雨入滲計算程序，將式(2.38)與式(2.39)重新排列後，分別可得： CD： 0 0 0 cos cos 1 2 sin 2 cos 1 2 sin 1 2 1 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 1 , 1 , 1 1 , 1 , 1 , 1 , 1 , 1 1 , 1 , 1 1 , 1 , 1 , 1 1 , 1 , 1 , 1 , 1 , 1 1 , 1 , 1 , 1 , n m j i n m j i n m j i r n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i n m j i I z K K x K K z K K x K K (2.40)
DF： 0 0 0 1 2 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 1 , 1 , 1 , 1 , 1 1 , 1 , 1 , n m N i n m N i n m N i n m N i n m N i r n m N i n m N i z K K I (2.41) 土層底部(GH)為零流量邊界條件，式(2.37)以差分式可表示為： GH： 1 2 0 1 , 1 1 , 1 , 1 2 , 1 , 1 1 , 1 , 1 2 , z K Kim n im n im n im n (2.42) 上式經整理後可得： GH： im,11,n1 im,21,n1z (2.43) 求解步驟 本研究以壓力水頭為唯一求解之變數，因此上述式中有關含水量以及水力傳導係 數，皆可藉由已求解之壓力水頭代入式(2.13)至式(2.16)分別求得。另外，比含水量由 式(2.19)定義其為含水量對壓力水頭之一階導數，因此，將所採用的土壤保水曲線(式 (2.13)至式(2.15))對壓力水頭一階微分後即可得比含水量與壓力水頭之關係如下：
0 0 0 1 1 v v v v m m n n v v v r s m n d d C (2.44) 圖 2.3 顯示土壤保水曲線其含水量隨壓力水頭上升而增加，因此，土壤保水曲線之斜 率(d / d )應為正值，亦即比含水量必大於 0，但由於式(2.14)中係取壓力水頭之絕對值，因此須將式(2.44)修正為
0 0 0 1 1 v v v v m m n n v v v r s m n d d C (2.45) 將上述邊界條件式(2.31)、(2.32)或(2.35)、(2.40)、(2.41)、(2.43)與(2.29)合併整理 後，可得每個計算格點上展開後之線性聯立方程組，再利用 generalized Thomas algorithm (Fletcher and Srinivas, 1991)求解。另為避免數值不易收斂的情形，在同一時 階內，在非邊界之計算點以下式進行更新並代入下一次的迭代進行計算：
1, 1 , 1 , , 1 , 1 , 1 m n j i n m j i n m j i w w (2.46) 式中，w 為權重因子，本研究預設其值為 0.5。另外，設定數值計算收斂條件( )如下： j i all for n m j i n m j i 10 , 5 1 , , 1 , 1 , (2.47) 當達收斂條件後即可進行下一時階之運算。 2.1.4 模式驗證 以下分別以一維固定水頭邊界以及二維流量邊界兩個案例進行驗證。案例 1 為 Celia et al. (1990)在一土層厚度為 100 cm，貣始條件為-1,000 cm，頂部邊界與底部邊界 分別給予-75 cm 與-1,000 cm 之固定水頭。未飽和土壤水力相關參數設定如下：Ks =0.922 m/hr，
r=0.102，
s=0.368，
v=0.0335，nv=2。數值計算空間間距z設為 1.0 m，時間間距 t 為 0.2 hr。垂向隨時間變化的壓力水頭模擬結果如圖 2.5 所示，圖中顯 示所模擬結果與 Celia et al. (1990)所得相吻合。案例 2 採用 Vauclin et al. (1979)的詴驗資料進行驗證，其研究為一寬高為 6m×2m 的土體，貣始地下水位高設定為 0.65 m，中央寬 1.0 m 處給予 0.148 m/hr 均勻入滲量， 而其餘邊界皆為零流量邊界，未飽和土壤水力相關參數設定如下：Ks=0.35 m/hr， r
=0.01，
s=0.3，
v=3.3，nv=4.1。由於該案例具對稱性，在此僅針對計算域右半部 進行模擬，模擬範圍寬高則設定為 3m×2m；空間間距 x 、z均為 0.1 m；時間間距 t 為 0.2 hr。模擬結果如圖 2.6 所示，圖 2.6 顯示與實驗數據相比對，可獲致合理之模擬 結果。2.2 河岸退縮機制
2.2.1 塊體破壞塊體破壞(mass failure)分析理論係基於極限平衡法(limited equilibrium analysis)，以 莫爾-庫倫破壞準則(Mohr-Coulomb failure criterion)定義破壞面上的正向應力與剪力強 度之關係，並利用安全係數(FS, factor of safety)做為判斷破壞發生與否之依據。FS 定 義為破壞面上抗剪力(FR, resisting force)與驅動力(FD, driving force)之比值，如下表 示： FD FR FS / (2.48) 當破壞面上之抗剪力小於驅動力或驅動力大於抗剪力，河岸邊坡將沿著已知破壞 面產生滑動破壞。根據FS 之定義，當FS 1表示河岸呈穩定狀態；FS 1表示河岸發 生破壞；FS 1則表示河岸處於臨界狀態。 土壤剪力強度
Terzaghi (1948)利用莫爾-庫倫破壞理論(extended Mohr-Coulomb failure theory)提 出飽和土壤中有效應力公式如下：
sat c uw tan (2.49)式中，
sat為飽和土壤之剪應力；c為土壤有效凝聚力(effective cohesion)；
為正向應 力(total normal stress)；uw為孔隙水壓(pore water pressure)；
uw
為有效正向應力 (effective normal stress)；為土壤有效摩擦角(effective friction angle)。Fredlund et al. (1978)則基於前述的莫爾-庫倫破壞理論，提出廣義莫爾-庫倫破壞理 論(extended Mohr-Coulomb failure theory)，將未飽和土壤之剪應力以淨正向應力(net normal stress)與基質吸力(matric suction)所組成之函數來表示：
b w a a uns c u u u tan tan (2.50)式中，
uns為未飽和土壤之剪應力；ua為孔隙氣壓(pore air pressure)；b為基質吸力造 成剪應力增加所對應之角度(angle expressing the strength increase rate relating to the matric suction)；
ua
為淨正向應力；
ua uw
為基質吸力。剪應力、淨正向應力與 基質吸力的破壞包絡線可用三維的曲面如圖 2.7 表示。當土壤達飽和時，其孔隙氣壓將趨近於孔隙水壓，因此式(2.50)中
ua uw
為 0， 則式(2.50)將會等同於式(2.49)，故式(2.50)同樣可用以描述飽和土壤的剪應力。依據該 式，在飽和與非飽和區域的剪力可分別表示為：Tsat Lsat
c
uw
tan
(2.51)
b
w a a uns uns L c u u u T tan tan (2.52) 式中，Tsat為沿破壞面之飽和土層剪力；Lsat為沿破壞面之飽和土層長度；Tuns為沿破壞 面之未飽和土層剪力；Luns為沿破壞面之未飽和土層長度。孔隙氣壓一般視為大氣壓力(ua 0)，將式(2.51)與式(2.52)合併後，則破壞面上之總剪力可表示為：
b uns w sat wL u L u L L c T tan tan (2.53) 式中，T 為總剪力，亦即為式(2.48)中破壞面之抗剪力。另外，uwLsat必為正值，而uwLuns 必為負值，則上式可進一步簡化表示為：
b S U N L c T FR tan tan (2.54)式中，L 為沿破壞面之土層總長度(LLsat Luns)；N 為總正向力(total normal force)；
U 為破壞面上因土體飽和孔隙壓力所引發之上舉力(hydrostatic uplift force)；S 為破壞面
上因土體未飽和孔隙壓力所引發之基質吸力(suction force)。另外，上述之總正向力、 上舉力與基質吸力均取正值，其計算方法分別敘述於下。 總正向力包含破壞土體重與河川水位的靜水壓力分別在破壞面法向分量之和，可 表示如下： ) cos( cos W P N (2.55)
式中，W 為破壞土體重(weight of failure block)；P 為河川水位靜水壓力(hydrostatic confining force)；為河岸坡面角度；β 為破壞面角度(angle of failure plane)。而上舉力 與基質吸力則可依據破壞面上孔隙水壓之正負值，分別計算如下： 0 0
w L wdL u u U sa t (2.56) 0 0 L
uwdL uw S u n s (2.57)式中，L 表示延破壞面方向。上式中，孔隙水壓可藉由 2.1 節所述內容求得。如圖 2.8 所示，沿著破壞面上任一點壓力水頭可由鄰近計算格點之壓力水頭值，依據距離倒數 內插計算可得，其計算公式分別如下所示：
w w u (2.58)
i i i i i l l 1 (2.59)式中，
w為水單位重(unit weight of water)；
i與li分別表示鄰近計算點之壓力水頭以 及與待求點之距離。另外，式(2.48)中，沿破壞面之驅動力為破壞土體重與河川水位的 靜水壓力分別沿破壞面方向分量之和： ) sin( sin W P FD (2.60) FS 之計算結合式(2.54)、式(2.55)與式(2.60)，則FS可表示為(Casagli et al. 1999；Rinaldi and Casagli, 1999)：
