國
立
交
通
大
學
土木工程學系
碩士論文
一維軟弱岩盤河道侵蝕數值模式之研發與應用
Development and Application of 1-D Soft Bedrock Incision
Numerical Model
學生：劉柏傑
指導教授：葉克家 博士
軟弱岩盤河道侵蝕數值模式之研發與應用
學生：劉柏傑 指導教授：葉克家
國立交通大學土木工程學系
摘要
國內目前對於河道沖蝕問題之探討，多侷限於沖積層河道，缺乏 對岩石河床、岩石河岸之數模研究。當岩床或岩岸屬地質年代較年輕 之岩層或弱面較發達之岩體，河道沖蝕問題往往甚為嚴重。國外對於 岩床之沖蝕機制已累積不少成果，但實際將岩床沖蝕機制與河道動床 數模整合者鮮少，目前常見一維動床數模如 HEC-6、GSTARS、 NETSTARS、CCHE1D 等，皆尚未有軟岩侵蝕模擬之功能。 本研究藉由現有軟岩河道侵蝕模型，建構軟岩河道沖蝕模組，包 括水力沖蝕與泥砂磨蝕兩種機制，由侵蝕模組計算獲得軟岩河床之侵 蝕率與侵蝕量，結合 EFA 動床模式，建構具有可模擬軟岩侵蝕特性 之一維動床數模。 本研究以假設案例測試軟岩侵蝕模組是否符合真實之物理特性， 同時選取台灣大安溪為應用對象，以不同年份斷面資料進行檢定驗證， 結果顯示 EFA 動床模式在結合軟岩河道侵蝕模組後，對於具有軟岩 河床侵蝕特性之河川，可有良好河道變遷趨勢上之模擬預測。關鍵字：軟岩侵蝕，有限解析法、數值模式、大安溪
Development and Application of 1-D Soft Bedrock
Incision Numerical Model
Student: Po-Chieh Liu Advisor: Keh-Chia Yeh
Institute of Civil Engineering
National Chiao Tung University
Abstract
In the present numerical study on channel migration, discussion
and handling with incision of river course is mostly limited to
alluvium rivers, lacking research for rock bed and banks. River
incision is severe when rock bed and banks are young rock stratums
or well developed joints. There are many achievements with erosion
mechanisms of rock beds, but integration of erosion mechanisms
with numerical models is few. Commonly seen numerical models
like HEC-6, GSTARS, NETSTARS, and CCHE1D, do not have the
capability to simulate soft bedrock erosion.
This study relies on the development of soft bedrock river
incision module, which include hydraulic scouring and sediment
abrasion mechanisms. Soft bedrock incision module calculates soft
bedrock incision rate and capacity. Combing the incision
mechanisms with the EFA1D, a one dimensional numerical is
developed to simulate soft bedrock incision phenomenon.
This study first tests the physical factuality of soft bedrock
incision model with fictitious cases. Furthermore, Ta-An River in
Taiwan is selected as the study reach, and the cross sections
measured in different years are used for the model calibration and
validation. Simulation results show that river migration tendency can
be well predicted after combining the soft bedrock incision module.
Keywords: Soft bedrock incision, EFA, numerical model, Ta-An
River
謝誌
本文承蒙吾師葉克家教授於研究期間不厭其煩地細心指導與諄 諄教誨，得以順利完成，在此致上最誠摯的謝意。亦感謝口試委員國 立成功大學蔡長泰教授、聯合大學理工學院院長許銘熙教授及國家高 速網路與計算中心蔡惠峰博士之細心匡正與建議，使得本論文德以更 趨完整。 感謝仲達、昇學、仁凱、詩廷、佑民學長以及曉萍學姊在研究與 生活上之指導與照顧；感謝好友冠曄、宇翔、歆淳、振家、銓謚、俊 宏的砥礪與扶持；感謝學弟們明儒、新詠、浚騰、靜宜的幫助與關心。 因為有大家的參與，使我的研究生活多采多姿，永生難忘。 最後，僅將此拙作獻給我親愛的家人，感謝你們不斷地支持與鼓 勵，使我求學的路上能無後顧之憂，才能成功的跨出這一步，謝謝你 們。目錄
摘要 ...I Abstract ... II 誌謝 ... III 目錄 ... IV 表目錄 ... VI 圖目錄 ... VII 第一章 緒論 ... 1 1-1 研究動機 ... 1 1-2 研究目的 ... 2 1-3 研究方法與流程 ... 2 第二章 文獻回顧 ... 4 2-1 岩床河道定義及特性 ... 5 2-2 岩床河道侵蝕模式 ... 7 2-3 岩床河道侵蝕速率之預測 ... 9 2-3-1 顆粒沖擊效應 (Foley, 1980) ... 9 2-3-2 剪應力效應 (Howard ad Kerby, 1983) ... 102-3-3 剪應力效應 (Seidl and Dietrich, 1992) ... 11
2-3-4 泥砂磨蝕效應 (Sklar and Dietrich, 1998) ... 12
2-3-5 水力沖蝕效應 (Annandale, 1995) ... 14 2-3-6 外在環境變動效應 (程紹平等 2004) ... 15 2-3-7 瀑布侵蝕效應 (Yuichi et al. 2007) ... 16 2-3-8 文獻回顧總結 ... 17 2-4 軟岩河道沖淤數模簡介 ... 18 第三章 模式理論基礎 ... 20 3-1 顯式有限解析法模式理論 ... 20 3-1-1 水理控制方程式 ... 21 3-1-2 輸砂控制方程式 ... 22 3-2 顯式有限解析法模式數值方法 ... 24 3-2-1 水理方程式之數值方法 ... 24 3-2-2 輸砂方程式之數值方法 ... 26 3-3 軟岩侵蝕模組模式理論 ... 30 3-3-1 水力沖蝕模式理論 ... 31 3-3-2 泥砂磨蝕模式理論 ... 34 3-4 軟岩侵蝕模組數值方法 ... 36 第四章 模式檢定與敏感度分析 ... 38 4-1 實驗室案例之檢定驗證 ... 38
4-2 軟岩侵蝕案例說明與參數設定 ... 42 4-3 模擬結果與分析 ... 43 4-4 軟岩侵蝕模組敏感度測試與分析 ... 45 4-4-1 敏感度測試條件 ... 45 4-4-2 敏感度測試分析 ... 46 第五章 現場案例模擬分析 ... 48 5-1 現場案例簡介與參數設定 ... 48 5-1-1 流域概述 ... 48 5-1-2 模擬範圍 ... 49 5-1-3 模擬所需參數 ... 50 5-2 模擬結果與分析 ... 55 5-2-1 檢定案例 ... 55 5-2-2 驗證案例 ... 56 5-2-3 軟岩侵蝕模組分析 ... 57 第六章 結論與建議 ... 59 6-1 結論 ... 59 6-2 建議 ... 60 參考文獻 ... 62 附錄 運動波-地貌瞬時單位歷線模式簡介與應用 ... 111
表目錄
表2- 1 岩石材料強度表(摘錄自 Annandale, 1995) ... 66
表2- 2 參考場址模式輸入及輸出值(摘錄自 Sklar and Dietrich, 2004) ... 67
表4- 1 各案例上游入流量與下游水位條件 ... 72 表4- 2 敏感度測試最大沖刷深度結果 ... 73 表5- 1 大安溪大斷面資料彙整 ... 74 表5- 2 模擬範圍之斷面樁號與曼寧 n 值 ... 75 表5- 3 大安溪各斷面河床質粒徑分組重量累積百分比 ... 76 表5- 4 沖蝕指數 Kh參數表 ... 77
圖目錄
圖2- 2 Bieniawski 岩石強度分類法 (摘錄自 Bieniawski, 1984)... 78
圖2- 3 河道縱向演變的兩種型態(摘錄自 Seidl and Dietroch, 1992)... 78
圖2- 4 岩床侵蝕的磨蝕機制示意圖(摘錄自 Whipple et al., 2000) ... 79
圖2- 5 岩床侵蝕的抽離機制示意圖(摘錄自 Whipple et al., 2000) ... 79
圖2- 6 Foley 河床載侵蝕模型概念圖(摘錄自 Foley, 1980) ... 80
圖2- 7 岩床侵蝕率與無床質供應量關係，粗實現為參考場址狀況 (摘錄自 Sklar and Dietrich, 2004) ... 80
圖2- 8 岩床侵蝕率與無因次相對剪應力關係，粗實線為參考場址狀況 (摘錄自 Sklar and Dietrich, 2004) ... 81
圖2- 9 岩床侵蝕率與河床質粒徑關係，粗實線為參考場址狀況 (摘錄自 Sklar and Dietrich, 2004) ... 81
圖2- 10 可侵蝕程度臨界線(摘錄自 Annandale, 1995) ... 82 圖3- 1 泥砂傳輸之示意圖 ... 82 圖3- 2 懸浮泥砂之移流特性軌跡 ... 83 圖4- 1 模型佈置圖(摘錄自 Suryanarayana 1969) ... 84 圖4- 2 非均質粒徑分布曲線圖(摘錄自 Suryanarayana 1969) ... 85 圖4- 3 Run21 淤積案例(t=1hr) ... 86 圖4- 4 Run21 淤積案例(t=2.5hr) ... 86 圖4- 5 Run21 淤積案例(t=4.5hr) ... 87 圖4- 6 Run21 淤積案例(t=7hr) ... 87 圖4- 7 Run21 淤積案例(t=10hr) ... 88 圖4- 8 Run22 沖刷案例(t=1.5hr) ... 88 圖4- 9 Run22 沖刷案例(t=4hr) ... 89 圖4- 10 Run22 沖刷案例(t=10hr) ... 89 圖4- 11 Run25 淤積案例(t=1.5hr) ... 90 圖4- 12 Run25 淤積案例(t=3hr) ... 90 圖4- 13 Run25 淤積案例(t=5hr) ... 91 圖4- 14 Run25 淤積案例(t=7hr) ... 91 圖4- 15 Run25 淤積案例(t=9hr) ... 92 圖4- 16 Run26 沖刷案例(t=1.2hr) ... 92 圖4- 17 Run26 沖刷案例(t=3.25hr) ... 93 圖4- 18 Run26 沖刷案例(t=6hr) ... 93 圖4- 19 Run26 沖刷案例(t=10hr) ... 94 圖4- 20 Run14 沖刷案例(t=1.9hr) ... 94 圖4- 21 Run14 沖刷案例(t=6.67hr) ... 95 圖4- 22 Run14 沖刷案例(t=12.67hr) ... 95
圖4- 23 測試案例渠道底床高程沿程變化 ... 96 圖4- 24 測試案例渠道水力沖蝕模擬結果高程沿程變化 ... 96 圖4- 25 測試案例渠道泥砂磨蝕模擬結果高程延程變化 ... 97 圖4- 26 測試案例渠道水力沖蝕各模擬時程底床沖淤變化圖 ... 97 圖4- 27 測試案例渠道泥砂磨蝕各模擬時程底床沖淤變化圖 ... 98 圖4- 28 測試案例渠道沖蝕指數 Kh變動模擬底床沖淤變化圖 ... 98 圖4- 29 測試案例渠道粗糙高度 ks變動模擬底床沖淤變化圖 ... 99 圖4- 30 測試案例渠道軟岩彈性模數 Y 變動模擬底床沖淤變化圖 ... 99 圖4- 31 測試案例渠道軟岩張力強度σT變動模擬底床沖淤變化圖 ... 100 圖4- 32 敏感度分析結果圖 ... 100 圖5- 1 大安溪劇烈沖蝕河段照片 ... 101 圖5- 2 模擬範圍流域示意圖(含各流量水位測站位置) ... 102 圖5- 3 流量-含砂量濃度率定曲線 ... 102 圖5- 4 斷面 44 各分區岩性示意圖 (摘自軟弱岩床具列沖蝕河段沖蝕行為之探討─以大安溪為例) ... 103 圖5- 5 斷面 44.1 各分區岩性示意圖 (摘自軟弱岩床具列沖蝕河段沖蝕行為之探討─以大安溪為例) ... 103 圖5- 6 斷面 45 各分區岩性示意圖 (摘自軟弱岩床具列沖蝕河段沖蝕行為之探討─以大安溪為例) ... 104 圖5- 7 現場檢定案例上游流量歷線(雙崎站)... 104 圖5- 8 現場檢定案例下游水位歷線(義里橋)... 105 圖5- 9 現場檢定案例底床高程沿程變化(斷面 25~42) ... 105 圖5- 10 現場檢定案例底床高程沿程變化(斷面 43~56) ... 106 圖5- 11 現場檢定案例實測值與模擬值之底床高程差比較 ... 106 圖5- 12 現場驗證案例上游流量歷線(雙崎站)... 107 圖5- 13 現場驗證案例下游水位歷線(義里橋)... 107 圖5- 14 現場驗證案例底床高程沿程變化(斷面 25~42) ... 108 圖5- 15 現場檢定案例底床高程沿程變化(斷面 43~56) ... 108 圖5- 16 現場驗證案例實測值與模擬值之底床高程差比較 ... 109 圖5- 17 現場檢定案例軟岩模組開啟與否比較圖 ... 109 圖5- 18 現場驗證案例軟岩模組開啟與否比較圖 ... 110
符號說明
A ：通水斷面積 B ：渠道寬 Ck ：粒徑k 之懸浮質濃度 D ：代表粒徑 E1 ：軟岩河床水力沖蝕速率 E2 ：軟岩磨蝕速率 Em ：作用層厚度 G ：重力加速度 H ：水深 Js ：地盤構造條件 Kb ：顆粒/塊體尺寸 Kd ：弱面/或顆粒間之抗剪強度 Kh ：沖蝕指數 Kp ：無因次沖蝕係數 ks ：相對粗糙高度 kv ：軟岩強度參數 Ls ：泥砂顆粒躍動長度 Ms ：為材料強度 n ：曼寧糙度係數 p ：孔隙率 Pcrit ：臨界流功 Q ：流量 Qbk ：粒徑k 之河床載通量 qs ：懸浮載單寬通量ql ：單位渠長之支流測流量 R ：水力半徑 S ：砂比重 Sf ：磨擦坡降 Sk ：粒徑k 之懸浮載源項 Sak ：粒徑k 於作用層底部源項 T ：時間 u ：水體速度 u* ：河床剪力速度 ws ：泥砂顆粒之沖擊速度 x ：主流方向之距離 Y ：軟岩彈性模數 Z ：水位 Zb ：底床高程 βw ：動量較正係數 βk ：粒徑k 之百分組成 ρ ：流體密度 σT ：軟岩張力強度 τ ：作用於軟岩河床之剪應力 τc ：軟岩河床之臨界剪應力
第一章 緒論
1-1 研究動機
由於自然環境變遷或人為因素影響，台灣一些主要河川已有部分 河段出現嚴重河床侵蝕下切現象。促使河道侵蝕行為加劇因素包括自 然環境變化及人為開發利用影響，例如車籠埔斷層抬升使大安溪、大 甲溪等部分河道產生嚴重之侵蝕下切現象，八掌溪、濁水溪、頭前溪 等部分河道上因水利設施、跨河構造物等，亦形成河道嚴重侵蝕下切 現象。河床面快速降低與河谷地形持續變化對於河川治理有相當大之 衝擊，諸如堤防及橋墩基礎裸露破壞、取水工程功能之喪失等工程上 之危機，均是未來所必須面對的問題。 國內目前對於河道侵蝕問題之探討，多侷限於沖積層河道，缺乏 對岩石河床、岩石河岸之數模研究。當岩床或岩岸屬地質年代較年輕 之岩層或弱面較發達之岩體，河道侵蝕問題往往甚為嚴重。近年來， 國內對於河床侵蝕或下切造成工程上的不良影響都有所體認，亦加以 探討及處理，但未對岩床侵蝕機制進行詳細探討及評估其影響。因此 本研究以探討軟岩侵蝕機制之特性，與侵蝕機制對河道變遷之影響為 研究主軸。1-2 研究目的
軟岩侵蝕所造成河川水流特性劇變致使河工構造物無法安全使 用，甚至損壞等河防安全問題，仍需透過相關模型試驗分析以及數值 模式模擬作為預測與防範之用，但模型試驗分析往往需耗費大量時間 、人力與經費，而數值模式則不需大量時間與人力即能快速提供相近 於模型試驗分析之結果，因此本研究以研發軟岩侵蝕數值模式為主要 研究方向。 本研究研發一可模擬軟岩侵蝕之數值模式，結合一維顯式有限解 析法(explicit finite analytic method, EFA)動床模式，使其具有可同時模 擬一般沖積層河道及軟岩河道侵蝕之功能，此外，透過敏感度分析探 討不同軟岩侵蝕機制參數之特性，以利於未來模式之模擬應用。1-3 研究方法與流程
本研究引用 Annandale (2006)提出之水力沖蝕機制以及 Sklar and Dietrich (2004)提出之泥砂磨蝕機制，建構軟岩河道侵蝕模組，由軟 岩河道侵蝕模組計算可獲得軟岩河床之侵蝕率與侵蝕量，結合 EFA 一維顯式有限解析法動床數值模式，建置具有可模擬軟岩侵蝕特性之 動床數值模式。另外，再以一假設實驗渠道進行模擬，探討軟岩侵蝕 模組中不同參數之敏感度特性及驗證物理現象之合理性，最後則以大
第二章 文獻回顧
本研究主要為探討於河道中軟岩(soft bedrock)對於河道變遷之影 響，依據國際岩石力學學會(ISRM)對軟岩之定義，軟岩為單壓強度介 於0.5~25 MPa 之岩體，國內泥岩、砂岩、頁岩、片岩、板岩大都可 歸類於軟岩。軟岩除強度較低之特性外，尚有膠結不良、高孔隙率、 變形性大、具潛變現象與異相性等特性。前述軟岩特性僅於文字上敘 述，在定義上較為模糊，由Bieniawski (1984)所提出依材料強度分類 法，如圖2-1 所示，則由簡單的岩石強度判斷分類則更將能明確定義。 但各學者由岩石強度提出之軟岩標準並不相同，所以本研究定義之軟 岩除參考此岩石強為標準外，另以Annandale (1995)文獻中所提出之 材料強參數為參考依據，如表2-1 所示，並輔以現地觀測等試驗，作 為岩石分類依據。 探討岩石侵蝕相關文獻最早可追溯至地質學之父 James Hutton 提出之岩石循環理論，他認為水為塑造地表大部分區域地形變化之主 要 影 響 因 素 ， 藉 由 侵 蝕 與 堆 積 行 為 改 變 原 有 之 地 貌 ， 此 概 念 被 Montgomery 延伸應用於河川發展過程之研究(Howard et al., 1994； Seidl and Dietrich, 1992)。Montgomery et al. (1993)認為河川之發展過程為複雜之水土互動 平衡結果，水藉由位能轉換為動能產生流動並消耗於河道，造成河床
泥砂等材料沉降或運移，此變化是河道產生變化最主要之動力來源。 Montgomery et al.認為河道局部輸砂能力 qc (local transport capacity)與 來自上游及河岸之輸砂量qs (bedload sediment supply)，影響河道之沖 淤特性分佈，岩床河段為區域輸砂能力超過輸砂量，而沖積河段則顯 示輸砂量高於區域輸砂能力或兩者平衡之狀態。Howard et al. (1994) 以Montgomery et al.之研究為基礎，對河道型態、產生原因及空間上 之分佈特性詳細探討，並將自然河道區分為細粒沖積河床、粗粒沖積 河床、岩石河床及沖積-岩石混合河床等。 Gilbert (1877)首先提出岩石下切速率之假說，他認為最主要影響 因素包括岩石抵抗能力、河川坡降、流量及泥砂供應量等，此假說之 概念被Howard 等學者應用發展於推估岩石河床速率之模式，其後學 者根據 Gilbert 假說之概念與前述岩石河床定義與特性相結合，陸續 發展出岩床河道侵蝕速率模式。關於岩床河道定義、特性與相關岩床 河道侵蝕模式分述如下。
2-1 岩床河道定義及特性
廣義上來說，岩床可以是任何具有凝聚性(cohesive)及阻抗性 (resistant)之材質，例如膠結沖積物或卵粒石、或是第三紀及第四紀黏 土，這些材料行為特性均可以類比於岩石。Tinkler and Wohl (1998)曾對岩床河道定義如下：整個河段沿線大 多均有岩床出露，或僅有極薄沖積層覆蓋，在高流量時此沖積層絕大 部分均會移動，位於底部岩床之幾何形狀對於水流特性及輸砂行為佔 有最重要之影響。 岩床河道許多特性均與沖積河床或礫石河床不同，Tinkler and Wohl (1998)提出以下五點岩床河道所具有之特性，包括： 1. 坡降 (gradient)：岩床河道之坡降常大於沖積河道之坡降；只 有在岩石之傾角正好接近水平時，才可能在局部區域有坡降 較平緩之狀況。 2. 變化 (change)：地形改變為單一方向性，岩石於岩床中移除 後，將降低該點以上之侵蝕基準面(base level)。 3. 變異性 (variability)：沿著岩床河道觀察到之水力特性變化， 通常可以反映底床材料在岩性上或構造上之變化。由於岩石 變化通常都很緩慢，岩床系統可以累積變化，岩石外觀即為 長時間之累積變化下之結果。 4. 阻抗能力 (resistance)：岩石對水流之阻抗能力隨著岩性不同 而改變，在岩石構造不發達時，岩性為主要影響因素；在構 造發達之水平岩層與厚層岩石，節理面等位態將影響水流對
岩石之侵蝕行為。
5. 岩性 (lithology)：岩性差異將造成不同河道外觀，如層面發達 之水平岩層與厚層岩石，河道剖面形狀與河道型態均有明顯 之差異。
Turowski et al. (2008)認為 Tinkler and Wohl (1998)等提出河道之 分類方式，應用於河道研究時，不論藉由現地觀測或是航拍圖觀察均 無法有效區分岩床河道與沖積河道。因此提出以河川斷面變化特性為 判斷基礎，作為分類岩床河道之方式，其定義如下：若河道斷面無法 以侵蝕岩床區域之外的方式增加河道斷面寬度、深度或是深槽變化等 河道變遷特性，則稱之為岩床河道。
此分類方式較先前 Tinkler and Wohl (1998)等提出之分類方式簡 便，更易應用於現地觀察區分河道，在對岩床河道越趨複雜及深入之 研究下，較能釐清與沖積河道不同之特性。
2-2 岩床河道侵蝕模式
岩床下切行為提供地殼構造及地形演變間之鏈節，岩床下切速率 決定地形演變之速率。Hancock et al. (1998)由現地河道對於硬岩(hard rock) 的 觀 察 ， 歸 納 主 要 之 侵 蝕 行 為 包 括 磨 蝕 (abrasion) 、 採 石 (quarrying)，而穴蝕(capitation)則需在異常之高流速下才可能產生。
Seidl and Dietrich (1992)以假設兩種侵蝕現象的示意圖來說明岩 床河道在縱向上的變化，圖2- 2(上)是由遷急點逐漸向上游傳遞所造 成，圖2- 2(下)不同於遷急點的傳遞，主要是由垂直於河道方向的侵 蝕行為所造成，如磨蝕作用。
Whipple et al. (2000)將侵蝕行為區分為磨蝕(abrasion)、抽離 (plucking)、及穴蝕(cavitation)等三種，並以示意圖說明磨蝕、抽離機 制。圖2- 3 為磨蝕機制示意圖，凸起岩床面前緣受到床質顆粒沖擊而 發生侵蝕，水流受到岩床面凸起的擾動，於其後產生紊流，造成壺穴 化(potholing)的侵蝕。 圖2- 4 為抽離機制示意圖，經由現地觀察及考量水流輸砂行為， 要使部分節理切割岩塊鬆動破裂並產生抽離機制侵蝕，至少包含四種 可能相關機制，其分述如下： 1. 沿著節理面進行物理性或化學性風化。 2. 水力楔型 (hydraulic wedging)，砂、卵岩或礫石經由水流擠入 楔型開口內，而使得開口破裂並逐漸延伸。 3. 大型塊石由水流帶動沖擊岩石面，瞬時之高動量傳遞造成不 均勻之應力作用，而使得垂直向及水平向之裂隙逐漸地延伸 與發展。
4. 由紊流造成之瞬時水壓擾動，造成裂隙發展。
2-3 岩床河道侵蝕速率之預測
2-3-1 顆粒沖擊效應 (Foley, 1980)
Bitter (1963a、1963b)分析砂粒沖擊表面所造成的磨蝕效應，區分 為低角度沖擊之切割行為，以及高角度沖擊之疲勞性破裂行為。Foley 將此研究應用於小區域河床面之河床載侵蝕模型(bed-load abrasion model)，其概念如圖 2- 5 所示，於單位河寬河床載傳輸率作用下，由 顆粒沖擊所造成之總侵蝕率Yt可表示成： c d t Y Y Y (2-1) 其中，由高角度沖擊造成的變形磨損侵蝕率Yd為
2 2 1 g v K Y s d (2-2) 由切割磨損所造成的侵蝕率Yc為
12 2 2 1 2 2 K U K C g Y s c (2-3) 圖 2-5 及上二式中，d 為河床載顆粒之粒徑，gs為單位河寬河床載傳 輸率，顆粒跳動高度為 n 倍粒徑，λ 為顆粒跳動長度，v 為垂直向速 度，U 為水平向速度，ε
與ξ 分別是由變形及切割磨損單位岩床體積所需的能量，K 為顆粒及岩床楊氏模數、顆粒密度、岩床彈性載重極 限值、顆粒及岩床普松比、及顆粒粒徑的函數，C 為公式係數。
2-3-2 剪應力效應 (Howard ad Kerby, 1983)
地形演變為相當長時間尺度之累積影響，對於人類時間尺度而言， 相當難以去做量測。自然或人為所造成之惡地地形(badland landscapes) 在缺乏植被保護及其軟弱岩性之本質下，可以是一個代表自然界大型 系統之縮影，保留了大型系統之許多特性。 基於對兩個惡地地形長達十五年之觀察，透過地形資料之統計分 析，Howard and Kerby 認為岩床河道下切速率隨剪應力而增加，所以 侵蝕速率是河道坡降S 及集水面積 AD的增函數，表示如下： n m DS KA dt dz (2-4)Howard and Kerby 以最小平方法迴歸估計常數 K 為-0.11，指數 m 為 0.11，指數 n 為 0.68。
Howard (1994) 將 前 期 之 模 式 歸 類 為 傳 輸 限 制 型 模 式 (transport-limited model)，此種模式應用於模擬集水區地形演變時，通 常會高估實際之傳輸率，透過現地觀察，植生之根莖葉、岩層或風化 土層之凝聚力均有保護地表及抗侵蝕之能力，因此對於集水區之演變
Howard 發展分離限制型模式並進行驗證，但此模式並非專注於岩床 河道之侵蝕過程，所以於此僅介紹其概念。
2-3-3 剪應力效應 (Seidl and Dietrich, 1992)
Seidl and Dietrich (1992)認為 Howard ad Kerby (1983)之模式對於 預測岩床侵蝕速率之參數中，並未考量上游之泥砂供應，所以該模式 只適用於岩床面。Seidl and Dietrich 並提出一個測試 Howard and Kerby 模式之方法，假設一支流與其主流之下切速率相等，則在會流 口處兩者應有相同之高程，假設二者之其他條件大致相同，則可得 到： n t m t n p m pS A S A (2-5a) 或 n m t p p t A A S S (2-5b)
上式中，下標p 及 t 分別代表主流及支流。Seidl and Dietrich 以美國
Coast Ranges of Oregon 集水區為應用區域，採用美國地質調查所 (USGS)之 1/24,000 地形圖，其分析結果顯示，低坡度河川之侵蝕與 水流造成之剪應力呈線性關係，此時m/n=1；但較陡河川之侵蝕則是
2-3-4 泥砂磨蝕效應 (Sklar and Dietrich, 1998)
此模式僅考慮由床質顆粒跳動(saltation)所造成之侵蝕，忽略其他 可能造成岩床侵蝕之因素或機制，如懸浮載造成之侵蝕，以及穴蝕 (cavitation)、溶解(dissolution)、岩塊抽離(plucking)等現象之影響。僅 考慮顆粒跳動所造成的侵蝕是因為此現象是水流傳遞動量最直接之 方式，且泥砂運移之行為也屬目前較能清楚理解之範疇。 岩床侵蝕速率假設與河床載顆粒撞擊岩床表面垂直方向之動量， 與超出將岩石顆粒自岩床表面移除所需動量成比例關係，其方程式表 示如下：
t s s s v t s s s Q Q W D Q v u E sin( ) 3 1 2 2 (2-6) 式中，E 為岩床侵蝕速率(m/s)，D 為均質顆粒粒徑(m)，us及 vs分別 是床質顆粒撞擊時之水平方向及垂直方向速度(m/s)，Qs輸砂量(kg/s)， Qt為輸砂能力(kg/s)，W 為河寬(m)，λ 為跳動距離(m)，α 為顆粒撞擊 之角度，ε
v 為移除單位岩床體積所需之能量(J/m3)，ε
t 為移除材料 的能量臨界值(J)。輸砂能力可採用各種輸砂經驗公式計算。Sklar and Dietrich 模式中定義岩床侵蝕速率主要為三個因子之 乘積，分別是：
1. 每一個顆粒沖擊岩床表面造成與岩床分離之材料體積，是顆 粒垂直於岩床面之動能函數，必須高於要分離材料之能量臨 界值
ε
t。 2. 每單位面積顆粒的沖擊率，與單位寬度總床質流量、粒徑尺 寸、及跳動距離相關。 3. 岩床面裸露比例，假設與輸砂量及輸砂能力有關。Sklar and Dietrich (2004)針對前期模式重新回顧，並將公式改寫 成更簡化之形式，以美國北加州之South Fork Eel River 為參考場址 進行測試，測試條件及模式輸出如表2-2 所示。圖 2-7 表示由 Sklar and Dietrich 模式預測之岩床侵蝕率與床質供應量在不同河道坡降下之關 係，岩床侵蝕率將隨著床質供應量增加而達到最大值，之後即隨著床 質供應量之增加而下降至最後停止侵蝕，此時表示岩床面已經完全被 河床質所覆蓋。另外，最大侵蝕率將隨著坡降之增加而增加，到達最 大值之後，隨著坡降增加而降低，顯示超量之輸砂能力將使得顆粒的 跳動距離大幅增加，因而降低顆粒沖擊岩床表面之頻率。 圖2-8 表示該模式預測之岩床侵蝕率與無因次相對剪應力之關係， 由圖可知於相同河床質供應量下，無因次相對剪應力須大於門檻值才 會發生侵蝕。圖2-9 則為該模式測試岩床侵蝕率與河床質粒徑之相對
關係，於相同河床質供應量下，侵蝕率隨著粒徑增加而增加，但達到 最大侵蝕率後，若粒徑再增加，沖刷率則下降。
2-3-5 水力沖蝕效應 (Annandale, 1995)
Annandale 針對沖刷行為提出之模式主要包含兩個部分，分別為 水 流 的 侵 蝕 能 量(erosive power of water) 以 及 材 料 抗 沖 蝕 能 力 (erodibility)，其關係式表示如下： ) (Kh f P (2-7 ) 其中，PqE，Kh為抗沖蝕指數，γ 為水的單位重，q 為單位寬度流 量，ΔE 為能量損失；Kh又可表示為： s a b s h M K K J K (2-8) 式中各參數可由對應之評分表獲得各參數之代表值，Ms 代表顆粒性 土壤/凝聚性土壤/岩石/風化物等材料之強度；Kb代表顆粒/塊體尺寸； Kd代表弱面/或顆粒間之抗剪強度；Js代表地盤構造條件。 Annadale 以該模式計算了 150 處不同河道之抗沖蝕指數。又依不 同河道流況如水躍、坡度突然改變、均勻明渠等，依明渠水力學理論 計算各處河道之能量消耗，並記錄現場是否發生沖蝕，將能量消耗與 抗沖蝕指數繪出雙對數關係圖，如圖 2-10 所示，進而推估對應特定
抗沖蝕指數條件下發生沖蝕時所需之能量消耗門檻值。
2-3-6 外在環境變動效應 (程紹平等, 2004)
程氏等針對河道中流水侵蝕岩盤的機制作一整合性的探討，整理 國內外文獻並歸納出三種主要流水侵蝕機制，分別為動力學機制、物 理侵蝕機制，及外在環境變動影響。動力學機制相關概念與前述 Howard and Kerby 模式概念相似，整理出一系列河川特性與侵蝕速率 之關係式；物理侵蝕機制概念與Sklar and Dietrich 所提出泥砂磨蝕效 應相似，此二種機制於此不再贅述。 外在環境變動方式又可分為三種：包含地盤構造變化、氣候改變 及侵蝕基準面變化，其內容分述如下： 1. 地盤構造變化：常見之地盤構造變化，如斷層錯動造成之地 表隆起、下陷等，造成河川坡度之增加而產生河流回春作用， 導致侵蝕基準面的下降，諸如此類等影響將改變河川之侵蝕 能力，使其侵蝕能力上升。 2. 氣候變化：此處所提氣候變化主要為冰河期與間冰期兩種極 端氣候現象所造成的侵蝕影響。於冰河期間，氣候不利於河 川流域內植物生長，而導致河床質增加，河川對岩盤侵蝕能 力降低；而間冰期間，溫濕條件利於植被生長，植被茂盛之
條件下，河川水流侵蝕能力增加。 3. 侵蝕基準面變化：侵蝕基準面為一假想之水平面，侵蝕作用 朝著侵蝕基準面進行。侵蝕基準面變化造成侵蝕效應之改變 主要取決於其變化之方向、大小、與速度。侵蝕基準面變化 方向影響著河川屬性為沖蝕或淤積型態；而侵蝕基準面變化 幅度若大，則侵蝕或淤積作用強度隨之增加，反之侵蝕淤積 作用強度減小；侵蝕基準面變化速度亦影響侵蝕作用特性， 若侵蝕基準面向下變化速度快，則河流展現垂直向下侵蝕特 性；若向下變化速度慢，則河川可能藉由側向侵蝕作用調整 河川坡度。
2-3-7 瀑布侵蝕效應 (Yuichi et al., 2007)
Yuichi et al. (2007)等認為在瀑布條件下所造成的各種流水下切 侵蝕模式中，當侵蝕速率過快時，軟岩侵蝕為主要侵蝕原因，可用包 含流水侵蝕能力與阻力(force/resistance)參數之關係式做為推估瀑布 條件下軟岩侵蝕速率之基礎，其流水侵蝕能力與阻力參數關係式如 下： c S WH AP FR (2-9)式中，A 為上游集水區面積，P 為該區域平均年降雨量，W 與 H 為瀑 布之寬度與高度，ρ 為水之密度，Sc為該區軟岩無圍單壓強度。推估 瀑布流水下切侵蝕速率經驗式如下： 73 . 0 73 . 0 99.7 7 . 99 c R S WH AP FR E (2-10) 式中，常數99.7 為一具因次之係數。Yuichi et al. 利用此經驗式估算 日本西南區 Aso 火山區內十一座瀑布，並與實測資料相互比較， Yuichi et al. 認為此經驗式可推估該區域瀑布之軟岩侵蝕平均速率。 但於火山區域中，地震等外部營力活動旺盛，各種外部營力為地形變 化之主要原因，若僅需概算瀑布區軟岩侵蝕速率則該經驗式仍可提供 一相當程度近似之估算。
2-3-8 文獻回顧總結
岩石河床侵蝕速率預測之文獻，以 Annandale (1995)之水力沖蝕 效應與Sklar and Dietrich (1998)之泥砂磨蝕效應較能量化岩石河床之 侵蝕速率及侵蝕量，且 Sklar and Dietrich 所提出之泥砂磨蝕效應與 Foley (1983)之顆粒沖擊效應相比，泥砂磨蝕效應適用性較為廣泛， 且考量因子較完整，較適於後續研究之應用。重於地表岩石侵蝕現象，如惡地地形演變等，是否適用於岩石河道侵 蝕研究仍有待考量。而Yuichi et al.提出之侵蝕模式，主要以其研究區 域觀測之數值迴歸分析而得，與前述等文獻相比，較無理論基礎，且 是否適用於其他區域仍有待考量。
因此本研究於岩床河道侵蝕理論方面，Annandale (1995)之水力 沖蝕效應與 Sklar and Dietrich (1998)之泥砂磨蝕效應為主，並結合 EFA 河道動床數值模式，作為軟岩河道侵蝕數值模式發展應用之研究 基礎。
2-4 軟岩河道沖淤數模簡介
具模擬軟岩侵蝕特性之數值模式，目前有美國墾務局 (Bureau of Reclamation, US Department of Iterior)之 SRH-1D，SRH-2D 模式，美 國 計 算 水 科 學 與 工 程 中 心 (National Center for Computational Hydroscience and Engineering)之 CCH2D 模式。為與 EFA1D 模式之功 能比較，茲簡介SRH-1D 模式如下。
SRH-1D 全 名 為 Sedimentation and River Hydraulics - One Dimension，其前身為 GSTARS 模式，為美國墾務局最新研發之一維 河道沖淤模式。SRH-1D 於水理演算方面，可進行定量流、變量流、 超亞臨界流等水理演算；動床模擬方面，具備模擬非均勻輸砂、凝聚
性沉滓、懸浮載與推移載分離運算等功能，模式亦能反映護甲層效應。 輸砂公式則有Meyer-Peter and Muller、Parker、Wilcock and Crowe 等 13 種公式可供選擇使用。SRH-1D 軟岩侵蝕模組之模式理論參考引用 Annandale (1995)之水力沖蝕理論與 Sklar and Dietrich (1998)之泥砂磨 蝕效應理論。
SRH-1D 前身為 GSTARS 模式，於數值模擬上，同樣具備流管特 性，應用於河道斷面之幾何處理上，可將河道細分為不同流管組成， 不同流管亦可設定不同底床特性，模式提供四種底床特性可供選擇： 河道河床質(channel bed material)、可沖刷河床質(erodible bed material)、 可沖刷裸露軟岩(erodible exposed bedrock)及軟岩岩盤(bedrock)四類， 利用此種方式分離一般動床輸砂模擬與軟岩侵蝕模擬。將此特性應用 於斷面變化處理上，則可使底床具側向變化模擬，使其具有擬似二維 動床模式之功能，此為SRH-1D 一大特色。 SRH-1D 模式曾應用於大甲溪與濁水溪河道之軟岩案例 (水利規 劃試驗所，2008)，但由於缺乏現地調查資料，軟岩參數之給定多採 用假設值，且僅有測試案例分析之結果。目前水利署仍有SRH-1D 相 關研究計畫，目的為尋求適當之沖刷公式納入SRH-1D 模式中，用以 預測台灣河川之沖淤行為。
第三章 模式理論基礎
3-1 顯式有限解析法模式理論
本研究之顯式有限解析法(explicit finite analytic, EFA)數值模式 可分成水理及輸砂兩大部分，於水理計算方面，在滿足de Saint Venant 之基本假設下，應用顯式有限解析法離散動量方程式，並利用沿特性 線積分概念，求解動量方程式，再配合適當斷面處理與差分形式，求 解水流連續方程式。顯式有限解析法有推導容易且程式撰寫較隱式法 簡易之優點，而此法在計算流力與水力計算領域之應用已證實有相當 不錯之成果。 至於輸砂計算方面，因天然河床係由多種不同粒徑之泥砂所組成， 且上下游河段之底床粒徑大小有所差異，因此所採用之模式須具有模 擬非均勻泥砂之特點，以反映河道泥砂部份之不同特性。此外，為考 慮懸浮載與河床載不同之運移機制，故將兩者予以分開計算，並考慮 泥砂在河道底床附近發生沉淤與再懸浮之情形；因此本模式引用許氏 (2002)對懸浮載與河床載間交換速率估算之研究，藉以推估水體中各 懸浮泥砂之濃度變化，以及河床上床質粒徑之組成。本研究於處理輸 砂控制方程式時，採取適用於雙曲線型方程式之特性法求解懸浮載質 量守恆方程式，並與河床載質量守恆方程式及整體河床輸砂之質量守
恆方程式進行結合演算，再利用Newton-Raphson 法疊代聯立求解。
3-1-1 水理控制方程式
水理演算係根據de Saint Venant 所推導之一維緩變非穩流控制方 程式，其基本假設如下： 1. 流速均勻分佈：流速均勻分佈在通水面積上，即每一個通水 斷面積僅存在一個流速，此即一維水流。 2. 靜水壓分佈：假設河道中水流之垂向流線曲率很小而且忽略 其垂直加速度，因此水深方向速度梯度為零，可忽略垂向加 速度，則假設成立。 3. 河道定量流摩擦損失估計：渠底摩擦與紊流效應對水流所造 成的損失，可以定量流摩擦律估算。 4. 底床坡度甚小：當假設成立時，重力沿河道所造成的分力將 會很小，甚至可忽略不計，亦即水深可以垂向水面與渠底高 程差表示。 5. 忽略柯氏力及風力的影響。 對於不可壓縮水流之控制方程式，包括水流連續方程式與水流動 量方程式，為如下形式。 水流連續方程式0 l q x Q t A (3-1) 水流運動方程式 0 2 l l f w gAS qu x Z gA A Q x t Q
(3-2) 式中，A 為通水斷面積；Q 為流量；t 為時間；x 為主流方向之距離； g 為重力加速度；Z 為水位；w為動量校正係數；ql為單位渠長之支 流側流量，q l為正屬合流之處理，ql為負屬分流之處理；ul為支流在 主流方向的速度分量； 3 4 2 2 R A n Q Q Sf 為摩擦坡降，其中 R 為水力半徑，n 為曼寧值。3-1-2 輸砂控制方程式
泥砂傳輸方式可依據運動機制不同而分為滾動、滑動、跳躍以及 懸浮等形態，通常可將滾動、滑動、跳躍這三種型態合稱為河床載(bed load)，或稱為推移載，一般來說，此類傳輸型態的共通點在於沈滓顆 粒明顯受到重力作用影響，並使其運移範圍僅限於河床上方有限高度， 即所謂參考高度內。然而，當紊流作用大於沈滓懸浮之啓始條件時， 部份沈滓顆粒即向上運移超越參考高度，形成懸浮載(suspended load) 傳輸型態；同樣的，懸浮泥砂在受重力的影響下亦會逐漸沈降至底床， 進而轉變為河床載傳輸型態，圖3- 1 為泥砂傳輸之示意圖。透過以上敘述可知，泥砂會因所受作用力的不同而改變其傳輸機制，因此可將 此種特性視為輸砂行為中極重要之物理機制。 EFA 模式將輸砂控制方程式中河道輸砂通量(即河床質載)分離為 非均勻之懸浮載與河床載兩部份，同時求解某一粒徑k 之懸浮載、作 用層質量守恆，及整體河床質載之質量守恆等控制方程式，分別表示 如下： ( ) ( ) ( ) k k k k C A C Q C A S t x x x k=1,2,...,TK (3-3) ( ) (1 ) k m bk 0 k ak BE Q p S S t x k=1,2,...,TK (3-4) 1 ( ) (1 ) b TK ( bk ) 0 k k BZ Q p S t x
(3-5) 上三式中，Ck 為某一代表粒徑 k 之懸浮載濃度； 為孔隙率；k 為 作用層內粒徑 k 之百分組成；Em為作用層厚度；Qbk為粒徑 k 之河床 載通量；Sk 為粒徑k 之懸浮載源項；Sak為粒徑 k 於作用層底部源項； b Z 為底床高程；B 為渠道寬；TK 為非均勻泥砂之代表粒徑數。 由上三式可知，尚有部份參數如Qbk、 、Sk及Sak須利用一些輔 助經驗公式來決定，而這些公式大都由許多學者藉由實驗及既有經驗 式修正求得。如前所述，因模式將懸浮載與河床載分開計算，故本研 究在經驗公式上以Van Rijn (1984)公式為主體，輔助相關學者及模式 之經驗加以調整並應用於模式中，詳細理論可參考許(2002)之論文。 p m E3-2 顯式有限解析法模式數值方法
由於河道水理計算之控制方程式為非線性聯立方程組，故本數值 模式沿用葉等(1996、1997)成果進行河道水理演算。水理模式採用顯 式有限解析法(EFA)進行水理控制方程式之離散化，此數值方法主要 係用以求解雙曲線型偏微分方程式，符合移流項之數學形式。輸砂模 式方面採用與水理分離演算(uncoupled)的計算方式，即在每一計算時 段內先求解水理條件，再以此推估輸砂量與底床沖淤量等，反之水理 條件受到輸砂行為的影響則在累進時間的過程中反應。3-2-1 水理方程式之數值方法
連續方程式保存保守型方程式之特性，並以控制體積觀念來差分 之，用以求得水位變化量。動量方程式則因其具有雙曲線型方程式之 特性，故針對移流項之部份採用顯式有限解析法予以處理。經離散後 之連續控制方程式如下：
0 2 1 2 1 1 1 1 1 1 1 1 x q Q Q x q Q Q t A A n li n i n i c n li n i n i c n i n i (3-6) 式中， n1 i A 為未知數，上標為(n+1)者，係先給定 n 時刻之量測值，經 反覆疊代後，再將(n+1)時刻所計算之值帶入；
c與以下的 分別m其範圍在［0,1］之間。qli為支流流量，合流時qli為負，分流時qli為 正。 經離散後之動量控制方程式如下：
0 1 1 1 1 1 1 1 1 1 1 1 n i l n i l n i f n i d n l n r n m d n l n r n i m d n l n r n m d n l n r n i m n n i u q S gA x n Z Z gA x n Z Z gA x n v v Q x n v v Q t Q Q (3-7) 式中， n1 i Q 為未知數，而結合特性線與有限解析法之觀念，依水流方 向及流況採用不同的差分方式。當流況為亞臨界流(-1<Fr<1)時，則 r =i+1，l=i-1，nd=2，代表中央差分；當流況為正（向下游）之超臨界 流(Fr＞1)，則 r =i，l =i-1，nd=1，代表後項差分；當流況為負（往上游）之超臨界流時(Fr＜-1)，r =i+1，l=i，nd=1，代表前項差分，i
為計算點位置。上標者為時間點，下標者為空間位置，
t
為時間間 距，
x
為二斷面之間距，下標符號 者為 n 時段上之特性線位置， 該特性線係由n+1 時段上之計算點向後（backward）投射，此為顯式 有限解析法求解特色之一。 當超亞臨界混合流發生時，超臨界流下游端水躍發生處，會因上、 下游水深之共軛特性自然形成一限制條件，對整個流場而言，此限制 條件即為一內部邊界條件，其存在使超臨界流區域得到多於之邊界條件，導致數值計算發生over-determined 現象。鍾(2009)利用內部相鄰 計算點之水位高程，透過預測-修正數值方法先求取超亞臨界混合流 區域內各斷面通水面積及流量，以減少地形驟變對模式演算之影響與 降低數值震盪幅度，並增加之後代入顯式有限解析法運算之穩定性。 顯式有限解析法利用已預測修正後之通水面積、流量等數值繼續運算， 得到真正所需超亞臨界混合流區域之通水面積、流量及超臨界流區域 之下游邊界水深。鍾(2009)以此方式改良修正顯式有限解析法於超亞 臨界混合流流場計算之限制，詳細理論可參考鍾(2009)之碩士論文。
3-2-2 輸砂方程式之數值方法
由於輸砂方程式中各物理量具有高度之相關性，如河床質與懸浮 質藉著懸浮載源連結，河床高程及作用層內之粒徑百分組成又因懸浮 載源而變動等，故有賴結合演算方式求解此三條方程式，所採用之方 法為半隱式法差分聯立求解。因懸浮載之質量守恆方程式依物理特性 可分割為移流及反應項(advection-reaction term)與擴散項(diffusion term)二部份，故首先將懸浮載質量守恆之移流及反應項與作用層質量 守恆方程式、整體河床質載守恆方程式，利用 Newton-Raphson 法疊 代聯立求解；然後，根據所獲得之變數值，再與懸浮載質量守恆方程 式之擴散項反覆疊代至收斂為止。各控制方程式離散後之形式如下：作用層質量守恆差分式： ] ) ( ) [( 2 ] ) ( ) [( ) 1 ( 1 1 1 1 1 n i bk n i bk n i m k n i m k Q Q x BE BE t p 0 ) )( 1 ( ) ( ] ) ( ) [( 2 ) 1 ( 1 1 1 n i a k n i a k n i bk n i bk Q S S S S Q x (3-8) 整體河床質載守恆差分式：
n k n i bk n i bk n i b n i b Q Q x BZ BZ t p 1 1 1 1 1 1 ( ) ( ) 2 1 ] ) ( ) [( ) 1 (
(1 )
0 ]} ) ( ) )[( 1 ( 1 1 1 1
n k n k n k n i bk n i bk Q S S Q (3-9) 圖3- 2 是以 一 維 空 間 為 例 示 意 沈 滓 之 移 流 軌 跡。就tn1計 算 時 刻 之 計 算 點 A 而 言 ， 移 流 軌 跡 存 在 於tn1與tn時 刻 之 間 ， 定 義tn1時 刻 之 端 點 A 為 到 達 點 (arrival point)，tn時 刻 之 端 點 D 為 離 開 點(departure point)。離 開 點 D 之 懸 浮 沈 滓 濃 度 可 由 初 始 條 件 求 得，但 由 於 模 式 採 用 固 定 格 點，故 離 開 點 D 並 不 會 剛 好 落 在 格 點 上 ， 因 此 該 點 之 濃 度 必 須 藉 由 鄰 近 格 點 濃 度 以 內 插 的 方 式 求 得 。 此 外 ， 在 一 般 情 況 下 ， 河 床 質 移 動 量 遠 小 於 懸 浮 質 移 動 量 ， 即 兩 者 在 時 間 尺 度 上 相 差 甚 大 ， 在 必 須 聯 立 求 解 的 前 提 下，懸 浮 載 方 程 式 須 使 用 較 大 之 可 蘭 數(Courant number)，此 將 導 致 移 流 軌 跡 穿 越 若 干 個 計 算 格 點 空 間 ， 因 而 必 須 採 分 段 處理 以 求 得 較 正 確 之 移 流 軌 跡。假 設 移 流 軌 跡 從 離 開 點 D 至 到 達 點 A 共 跨 越 LNS 個 計 算 格 點 空 間 ， 且 將 該 軌 跡 進 入 及 離 開 各 計 算 格 點 空 間 之 座 標 依 序 編 號 為(LNS+1)個 節 點 ， 則 各 節 點 間 的 相 對 位 置 可 以 表 示 為 ： 1 1 1 ( ) ( ) 2 l l l l l l u u x x t t l=1,2…, LNS (3-10) 其 中 ， 下 標l為 節 點 編 號 ，l=LNS+1 代 表 到 達 點 A，l=1 代 表 離 開 點 D。利 用 式 (3-10)推 求 各 節 點 位 置，必 須 要 先 知 道 各 節 點 上 的 移 流 速 度 ， 但 移 流 速 度 又 與 節 點 位 置 有 關 ， 可 利 用 疊 代 收 斂 的 方 式 來 推 求 一 正 確 的 移 流 軌 跡 。 1. 懸浮載質量守恆特性方程式： 當運動軌跡已知後，吾人即可積分懸浮載質量守恆控制方程式， 其離散化之方程式為： 1 1 1 ( ) ( ) [( ) ] 2 LNS l l k arr k dep l l l t t C A C A S S
+ C C C t x t C C C x n i k n i k n i k n i k n i k n i k [( ) 2( ) ( ) ] ) ( A ) 1 ( ] ) ( ) ( 2 ) [( ) ( A 1 1 2 1 1 1 1 1 2 x u u A C i i i 2 1 1 (3-11) 由 上 式 可 清 楚 瞭 解 到，在 水 深 平 均 模 式 中 沈 滓 交 換 速 率S
所 造 成 之 影 響 會 直 接 反 應 在 水 深 平 均 濃 度 的 改 變 上 ， 應 被 視為 一 源 項 ， 為 影 響 沈 滓 濃 度 分 佈 的 重 要 物 理 量 。 經由以上數值離散處理後，在非均勻沈滓共區分為TK 個粒徑區 間的情況下(TK 2)，任一計算點共可得(2TK 1)條代數關係式，包括 TK 條懸浮載質量守恆離散式、TK 條作用層質量守恆離散式及 1 條整 體河床輸砂質量守恆離散式。但在考慮均勻沈滓的情況下，則僅存懸 浮載與整體河床輸砂質量守恆離散式各1 條，至於作用層質量守恆離 散式則退化成 1的恆等式，符合單一粒徑時之情況。 為方便說明起見，在計算點上之未知量可以向量形式表式如下： 1 1 1 1 ( , , ,..., , ,..., , ) n TK TK k k b n z c c c s (3-12) 或更簡潔地寫成： 1 1 2 2 1 1 ( , , ) n k k n s s s s k =1,2,…,TK (3-13) 其中，k 為粒徑區間之標號。則整體河床輸砂質量守恆離散式、作用 層質量守恆離散式與式(3-11)可分別寫成： 0 ) ( 1 1 n s F (3-14) 0 ) ( 1 2 n k s F k =1,2,…,TK (3-15) 0 ) ( 1 1 2 n k s F k =1,2,…,TK (3-16) 以上三式為非線性代數式，可線性化後利用 Newton-Raphson 法疊代
求解： ) ( ] [ 1 1 1 s F msn s F (3-17) ) ( ] [ 1 2 2 m n k k s F s s F (3-18) ) ( ] [ 1 1 2 1 2 m n k k s F s s F (3-19) 式中，F s為Jacobian 係數矩陣中之列向量；msn1 為前一次疊代所 得 之 向 量 ， 上 標 為 疊 代 計 數 ； s 為 疊 代 修 正 向 量 ， 可 表 為 ) , , ( 1 2 2 1 s s s k s k 。解得修正向量s後，可得新的m1sn1向量： s s sn m n m1 1 1 (3-20) 當 s 小於某一收歛容許值時，疊代得以結束。
3-3 軟岩侵蝕模組模式理論
目前國內外於文獻上發表之河道軟岩侵蝕機制，均將複雜的交互 影響機制簡化成一個整體侵蝕速率關係式，諸如平均河床剪應力、或 單位河川流功率。根據前述之文獻回顧，以Annandale (1995)之水力 沖蝕理論與Sklar and Dietrich (2004)之泥砂磨蝕理論較能定量描述河 道中軟岩侵蝕問題。而軟岩受到水流作用之塊狀崩落或風化沖洗等其 餘侵蝕現象，目前仍未有較具明確之經驗公式或數學方程式可描述， 多為實驗或物理過程之觀測研究，數模應用上較為困難。有鑑於此，本研究挑選 Annandale (2006)之水力沖蝕機制與 Sklar and Dietrich (2004)之泥砂磨蝕機制作為軟岩侵蝕模式之理論基礎，本研究定義水 力作用於軟岩表面之沖刷為沖蝕，泥砂作用於軟岩表面之沖刷為磨蝕， 以下則針對相關理論加以說明。
3-3-1 水力沖蝕模式理論
一般來說，水流作用力常對軟岩節理處產生瞬間或逐漸沖蝕破壞， 或對單獨軟岩塊體表面產生沖擊性之沖蝕。此二種機制均屬水力沖蝕 之範疇。此外在紊流條件下，瞬間剪應力與流速將隨著河川平均流功 而改變，與穩流條件下之穩定作用力相比，此瞬間作用於岩床之剪應 力或壓力差為造成岩床破壞之主要原因(Annandale, 2006)。 Annandale (1995)曾以水流沖蝕能量為概念，探討河川流功與沖 蝕指數之相對關係。沖蝕指數之概念最早係由Kirsten (1982)所提出， 沖蝕指數可以以下形式所表示： s d b s h M K K J K (3-21) 式中，Ms為材料強度；Kb為顆粒/塊體尺寸；Kd為弱面/或顆粒間之抗 剪強度；Js 為地盤構造條件。各指數之物理特性及計算方式如下所 述： 材 料 強 度 參 數 Ms 主 要 代 表 岩 石 密 度 以 及 無 圍 單 壓 強 度 (unconfined compressive strength, UCS)等特性，無圍單壓強度可由ASTM D-7012 標準試驗獲得。Kirsten 對於計算 Ms之方式表式如下：
MPa UCS UCS C MPa UCS UCS C M r r s 10 , 10 , 78 . 0 1.05 (3-22) 式中 7 . 2 r r S C ，Sr為該材料之比重。 顆 粒/ 塊 體 尺 寸 參 數 Kb 主 要 由 岩 石 品 質 指 標(rock quality designation, RQD)與節理數量 Jn 所決定，岩石品質指標 RQD 可由 ASTM D6032 標準試驗獲得；而節理數量 Jn則依靠現地觀測後，再對 照參數表所決定，參數表如表3-1 所示。其關係式如下： n b RQD J K / (3-23)若該材料為非黏性粒狀材料(non-cohesive granular material)，則可簡化 為Kb=1000d3，d 為材料粒徑大小。 弱面/或顆粒間之抗剪強度 Kd，其參數代表不連續材料接觸面的 相互阻力大小，Kirsten 對於 Kd的計算方式如下所示： a r d J J K / (3-24) 式中，Jr與Ja分別代表材料接觸面之粗糙程度與改變參數，粗糙程度 Jr受到不連續面形狀之影響，判斷節理面是否分離後，再依接觸面不 同形狀後即可得此參數，如表3-2 所示。節理改變參數 Ja與不連續面 空隙之填充材料有關，如植生、凝聚性或非凝聚性材料填充，會影響 不連續面間之磨擦力，判斷節理分離程度後，再依空隙填充材料即可
地盤構造條件參數 Js，其代表岩石材料於地面顯露時，抵抗侵蝕 的能力，參數值與水流流向、岩床較密節理傾向、節理傾角、以及岩 塊形狀有關。岩塊形狀則利用節理間距比 (ratio of joint spacing, r)代 表岩塊形狀對抗侵蝕程度之影響，從水流過岩床之縱剖面觀察兩個方 向之岩塊長度比值y/x 作為 r 之計算，其中 y/x 最大值為 8。r 參數能 反映出岩石河床發生侵蝕時，受長塊體較等邊塊體抗侵蝕能力高之行 為。判定河流流向與岩層位態關係─岩層位態瞬向或逆向河向將影響 地盤構造條件參數Js，其相關參數如表3-4 所示。 根據 Annandale (2006)對河川流功與沖蝕指數關係之研究，沖蝕 指數 Kh之數值可將河川流功作用之方式區分為兩種類型，分別為顆 粒狀材料與岩盤材料兩種類型，其臨界流功 Pcrit (kW/m2)與沖蝕指數 Kh之關係式如下： 1 . 0 , 1 . 0 , 48 . 0 75 . 0 44 . 0 h h h h crit K K K K P (3-25) 綜合上述，軟岩河道水力沖蝕速率可藉由結合河川水流之臨界 流功、平均流速、河床剪應力等物理參數估算之。水力沖蝕速率經驗 公式型式可表示如下(Greimann, 2008)： 1 1 c pU K E (3-26) 式中，E1為軟岩河床水力沖蝕速率(m/s)；Kp為無因次沖蝕係數，該
係數需由試驗結果以及現場資料進行檢定；U 為河川斷面平均之水流 流速(m/s)；τ 為作用於軟岩河床之剪應力(N/m2)；τc為軟岩河床之臨 界剪應力(N/m2)，該參數為臨界流功之函數，其關係式如下： 3 2 6 1 66 . 7 R k Pcrit s c (3-27) 其中，ρ 為流體密度(kg/m3)；R 為水力半徑(m)；ks 為相對粗糙高度 (roughness height)。
3-3-2 泥砂磨蝕模式理論
一般於河道表面以河床載形式運移之泥砂顆粒，對軟岩河道同時 具有具侵蝕特性的磨蝕效應與覆蓋保護特性的工具效應(tool effect)。 最大軟岩磨蝕量則受到河床載之臨界供砂量(critical sediment supply) 限制；若運移之河床載超越臨界供砂量，則泥砂將落淤於軟岩河道表 面形成覆蓋，並對河道表面產生具保護特性之工具效應，使其不受磨 蝕。依據 Sklar and Dietrich (2004)之實驗結果，發現軟岩磨蝕速率與 軟岩表面張力強度成反比，但與泥砂顆粒質量成正比，並在達到最大 泥砂傳輸率後，軟岩磨蝕速率將下降。當泥砂顆粒質量較小時，磨蝕 作用將會增加更多顆粒磨蝕軟岩表面；而超過最大值時，泥砂顆粒則
Sklar and Dietrich (2004)修正 Foley (1980)理論後，提出泥砂磨蝕 速率公式如下： t s s T v s s q q L k Y w q E 2 1 2 2 (3-28) 式中，E2 為軟岩磨蝕速率(m/s)；qs 為單位河寬之供砂量(kg/m/s)；qt 為單位河寬之輸砂能力(kg/m/s)；ws 為泥砂顆粒之沖擊速度(m/s)；Y 為軟岩彈性模數(Pa)；kv為軟岩強度參數，該參數值介於 1012~1013， 需進行檢定；
σ
T為軟岩張力強度(Pa)；Ls為泥砂顆粒躍動長度(m)。 式(3-28)中軟岩彈性模數 Y 與軟岩張力強度σ
T 之值需由材料試 驗獲得。而泥砂顆粒之沖擊速度ws與泥砂顆粒躍動長度 Ls，Sklar and Dietrich (2004)藉由理論與試驗分析之結果，提出了 ws與 Ls之計算方 式，如下所示： 5 . 0 2 * 18 . 0 * 1 1 8 . 0 f c s w u u w (3-29) 5 . 0 2 * 88 . 0 1 1 8 f c s s w u D L (3-30) 式中，u*為剪力速度(m/s)， ⁄ ；τ 為作用於河床之剪應力(N/m2)；τ
c為泥砂顆粒臨界剪應力(N/m2)，wf為顆粒沉降速度(m/s)；Ds為泥 砂粒徑(m)。 式(3-28)中之 1 q /q 項反應軟岩表面之覆蓋效應，qs、qt分別 為泥砂非平衡狀態與平衡狀態下之輸砂量，前者為輸砂率(sedimenttransport rate)，後者為輸砂能力(sediment transport capacity)。若軟岩 表面覆蓋泥砂顆粒，則由顆粒跳動所引起之軟岩磨蝕量將減少；若輸 砂量等於輸砂能力，磨蝕作用將停止。 1 / 項反應懸浮效應， 當水流強度增大時，屬懸浮狀態之泥砂顆粒比例增高，此時泥砂顆粒 於軟岩表面跳動行為減少，磨蝕量則下降。同時，泥砂磨蝕率亦與 ⁄ 1 項成反比，反應其物理特性。
3-4 軟岩侵蝕模組數值方法
結合軟岩侵蝕模組之顯式有限解析法於模擬運算前，需先設定具 軟岩侵蝕特性之斷面，大安溪河段現場案例之設定區域以現地觀測及 航拍圖為主要設定依據。其後須給定軟岩侵蝕模組參數之計算方式， 大安溪河段現場案例以現地試驗實測各軟岩參數計算之沖蝕指數 Kh 為參數設定依據，該參數於模式中亦可給定其相關計算參數，由模式 於運算過程中求取各時刻之代表數值。至於軟岩參數如軟岩彈性模數 Y、軟岩張力強度σ
T，則以該區域之觀測資料為設定依據。 模式運算時，水理模組先進行運算，求取河道中各斷面之流量、 水位等水力特性參數，提供軟岩侵蝕模組計算所需之水理條件數值。 由於顯式有限解析法動床模式採非耦合運算，水理模組運算與動床模 組運算間隔不同，當時間增加至啟動動床模組運算間隔時，動床模組蝕模組計算所需之動床條件數值。軟岩侵蝕模組運算間隔與動床模式 相同，動床模組運算終了後，軟岩侵蝕模組開始運算，利用水理模組 與動床模組提供之數值為依據，計算該時刻之軟岩侵蝕率與侵蝕量， 與動床模組結合反應沖淤量於河道演變過程中，其後反覆運算至最終 模擬時刻。結束後動床模組與軟岩侵蝕模組之沖淤量即為顯式有限解 析法動床模式河道變遷最終結果。 模式於計算底床變化方面，在具備軟岩侵蝕特性斷面上，軟岩侵 蝕模組運算之侵蝕深度與動床模組運算之沖淤深度，因目前並無相關 文獻探討其相對關係，因此模式於同一時間間距下同時考量二模組之 沖淤量。側向斷面變化上，模式具備分離具軟岩特性區域與一般動床 特性區域之功能，可藉由參數設定分離深槽與灘地。但應用於大安溪 現場案例時，經由現地觀察，軟岩特性範圍遍佈整個斷面，屬一般沖 積特性之高灘地遠離河道且高程與河道高程差距甚大，模式模擬之水 深均無法溢淹至高灘地，因此於現場案例模擬時，軟岩特性斷面設定 為整個斷面。
第四章 模式檢定與敏感度分析
本研究探討軟岩侵蝕模組與河道動床數模結合後，是否優於一般 河道動床模式對於具軟岩侵蝕特性之河川模擬結果。但於比較軟岩侵 蝕模組加入前後之差異前，須先進行模式檢定與驗證，方可對其差異 性進行比較，以下說明檢定驗證案例。4-1 實驗室案例之檢定驗證
EFA 模式於實驗室案例模擬方面，採 Suryanarayana (1969)之動 床水槽實驗案例。陳(2006)改良並驗證 EFA 模式模擬此動床水槽案例 有相當精確度，採均質粒徑淤積案例 Run21 及均質粒徑沖刷案例 Run22 作為模式檢定案例；且根據檢定案例之各項參數，加以模擬驗 證均質粒徑淤積案例Run25、均質粒徑沖刷案例 Run26 以及非均質粒 徑沖刷案例Run14。於此簡介陳(2006)改良 EFA 動床模式之模擬成果， 並比較模擬與實測值之差異性。模式各項參數設定分述如下： 1. 初始渠道幾何資料 模擬之渠道為一矩形水槽，長18.3 公尺、寬 0.6 公尺，模型配置 如圖4-1 所示。數值模擬區域乃採水面線高程表格記錄資料之起點（第 五點）及終點（第 55 點）以作為模擬區域之上下游邊界。數值模擬以實驗起始時間所量測之底床高程作為初始底床高程。 2. 初始底床資料 初始底床質採Suryanarayana (1969)實驗記錄值，分布如圖 4-2 所 示，其中Sand 2 可視為均質粒徑，D=0.45mm；而 Sand 3 則為非均質 粒徑，採用三種代表粒徑D1、D2、D3 分別為 0.4、0.9、及 1.6mm， 來進行沖淤變化之模擬。此三種粒徑於各斷面之初始組成百分比，則 分別設定為0.34、0.33、0.33。 3. 粗糙係數 粗糙係數曼寧n 值以數值試驗檢定之，均質案例之曼寧 n 值研採 0.013，非均質案例之曼寧 n 值研採 0.016。 4. 上游入砂濃度 淤積案例中之上游入砂濃度，將記錄之輸砂量與流量換算可得 Run21 為 409 PPM，Run25 為 740 PPM。沖刷案例的部分則皆為清水 沖刷，上游入砂濃度為0 PPM。 5. 參考高度 Van Rijn (1984b)指出此參數可以是底床砂丘高之一半，亦可利用 糙度高度(roughness height)給定，其最小值為水深之百分之一。檢定
案例水深約為5-10cm 左右，檢定結果顯示，當此值採用 1mm 時，模 擬結果較符合實驗值，此值約為實驗水深之1/50~1/100，合於原物理 模型之假設。 6. 作用層厚度 根據數值經驗，設定此值介於0.05~0.1 倍水深間，以避免作用層 厚度過大，而造成粒徑變化無法反應真實之改變量。 7. 各案例上游入流量、下游水位資料 整理如表 4-1 所示。 8. 模擬結果與分析 8-1 檢定案例 (Run21、Run22) 圖4-3~圖 4-7 與圖 4-8~圖 4-10 分別為 Run21 與 Run22 案例模擬 之底床變化示意圖，圖中右端為上游端。觀察 Run21 結果圖可知， 底床變化與淤積波傳遞之過程相關，淤積速率與底床坡度相關。觀察 Run22 結果圖可知，模擬初期階段沖刷發生於上游段，隨時間增加往 下游段增加，使整體渠道高度逐漸下降，當平均坡度越趨平緩時沖刷 速度亦減緩。就上述模擬成果而言，EFA 模式預測底床淤積、沖刷結 果與實驗量測值相當接近。
8-2 驗證案例 (Run25、Run26、Run14)
圖4-11~圖 4-15、圖 4-16~圖 4-19 與圖 4-20~圖 4-22 分別為 Run25、 Run26 與 Run14 模擬各時程底床變化示意圖。觀察 Run25 結果圖可 發現，在模擬時間結束時，渠道上仍具有淤積之鋒面，此現象是因為 本案例之上游入砂濃度為740 PPM，較前述檢定案例之 409 PPM 大， 且上游入流量又較小，流速相對較小之故。觀察 Run26 結果圖可發 現，整體沖刷程度不如前述檢定案例，模擬出其甚至有淤積現象。此 二案例藉由與檢定案例不同入流量與下游邊界探討不同水流狀況之 下，底床變化現象，對於不同入砂濃度、水流條件、水位高程等變化， EFA 模式具相當不錯之模擬結果。非均質沖刷案例 Run14 除模擬初 始階段推估值較低之外，隨時間增加底床變化趨勢與實測值相當一致。 而非均勻案例之沖刷現象不只根據局部水理條件產生變化，還需考量 護甲或水力篩選作用等效應，透過此案例模擬結果，顯示 EFA 模式 於非均勻泥砂模擬過程，仍具有相當精確度。 綜合上述分析可知，利用本實驗渠道之實驗資料對 EFA 模式模 擬結果作定量之比較，可發現無論對於一般因水理流況改變、非均勻 河床質分佈或入砂條件改變等不同情境所引起之輸砂行為來說，EFA 模式皆能合理的模擬，詳細理論及模擬結果可參考陳(2006)之碩士論 文。
