一維顯式有限解析法應用於石門水庫水力排砂之研究

國 立 交 通 大 學

土木工程研究所

碩士論文

一維顯式有限解析法應用於石門水庫水

力排砂之研究

Study on Shihmen Reservoir Sediment Flushing

Using 1D Explicit Finite Analytic Method

研 究 生:康新詠

指導教授:葉克家 博士

中華民國九十九年九月

I 一維顯式有限解析法應用於石門水庫水力排砂之研究 學生：康新詠 指導教授：葉克家 國立交通大學土木工程研究所

摘要

台灣的水庫淤積嚴重，故水庫於颱洪期間進行水力排砂以維持庫 容乃當今水資源運用重要的課題。在各項水力排砂策略中，又以異重 流排砂策略最適用於石門水庫。本論文延續前人研究，將已植入凝聚 性沉滓傳輸機制的一維顯式有限解析法數值模式(許，2008)應用於石 門水庫進行異重流排砂之研究，以符合石門水庫凝聚性沉滓的特性。 根據前人研究，低高程排水口的洩水能力越大，則水庫進行異重流排 砂的效率越好。模式中應用 van Rijn (1984)所提之方法決定壩址處水 深方向的泥砂濃度，配合水庫水位流量率定曲線求得不同設施放流量 作積分，可算得其排砂量；將各設施之排砂量與上游入砂量相除可得 各設施之排砂效率。 本研究模擬 2008 年辛樂克颱風與 2009 年莫拉克颱風，使用數值 模式計算石門水庫異重流排砂策略的排砂效率，並與全洪程觀測資料 作比較。模式並以 2004 年艾利颱風為案例，評估石門水庫既有設施 防淤功能改善工程前後，不同的水庫低高程排水口洩水能力下，異重 流排砂的效率。 關鍵詞：顯式有限解析法、石門水庫、凝聚性沉滓、水力排砂

II

Study on Shihmen Reservoir Sediment

Flushing Using 1D Explicit Finite Analytic

Method

Student: Sheen-Yeong Kang Advisor: Keh-Chia Yeh

Abstract

The reservoirs in Taiwan have many problems related to sedimentation, so the sediment flushing during flood season to maintain the reservoir capacity is an important topic of water resources. Among the flushing strategies, venting density current is the most suitable one for Shihmen Reservoir. This study extends the previous research of Hsu (2008), which considered one-dimensional explicit finite analytic numerical model with cohesive sediment in reservoir, and applied the numerical model to the Shihmen Reservoir. According to previous studies, the larger capacity of the low-level outlet is, the better reservoir sediment flushing efficiency will be. The mass of flushing sediment of each outlet could be obtained by integrating the sediment concentration along water depth at the dam determined by the method of van Rijn (1984) and the outflow of each outlet determined by water stage-outflow rating curve. Therefore, the sediment flushing efficiency could be calculated by dividing the mass of flushing sediment to the mass of the upstream inflow sediment.

III

typhoon in 2008 and Morakot typhoon in 2009, the flushing efficiencies of simulation results were compared by the observation data. And by simulating Aere typhoon in 2004, the flushing efficiencies of venting density current with different low-level outlet outflow capacities before and after the modification of the flushing facilities were also evaluated.

Keyword: explicit finite analytic numerical model, Shihmen Reservoir, cohesive sediment, sediment flushing

IV

謝誌

研究所期間承蒙恩師 葉教授克家學業上的諄諄教誨以及生活上 的關心，使得兩年的研究生涯過得充實且有意義。感謝成功大學的蔡 教授、逢甲大學許教授和水規所的陳課長，願意將各位寶貴的時間及 心力用在我的論文審查上，使得拙作更臻完善。 感謝學長姊平常不吝於回答我的問題，以及生活上的照顧；感謝 同學們，明儒、浚騰、勁頤，我永遠不會忘記一起在風雨中拉鉛魚的 事蹟，都是共患難的兄弟啊!大頭、小邱、琦雯、群玲，這兩年來若 少了各位將失色不少。還有學弟妹們都主動詢問有沒有要幫忙的地方， 實在令人非常感動。 謝謝瑋琳的陪伴，也謝謝妳給我的鼓勵。 感謝大哥三不五時的加油打氣，還有你的經驗分享讓我更有勇氣 與信心面對挑戰。最後要感謝我最愛的父母，給我溫暖與支持。

V

目錄

摘要 ... I Abstract ... II 謝誌 ... IV 目錄 ... V 表目錄 ... VIII 圖目錄 ... IX 符號說明 ...XII 第一章 緒論 ... 1 1.1 前言... 1 1.2 文獻回顧 ... 4 1.2.1 水庫水力排砂相關研究 ... 4 1.2.2 石門水庫防淤相關研究報告 ... 7 1.3 研究內容 ... 16 第二章 數模理論基礎 ... 17 2.1 顯式有限解析法數值模式概述 ... 17

VI 2.2 水理控制方程式 ... 18 2.3 輸砂控制方程式 ... 19 2.4 輸砂輔助方程式 ... 20 2.4.1 河床載通量 ... 20 2.4.2 作用層厚度 ... 21 2.4.3 非凝聚性懸浮載源 ... 22 2.4.4 作用層源 ... 23 2.4.5 凝聚性沉滓輸砂輔助關係式 ... 23 2.4.6 凝聚性沉滓起動條件 ... 23 2.4.7 凝聚性沉滓沖刷公式 ... 25 2.4.8 凝聚性沉滓沉降速度 ... 25 2.4.9 凝聚性沉滓沉淤公式 ... 26 2.5 水體承載濃度剖面 ... 26 第三章 數值方法 ... 29 3.1 水理方程式 ... 29 3.2 輸砂方程式 ... 30

VII 3.3 一維顯式有限解析法應用於水庫模擬之假設 ... 34 第四章 模式檢定驗證 ... 37 4.1 石門水庫概述 ... 37 4.2 模式檢定 ... 38 4.2.1 民國 97 年檢定案例模擬條件 ... 38 4.2.2 民國 97 年檢定模擬結果分析 ... 42 4.3 模式驗證 ... 50 4.3.1 民國 98 年莫拉克颱風驗證案例模擬條件 ... 50 4.3.2 民國 98 年莫拉克颱風驗證案例模擬結果分析 ... 52 第五章 水力排砂效能評估 ... 59 5.1 石門水庫既有設施防淤功能改善工程規劃簡介 ... 59 5.2 水力排砂策略模擬 ... 62 第六章 結論與建議 ... 73 6.1 結論... 73 6.2 建議... 74 參考文獻 ... 76

VIII

表目錄

表 1-1 石門水庫淤積一覽表 ... 2 表 4-1 石門水庫水庫底質粒徑分析表 ... 41 表 4-2 辛樂克颱風期間排砂效率_模擬結果 ... 48 表 4-3 辛樂克颱風期間排砂效率_實測資料 ... 48 表 4-4 石門水庫淤積物清除數量統計表 ... 50 表 4-5 莫拉克颱風期間排砂效率_模擬結果 ... 57 表 4-6 莫拉克颱風期間排砂效率_實測資料 ... 57 表 5-1 模擬艾利颱風入庫各設施排砂量百分比_改善前 ... 60 表 5-2 模擬艾利颱風入庫各設施排砂量百分比_方案一 ... 61 表 5-3 模擬艾利颱風入庫各設施排砂量百分比_方案二 ... 62 表 5-4 艾利颱風期間排砂效率_改善前 ... 67 表 5-5 艾利颱風期間排砂效率_方案一 ... 69 表 5-6 艾利颱風期間排砂效率_方案二 ... 71

IX

圖目錄

圖 1-1 石門水庫累計淤積_民國 52 年至 97 年 ... 3 圖 1-2 Parker et al. 水工詴驗模型水槽 ... 7 圖 1-3 石門水庫底孔洩水能力與排砂比關係 ... 8 圖 1-4 石門水庫上游規劃排砂隧道工程佈置 ... 12 圖 1-5 分洪防淤隧道規模與排砂效率評估 ... 13 圖 1-6 分洪防淤隧道啟動水量與減供水量之關係 ... 13 圖 1-7 分洪防淤隧道啟動流量與年排砂量評估 ... 14 圖 1-8 石門水庫模擬結束底床變化_水位圖(艾利颱風) ... 15 圖 1-9 石門水庫模擬結束底床變化量圖(艾利颱風) ... 16 圖 2-1 凝聚性沉滓力平衡示意圖 ... 24 圖 2-2 不同粒徑及濕密度下之啟動剪應力 ... 25 圖 2-3 懸浮載泥砂平衡濃度剖面 ... 28 圖 3- 1 懸浮沈滓之移流特性軌跡 ... 31 圖 3-2 石門水庫下游邊界幾何斷面 ... 35 圖 4-1 石門水庫既有洩水結構平面佈置圖 ... 37 圖 4-2 石門水庫各進水口高程示意圖 ... 38 圖 4-3 石門水庫庫區斷面範圍 ... 40 圖 4-4 檢定案例上游邊界條件_民國 97 年颱洪 ... 40 圖 4-5 檢定案例下游邊界條件__民國 97 年颱洪 ... 41 圖 4-6 鳳凰颱風羅浮、電廠與石門大圳懸浮載粒徑分布 ... 42 圖 4-7 石門水庫下游邊界幾何斷面 ... 42 圖 4-8 辛樂克颱風石門水庫平均排砂濃度檢定結果 ... 43

X 圖 4-9 平均排砂濃度模擬值與實測值之比較 ... 43 圖 4-10 泥砂濃度剖面_辛樂克颱風 16hr ... 44 圖 4-11 泥砂濃度剖面_辛樂克颱風 24hr ... 45 圖 4-12 泥砂濃度剖面_辛樂克颱風 48hr ... 45 圖 4-13 各設施實測排砂濃度_辛樂克颱風 ... 46 圖 4-14 各設施模擬排砂濃度_辛樂克颱風 ... 46 圖 4-15 各設施出流量_辛樂克颱風 ... 47 圖 4-16 石門水庫模擬結束底床變化暨水位 ... 48 圖 4-17 石門水庫模擬結束底床變化量 ... 49 圖 4-18 驗證案例上游邊界條件_莫拉克颱風 ... 51 圖 4-19 驗證案例下游邊界條件_莫拉克颱風 ... 52 圖 4-20 莫拉克颱風石門水庫平均排砂濃度驗證結果 ... 53 圖 4-21 泥砂濃度剖面_莫拉克颱風 16hr ... 53 圖 4-22 泥砂濃度剖面_莫拉克颱風 24hr ... 54 圖 4-23 泥砂濃度剖面_莫拉克颱風 48hr ... 54 圖 4-24 各設施實測排砂濃度_莫拉克颱風 ... 55 圖 4-25 各設施模擬排砂濃度_莫拉克颱風 ... 55 圖 4-26 各設施出流量_莫拉克颱風 ... 56 圖 4-27 石門水庫模擬結束底床變化暨水位 ... 58 圖 4-28 石門水庫模擬結束底床變化量 ... 58 圖 5-1 石門水庫設施改善前示意圖 ... 59 圖 5-2 石門水庫設施改善方案一示意圖 ... 60 圖 5-3 石門水庫設施改善方案二示意圖 ... 61 圖 5-4 上游邊界條件_民國 93 年艾利颱洪 ... 62 圖 5-5 下游邊界條件_民國 93 年艾利颱洪 ... 63

XI 圖 5-6 水庫下游邊界幾何斷面_改善前 ... 64 圖 5-7 水庫下游邊界幾何斷面_方案一 ... 65 圖 5-8 水庫下游邊界幾何斷面_方案二 ... 65 圖 5-9 各設施排放泥砂濃度_改善前 ... 66 圖 5-10 各設施出流量_改善前 ... 67 圖 5-11 各設施排放泥砂濃度_方案一 ... 68 圖 5-12 各設施出流量_方案一 ... 68 圖 5-13 各設施排放泥砂濃度_方案二 ... 70 圖 5-14 各設施出流量_方案二 ... 70 圖 5-15 數模與物模排砂效率之比較 ... 71 圖 5-16 不同方案水庫底床變化量_艾利颱風 ... 72

XII

符號說明

一般符號 A :通水斷面積 a :砂丘高度之一半 B :渠道寬 k C _{:顆粒 k 之平均濃度} c C _{:代表凝聚性沉滓之懸浮載濃度}

c

:顆粒蔡司係數 w c :凝聚性沉滓含水率 c(z) :水深方向水體承載濃度剖面 :斷面平均泥砂濃度 D _{:凝聚性沉滓淤積造成之懸浮量} Dk :不同沉滓代表粒徑 Dm :不動之最小顆粒粒徑 d :顆粒粒徑 E _{:凝聚性沉滓沖刷造成之懸浮通量} m E _{:作用層厚度} g :重力加速度 h :水深 H1 :清水面以下水深 H2 :渾水水深 k :凝聚性沉滓強度常數 K1 :係數

XIII K2 :係數 K3 :係數 K4 :係數 w m :凝聚性沉滓之濕比重 d m :凝聚性沉滓之乾比重 n :曼寧值 nce(z) :無因次平衡濃度剖面 P _{:凝聚性沉滓沉淤機率} :孔隙率 Q _:流量 Q1 :一場洪水的平均入流量 Q2 :一場洪水的平均出流量 s Q _{:懸浮載之擴散通量(flux)} b k Q _{:粒徑 k 之河床載通量} l q :單位渠長之支流側流量 R :水力半徑 s :砂比重 c S :凝聚性沉滓之懸浮載源 Sf :摩擦坡降 k S _{:粒徑 k 之懸浮載資源項} a k S _{:粒徑 k 於作用層底部資源項} s  _{:疊代修正向量} T : 入庫洪水不排砂歷時 k T :輸送參數 TK _{:非均勻沉滓之代表粒徑數} p

XIV TC _{:凝聚性沉滓之代表粒徑數} t :時間 t1 : 入庫洪水持續時間 t  _{:時間間距} ' t  _{:模擬輸砂時間間距} U :水體流速 l u :單位渠長之支流側流速 * u :有效河床剪力速度 *c u _{:臨界剪應力} V :渾水庫容 W :Stoke 定理之沉降速度 c W :凝聚性沉滓沉降速度 lk w _{:顆粒向上躍起之速度} k w _{:非凝聚性沉滓 k 之沉降速度} x :沿渠道中心線之距離 x  _{:兩斷面間間距} Z :水位 b Z _{:底床高程} z: :水深方向座標 希臘符號  _{:沖刷權重因子}  _{:底床粒徑百分比} k  _{:作用層內粒徑 k 之百分組成} c  _{:作用層內凝聚性沉滓之百分組成} w  _{:動量校正係數}

XV γ :懸浮載與河床質載之比值 δa :參考高程 ζ :隱藏因子  :水體密度 b  _{:凝聚性土壤濕密度} :非凝聚性沉滓密度  _{:運動滯度}  : von Karman 係數 λH :H1/H2 λQ :Q1/Q2 λt :(T+V/3600Q2)/t1 b



:底床剪應力 τc :凝聚性沉滓起動剪應力 cd  :凝聚性沉滓沉降剪應力  :水黏滯係數 ηs :排砂比  _{:運動黏滯係數}  _{:絕對黏滯度} c  _{:連續方程式之時間加權因子} m  _{:動量方程式之時間加權因子} s 

1

第一章 緒論

1.1 前言 經濟部水利署於民國 99 年 2 月 4 日新聞稿中指出，台灣年平均 用水量為 180 億立方公尺，至 110 年將成長至 200 億立方公尺；雖然 一年降雨量為 2,510 公厘，是世界平均值的 2.6 倍，但可用水量僅為 世界平均值的 1/6，使得台灣成為名列第 19 的缺水國家。因台灣降雨 不均，使得水庫在水資源的調度上扮演了重要的角色。 石門水庫位於淡水河最大支流大漢溪上，水庫每日平均調蓄供水 量約 80 萬立方公尺公共用水，主要供應台北縣板新地區、桃園縣及 新竹縣湖口鄉之公共用水。每年用水量總計 9 億 5,000 萬立方公尺， 以民國 97 年 12 月有效容量 2 億 969 萬立方公尺推算，石門水庫每年 要滿庫運轉 4~5 次方可達到供水需求。 石門水庫於民國 52 年 5 月開始蓄水，同年 9 月即遇上葛樂禮颱 風；民國 53 年 3 月首次作淤沙測量時，庫內淤沙量已經達到 1,947 萬立方公尺，約為原估計年淤沙量 80 萬立方公尺之 24 倍。民國 85 年賀伯颱洪過後，庫區新增淤積量 867 萬立方公尺，水庫上游攔 砂壩大致淤滿，入庫泥砂量逐漸增大，尤其民國 93 年艾利颱洪過後， 水庫淤積量增加 2,788 萬立方公尺，民國 96 年韋帕颱洪期間巴陵壩 損毀，庫區淤積 962 萬立方公尺，累計至民國 97 年底水庫總淤積量 9,428 萬立方公尺。歷年石門水庫淤積一覽如表 1-1 所示，根據表 1-1 繪圖如圖 1-1 所示，可清楚看出極端颱洪事件(如 53 年葛樂禮與 93 年艾利颱風)是造成石門水庫淤積的主要來源。

2

表 1-1 石門水庫淤積一覽表

(資料來源：北水局，2009) 期別 起迄年月 間隔 _(年) 呆容量 10 3_M3 _{有效容量 10}3_M3 _{合 計 10}3_M3 _{年平均淤積 10}3_M3 剩餘 淤積 剩餘 淤積 剩餘 淤積 累計淤積 每期 累計 1 52.05-53.03 0.8 43,390 13,850 246,260 5,620 289,650 19,470 19,470 2 53.03-54.04 1.1 43,220 170 245,660 600 288,880 770 20,240 700 10,653 3 54.04-55.05 1.1 42,080 1,140 244,810 850 286,890 1,990 22,230 1,809 7,410 4 55.05-56.05 1.0 40,450 1,630 244,250 560 284,700 2,190 24,420 2,190 6,105 5 56.05-57.06 1.1 38,630 1,820 244,890 -640 283,520 1,180 25,600 1,073 5,020 6 57.06-58.05 0.9 37,840 790 244,350 540 282,190 1,330 26,930 1,478 4,488 7 58.05-59.06 1.1 34,680 3,160 242,480 1,870 277,160 5,030 31,960 4,573 4,501 8 59.06-60.06 1.0 34,610 70 241,140 1,340 275,750 1,410 33,370 1,410 4,120 9 60.06-61.12 1.5 30,000 4,610 240,520 620 270,520 5,230 38,600 3,487 4,021 10 61.12-63.08 1.7 29,767 233 240,518 2 270,285 235 38,835 138 3,437 11 63.08-64.11 1.2 28,877 890 241,108 -590 269,985 300 39,135 250 3,131 12 64.11-65.11 1.0 29,480 -603 238,477 2,631 267,957 2,028 41,163 2,028 3,049 13 65.11-66.11 1.0 28,440 1,040 238,730 -253 267,170 787 41,950 787 2,893 14 66.11-67.11 1.0 28,534 -94 238,130 600 266,664 506 42,456 506 2,739 15 67.11-68.11 1.0 28,446 88 238,034 96 266,480 184 42,640 184 2,584 16 68.11-69.11 1.0 27,296 1,150 238,357 -323 265,653 827 43,467 827 2,484 17 69.11-70.11 1.0 27,155 141 238,375 -18 265,530 123 43,590 123 2,356 18 70.11-71.11 1.0 26,199 956 238,439 -64 264,638 892 44,482 892 2,281 19 71.11-72.11 1.0 25,554 645 238,986 -547 264,540 98 44,580 98 2,175 20 72.11-73.11 1.0 25,515 39 238,956 30 264,471 69 44,649 69 2,077 21 73.11-74.11 1.0 24,320 1,195 236,456 2,500 260,776 3,695 48,344 3,695 2,149 22 74.11-75.11 1.0 24,615 -295 236,036 420 260,651 125 48,469 125 2,063 23 75.11-76.11 1.0 24,868 -253 235,560 476 260,428 223 48,692 223 1,987 24 76.11-77.11 1.0 25,171 -303 235,669 -109 260,840 -412 48,280 -412 1,893 25 77.11-78.11 1.0 25,330 -159 235,337 332 260,667 173 48,453 173 1,828 26 78.11-79.11 1.0 24,367 963 235,607 -270 259,974 693 49,146 693 1,787 27 79.11-80.11 1.0 24,786 -419 236,590 -983 261,376 -1,402 47,744 -1,402 1,675 28 80.11-81.11 1.0 24,708 78 236,670 -80 261,378 -2 47,742 -2 1,618 29 81.11-82.11 1.0 24,972 -264 237,387 -717 262,359 -981 46,761 -981 1,533 30 82.11-83.11 1.0 22,751 2,221 239,075 -1,688 261,826 533 47,294 533 1,501 31 83.11-84.11 1.0 22,328 423 239,599 -524 261,927 -101 47,193 -101 1,452 32 84.11-85.11 1.0 18,762 3,566 234,495 5,104 253,257 8,670 55,863 8,670 1,668 33 85.11-86.11 1.0 17,954 808 233,826 669 251,780 1,477 57,340 1,477 1,662 34 86.11-87.11 1.0 17,682 272 233,086 740 250,768 1,012 58,352 1,012 1,644 35 87.11-88.11 1.0 17,544 138 232,766 320 250,310 458 58,810 458 1,611 36 88.11-89.11 1.0 17,476 68 232,560 206 250,036 274 59,084 274 1,576 37 89.11-90.11 1.0 15,407 2,069 238,251 -5,691 253,658 -3,622 55,462 -3,622 1,441 38 90.11-92.04 1.5 15,651 -244 237,601 650 253,252 406 55,868 271 1,397 39 92.04-93.03 1.0 15,216 435 237,760 -159 252,976 276 56,144 276 1,369 40 93.03-93.12 0.8 6,216 9,000 218,876 18,884 225,092 27,884 84,028 34,855 2,010 41 93.12-94.12 1.0 4,925 1,291 217,835 1,041 222,760 2,332 86,360 2,332 2,018 42 94.12-95.12 1.0 5,943 -1,017 217,825 9 223,768 -1,008 85,352 -1,008 1,949 43 95.12-96.12 1.0 5,066 877 209,078 8,747 214,144 9,624 94,976 9,624 2,120 44 96.12-97.12 1.0 5,150 -84 209,692 -614 214,842 -698 94,278 -698 2,058 合計 45.8 52,090 42,188 94,278 容量損失 91.00% 16.75% 30.50%

3

圖 1-1 石門水庫累計淤積_民國 52 年至 97 年

0 10 20 30 40 50 60 70 80 90 100 52.05 -53.03 53.03 -54.04 54.04 -55.05 55.05 -56.05 56.05 -57.06 57.06 -58.05 58.05 -59.06 59.06 -60.06 60.06 -61.12 61.12 -63.08 63.08 -64.1 1 64.1 1 -65.1 1 65.1 1 -66.1 1 66.1 1 -67.1 1 67.1 1 -68.1 1 68.1 1 -69.1 1 69.1 1 -70.1 1 70.1 1 -71.1 1 71.1 1 -72.1 1 72.1 1 -73.1 1 73.1 1 -74.1 1 74.1 1 -75.1 1 75.1 1 -76.1 1 76.1 1 -77.1 1 77.1 1 -78.1 1 78.1 1 -79.1 1 79.1 1 -80.1 1 80.1 1 -81.1 1 81.1 1 -82.1 1 82.1 1 -83.1 1 83.1 1 -84.1 1 84.1 1 -85.1 1 85.1 1 -86.1 1 86.1 1 -87.1 1 87.1 1 -88.1 1 88.1 1 -89.1 1 89.1 1 -90.1 1 90.1 1 -92.04 92.04 -93.03 93.03 -93.12 93.12 -94.12 94.12 -95.12 95.12 -96.12 96.12 -97.12 累 計 淤 積量 (10 6m 3) 年份 石門水庫累計淤積圖 葛樂禮颱風淤積1,947萬立方公尺 賀伯颱風淤積867萬立方公尺 艾利颱風淤積2,788萬立方公尺 韋帕颱風造成巴陵壩潰壩淤積962萬立方公尺 水庫設計總容量309*106_m3 有效設計容量252*106_m3 呆容量57*106_m3 總容量損失30.5% 有效容量損失16.75% 呆容量損失91%

4 統計至民國 97 年 12 月，石門水庫總庫容已損失三成左右，呆庫 容更已淤積九成以上，意味著有效庫容將因入庫泥砂而損失得更快， 減少可供調度的水資源。因此減少石門水庫淤積、使水資源可充分利 用乃當務之急。 經濟部水利署水利規劃詴驗所民國 94 年至 96 年間觀測海棠、馬 莎、泰利、珊珊、聖帕、韋帕與柯羅莎等 7 場颱洪期間入庫泥砂運移 特性，發現通過電廠排放泥砂量隨著入庫尖峰流量與入庫水量增加而 增加。此 7 場颱洪入庫前水庫均處於高水位情況下(接近 240 公尺)， 入庫尖峰流量低於 1,600cms 時，水庫異重流濃度明顯沿程衰減，不 易運移至大壩排出；當入庫尖峰流量介於 1,600cms 至 2,800cms 之間， 水 庫 異 重流 在洪峰 後 可 運移 至大壩 排 放 ；當 入庫尖 峰 流 量超 過 2,800cms 時，水庫異重流在洪峰前運移至大壩排放。 據水利署民國 99 年 2 月統計，石門水庫入砂量平均每年達 353 萬立方公尺，人工清淤約 40 萬立方公尺，而由既有設施水力排砂量 則約為 73 萬立方公尺；且石門水庫入砂來源多為颱洪期間上游集水 區因暴雨而沖蝕或崩塌的土壤，隨著洪流帶入庫區，因此水庫防淤的 效率是颱洪期間進行水力排砂高於平時之人工清淤。 本論文承襲前人研究，以具有凝聚性沉滓傳輸機制的水理輸砂數 值模式，即一維顯式有限解析法(EFA)進行颱洪期間水庫水力排砂的 數值模擬，研究在不同的水砂條件下水力排砂的情況；並將模式應用 於石門水庫，評估不同水庫低高程洩水口放流能力對異重流水力排砂 效率之影響。 1.2 文獻回顧 1.2.1 水庫水力排砂相關研究

5 (1992)針對中國的三門峽水庫、官廳水庫、雙江水庫和劉家峽水庫庫 區內的泥砂沖刷及淤積情況，依據現地觀測資料進行定性及定量之探 討。泥砂在水庫中的行為主要可分為三項：(1)粗顆粒多集中沉積於 庫區上游，於水庫迴水區起點處形成三角洲。(2)隨水流進入庫區的 細顆粒沉滓在壩前沿程淤積。(3)洪水期間形成異重流，使清水在上 而渾水在下沿庫底往壩址移動。 Wang et al. (2009)對水庫泥砂所造成的問題進行討論，包括在迴 水區泥砂淤積造成底床抬升增高洪水位、泥砂進入水力發電設施造成 管路堵塞以及減少水庫蓄水量等；並且對於水庫泥砂的控制，提出了 以下策略： (1) 蓄清排渾：在非汛期，入庫泥砂濃度不高時水庫採高水位操 作，儲蓄清水；而在汛期間採低水位操作，降低水庫水位以 增加入庫水流挾砂能力，連水帶砂一併排至庫區下游。 (2) 異重流排砂：於颱洪期間，大量入庫水砂在水庫庫區中會形 成異重流現象，使得庫區中接近底床的水流挾帶了高濃度泥 砂，可利用水庫低高程洩水孔將運移至壩前的高含砂水流排 出，可有效降低水庫泥砂濃度以及減少泥砂淤積。 (3) 空庫排砂：將水庫水位降至底孔高程，使入庫水流可直接藉 底孔排出，除了將入庫泥砂隨水流排出，甚至可擾動水庫已 沉澱之淤泥，隨水流排出庫區，增加庫容。然而因為此種水 力排砂方式會於短時間內，將大量水庫沉泥帶到庫區下游， 使沉泥內富含之重金屬連砂帶水流入下游河道，為生態環境 帶來極大的負面影響。 (4) 人工清淤：使用機械挖泥、抽泥將淤積泥砂移出庫區，或擾 動庫底使淤泥可藉水流排出庫區。與前三項水力排砂策略比

6 較起來，人工清淤相對較無效率。 Khosronejad (2009)使用一維水理輸砂模式模擬伊朗 Sefid-Roud 水庫之水力排砂操作，並分別以 1980 年和 1981 年的現場觀測資料作 檢定驗證。確認此數值模式適用於 Sefid-Roud 水庫後，模擬颱洪事件 中不同的水庫底孔操作方式(水庫放流量)，提出最有效率的排砂操作 策略。 位於巴基斯坦的 Tarbela 水庫亦受到泥沙淤積的威脅，一旦泥沙 將壩址處的進水口堵塞，將使水電站及灌溉供水系統失去功能。Emma et al. (2000)應用”Hydro”套裝軟體，模擬 Tarbela 水庫於不同水砂條件 下，未來六十年的發電與灌溉效益；模擬結果顯示，設置潛沒堰保護 進水口避免泥砂淤積，以及低高程的排砂設施對水電系統供應無虞是 最有效的。 Sloff (1997)以 Tarbela 水庫及位於相同河道上的其它計畫興建水 庫為研究對象，同時考慮了泥砂模式與水資源調度模式，以達成放水 灌溉和水力發電最佳化之水庫操作策略。 Parker et al. (2007)進行水工模型詴驗(實驗水槽布置如圖 1-2)，模 擬水庫中的異重流過程，並以其觀測資料為基礎，發展了一套具異重 流機制的數值模式，用來推求水庫的囚砂率。

7

圖 1-2 Parker et al. 水工詴驗模型水槽

(資料來源：Parker et al. 2007) 賴(2004) 建立入庫泥砂觀測系統，實際調查水庫現場渾水異重 流垂直濃度分佈及評估設置減淤通道之可行性，並配合水理模型詴驗 與數值模式，探求適合於台灣水庫之減淤操作方式及減淤通道型式.。 1.2.2 石門水庫防淤相關研究報告 水利署為改善石門水庫淤積問題，主導多項研究計畫，以下摘要 相關計畫之研究成果。並且將許(2008)使用一維顯式有限解析法應用 於石門水庫之成果作一簡單說明。 (1) 石門水庫既有設施防淤功能改善工程計畫 (2006) 報告引用凌來水等 (1991)排砂比函數公式： 4 3 2 K t K Q K H 1 s



K















(1-1) 估算石門水庫的排砂比ηS。式中，ηS= 入砂量 排砂量 ，H= 2 1

H

H

，Q= 2 1

Q

Q

，

8 t= 2 1 3600 V T Q t  ，H1為清水面以下水深(公尺)，H2為渾水水深(公尺)， V 為渾水庫容(萬立方公尺)，t1為入庫洪水持續時間(hr)，T 為入庫洪 水不排砂歷時(hr)，Q1 為一場洪水的平均入流量(cms)，Q2 為一場洪 水的平均出流量(cms)；K1、K2、K3、K4為待迴歸參數；該報告引用 中國恒山、汾河水庫排砂資料驗算及對黑松林水庫排砂驗算結果，取 K1＝1.175，K2＝-1.026，K3＝-0.655，K4＝-0.164。 以 89 年碧利斯、90 年納莉、94 年馬莎、85 年賀伯及 93 年艾利 颱風分別代表 2 年、5 年、10 年、20 年及 100 年重現期洪峰流量， 推估不同底孔洩水能力之排砂比如圖 1-3。

圖 1-3 石門水庫底孔洩水能力與排砂比關係

(資料來源：水規所，2006) 藉由改建的底孔(永久河道放水口、電廠和排洪隧道)洩水量將來 砂排洩，降低水庫渾水的濁度及淤積效應。改善設施整體而言並不影

9 響石門水庫的洩洪能力，改建工程亦確保輸水路排除沉木的能力。永 久河道放水口、電廠設施及排洪隧道防淤功能改善後(底孔洩水能力 達 2,400 秒立方公尺)，預計排砂比將由原來之 16%增加至約 71%， 具有相當高之排砂效益。 (2) 石門水庫相關設施水工模型詴驗九十八年度期初報告 (2009) 水規所於民國 98 年以石門水庫為對象，進行水工模型詴驗。包 括電廠改建模型詴驗、增設取水工程詴驗、羅浮防砂設施詴驗、排洪 隧道改建詴驗以及水庫既有設施改善排砂減淤功能詴驗，分別簡介如 下： i. 電廠改建模型詴驗：因低高程洩水口排放能力不足，異重流來不 及排放容易造成渾水累積，水庫下層渾水置換水庫上層清水溢洪排 出。因此，石門水庫急需擴增下層通道之排放高濃度渾水能力，提 升水庫防淤功能，達到「蓄清排渾」之目的。石門水庫電廠更新改 善工程主要目標，即配合颱洪期間入庫渾水異重流抵達電廠進水口 ，擴大電廠排放流量由原本之 110cms 增為 462cms，以提升過庫泥 砂量。 ii. 增設取水工程詴驗：水庫異重流到達大壩附近來不及由下層排放 時，會累積造成下層取水之原水濁度過高，超過自來水公司處理能 力，導致無法供應民生用水與工業用水。因此，石門水庫急需增設 分層取水工工程，以利颱洪期間能分層取得濁度較低原水。本項增 設取水工工程包括上中下三層取水隧道，各層入口底部高程分別為 236 公尺、228 公尺及 220 公尺。 iii. 羅浮防砂設施詴驗：目前每年義興壩清淤量 25 萬立方公尺及羅 浮國軍清淤區 15 萬立方公尺，仍無法應付颱風豪雨沖蝕而下之土 砂量。因此，計畫於羅浮地區再增設攔砂設施，以有效攔阻入庫泥

10 砂及擴增羅浮國軍清淤量，以延長水庫使用壽命及維持水庫正常營 運功能。本項規劃設計目的為增加沉砂功能與擴增清淤量，經考量 施工、生態及經濟效益等，評估結果以「樁體丁壩」型式掛淤為較 佳方案。 iv. 排洪隧道改建詴驗：石門水庫排洪隧道改善工程初步規劃案，其 中主要目標係配合颱洪期間入庫渾水異重流抵達排洪隧道進水口 前，藉由入口高程 220 公尺向下延伸管線至標高 190 公尺，排放下 層異重流與渾水潭，擴增下層排放流量提升過庫泥砂量。 v. 水庫既有設施改善排砂減淤功能詴驗：評估現況既有設施改善功 能，包括發電進水口與永久河道放水口，以及排洪隧道延伸下層入 口之防淤功能，藉由現場觀測分析及水工模型詴驗驗證水庫異重流 防淤功能與相關設施水理功能，以及評估水庫更新改善工程對過庫 泥砂百分比之提升，提供相關研究、工程規劃設計與營運管理之參 考。詳細內容可參照 5.1 節。 (3) 石門水庫上游主河道分洪防淤工程初步規劃及水工模型詴驗研 究-初步規劃報告 (2010a) 計畫目標係藉水源量影響分析、分洪防淤功能探討、地質地形 地貌調查評估選線、工程方案研擬與配置、相關計畫之規劃設計案， 探討石門水庫增設防砂設施工程攔截粗顆粒泥砂，評估水庫上游分洪 防淤通道水理流況與分洪防淤功能，並結合水庫暨有設施更新改善工 程，提升過庫泥砂百分比，提供後續規劃設計參考。圖 1-4 為四個工 程規劃布置方案示意圖，水規所(2010a)就 A 案以物理模型詴驗進行 分洪防淤效能分析，其分析結果摘要如下。 圖 1-5 為模擬歷年來在石門庫區具代表性的颱洪事件，在不同規 模隧道的排砂效率。考慮到隧道規模成本與排砂效率，以分洪防淤隧

11 道規模設計為 1,600 秒立方公尺為最佳值。若是以民國 55~97 年共 35 場颱洪事件之平均值來看，則分洪防淤隧道在規模達 1,600 秒立方公 尺時，排砂效率為 61%。 若啟動時機為 300 秒立方公尺，則減供水量為每日九萬立方公尺 (圖 1-6)。圖 1-7 是設定不同隧道規模時，不同啟動流量可排出庫區的 泥砂量。在隧道規模為 1,600 秒立方公尺時，若啟動流量設定為 300 秒立方公尺，則年排砂量可達 129 萬噸。

12

圖 1-4 石門水庫上游規劃排砂隧道工程佈置

13

圖 1-5 分洪防淤隧道規模與排砂效率評估

(資料來源：水規所，2010a)

圖 1-6 分洪防淤隧道啟動水量與減供水量之關係

14

圖 1-7 分洪防淤隧道啟動流量與年排砂量評估

(資料來源：水規所，2010a) (4) 凝聚性沉滓傳輸機制之模擬應用研究 (2008) 利用一維顯式有限解析法數值模式模擬石門水庫庫區羅浮至壩 址處，分析泥砂在庫區中之運移、沖淤和分布規律，選擇民國 93 年 8 月之艾利颱洪事件，模擬庫區水體底床、濃度及輸砂傳輸的變化情 形。對於下游石門大壩，假設為虛擬結構物，即水流輸砂量均勻分佈 排出庫區。 圖 1-8 可見模擬範圍羅浮至壩址前整體河床及水位變化情形，圖 1-9 為底床變化量與實測結果之比較，與實測相比，整體趨勢也有相 當不錯的結果。而鑑於淤積面日漸升至發電進水口擋泥牆頂，石門水 庫管理局於民國 67 年起開始規劃壩前庫區之清淤計畫，因此民國 74 年起在下游段庫區以水力抽泥船清淤壩前 1.5 公里範圍之淤泥，每年 抽泥 300,000 至 450,000 m3_{，淤積面亦平均下降 7.2m，故為造成模擬} 結果在壩址前 1.5 Km 庫區差異之可能原因。在庫區下游部分，亦有

15 2 至 5 公尺的低估誤差結果，主要原因為斷面測量時間間距為一年左 右，而本案例模擬時間僅挑選民國 93 年造成石門水庫入流流況衝擊 最大的颱洪事件，相對於當年其他颱洪暴雨事件，假設對凝聚性沉滓 淤積沖刷影響力不大，以及凝聚性沉滓的沉淤速度極緩慢。經參考周 (2005)對於石門水庫中凝聚性沉滓沉淤詴驗歸納出石門水庫沉降速 度迴歸式，而推算每天石門水庫沉淤深度約數公分左右，因此當石門 水庫下游為靜水狀態時，年估沉淤深度也將達到三至四公尺，此為模 擬結果在石門下游沉淤低估的主因之一。對於水庫庫區長期下來皆處 於靜水狀態，而凝聚性沉滓在高濃度狀況下，長時間沉澱過程亦是造 成水庫淤積原因之一。

圖 1-8

石門水庫模擬結束底床變化_水位圖(艾利颱風) (資料來源：許，2008)

16

圖 1-9

石門水庫模擬結束底床變化量圖(艾利颱風) (資料來源：許，2008) 1.3 研究內容 由文獻回顧可知，蓄清排渾或空庫排砂等水力排砂方式，皆需要 降低水庫水位，影響水資源利用甚鉅；石門水庫扮演北台灣供水重要 角色，因此現行水庫排砂策略以異重流排砂為主，本研究探討之水力 排砂策略亦以異重流排砂為對象。 本研究延續已植入凝聚性沉滓傳輸機制之一維顯式有限解析法 動床數值模式，將其應用於石門水庫；藉由對水庫水理輸砂情況做各 種假設並加以簡化，可以一維水理輸砂數值模式模式推估水庫進行異 重流排砂之排砂比，可供水庫進行水力排砂規劃之參考。

17

第二章 數模理論基礎

2.1 顯式有限解析法數值模式概述

顯式有限解析法數值模式可分成水理及輸砂兩大部分。於水理計算方 面，在滿足 de Saint Venant 之基本假設下，應用顯式有限解析法（EFA）離 散動量方程式，並利用沿特性線積分的概念，求解動量方程式，再配合適 當的斷面處理與差分式，求解水流連續方程式。顯式有限解析法具有推導 容易，且程式撰寫較隱式法簡易之優點，而此法在計算流力與水力計算領 域之應用上已證實有相當不錯的成果。 在沉滓運移計算方面，因天然河床由多種不同粒徑之沉滓所組成，且 由於水庫迴水區甚長，能流淤在壩址前之沉滓一般均為細顆粒者，且經由 水庫防淤操作排放至下游河道之沉滓多屬細小之黏土或粉土，其特性將與 其下游河床上之床質有所差異，因此所採用之模式頇具有模擬非均勻沉滓 之特點，以反映水庫內淤積及排淤與河道原有沉滓之不同特性。此外，為 考慮懸浮載與河床載不同之運移機制，故將兩者予以分開計算，並考慮沉 滓在渠道底床附近發生沉淤與再懸浮之情形，因此模式引入懸浮載與河床 載間之交換機制，藉以推估水體中各懸浮沉滓之濃度變化，以及河床上床 質粒徑之組成。本研究採用適於雙曲線型方程式之特性法求解懸浮載質量 守衡方程式，並與河床載質量守恆方程式及整體河床輸砂之質量守恆方程 式進行結合演算，利用 Newton-Raphson 疊代聯立求解。 台灣水庫集水區內土壤質地多屬凝聚性沉滓，在水體中之運移機制較 粗顆粒複雜。其主要原因在於凝聚性沉滓具有凝聚結合力性質，使得顆粒 間產生碰撞而結合，或因重力效應而沈降，或是受底床剪應力效應分解而 懸浮。在如此動態過程中，很難再用以往粗顆粒輸砂方程式來適切描述此 物理現象。凝聚性沉滓運移機制包括啟動條件、沖刷經驗式、沉降速度公

18

式及沉淤經驗式；根據許 (2008)之研究成果，最適用於水庫之凝聚性沉滓 運移機制分別為 Wilbert et al. (2004)啟動機制、Krone (1999)沖刷公式、You (2004)沉降公式以及 Krone (1962)沉淤公式。 2.2 水理控制方程式 水理演算部分引用自葉等 (1996)，係根據 de Saint Venant 所推導之一 維緩變非穩流控制方程式，其基本假設如下： (1) 流速均勻分佈：流速均勻分佈在通水面積上，即每一個通水斷面積 僅存在一個流速，此即一維水流。 (2) 靜水壓分布：假設渠道中水流之垂向流線曲率很小而且忽略其垂直 加速度，因此水深方向速度梯度為零，可忽略垂向加速度，則假設成立。 (3) 可以類似定量流中之阻力公式 (resistance laws)考慮變量流邊界摩 擦及紊流現象。 (4) 底床坡度甚小：當假設成立時，重力沿渠道所造成的分力將會很小， 甚至可忽略不計，亦即水深可以垂向水面以水位高程及渠底高程差表 示。 (5) 忽略柯氏力及風力的影響。 (6) 水體密度假設為均勻分佈。 對於不可壓縮水流之控制方程式，包括水流連續方程式與水流動量方 程式，為如下形式。 水流連續方程式：

19  0      l q x Q t A (2-1) 水流運動方程式： 0 2                  l l f w gAS qu x Z gA A Q x t Q _ (2-2) 式中：A =通水斷面積；Q =流量；t =時間；x =主流方向之距離；g = 重力加速度；Z =水位；w=動量校正係數；ql=單位渠長之支流側流量(q l 為正屬合流之處理，q l為負屬分流之處理)；ul=支流在主流方向的速度分量； 3 4 2 2 R A n Q Q Sf  =摩擦坡降，其中 R =水力半徑，n =曼寧 n 值。 2.3 輸砂控制方程式 沉滓傳輸的方式可依據運動機制不同而分為滾動、滑動、跳躍以及懸 浮等形態，通常可將滾動、滑動、跳躍這三種型態合稱為河床載(bed load)， 或稱為推移載，一般來說，此類傳輸型態的共通點在於沈滓顆粒明顯受到 重力作用的影響，並使其運移範圍僅限於河床上方有限的高度，即所謂參 考高度內。然而，當紊流作用情況達到沈滓懸浮的起動條件時，部份沈滓 顆粒即向上運移超過參考高度，形成懸浮載(suspended load)傳輸型態；同樣 的，懸浮沈滓在受重力的影響下亦會逐漸沈降至底床，進而轉變為河床載 傳輸型態。沈滓會因所受作用力的不同而改變其傳輸機制，因此可將此種 特性視為輸砂行為中重要的物理機制。本研究中使用之輸砂模組部分引用 自許 (2002)，而凝聚性沉滓傳輸機制(包括啟動、沉降速度、淤積和沖刷等 機制)則引用自許 (2008)。 將輸砂控制方程式中河床質載分離為非均勻之懸浮載與河床載兩部份， 同時求解某一粒徑 k 之懸浮載、作用層質量守恆，及整體河床質載之質量 守恆等控制方程式，分別表示如下：

20 ( k ) (k ) ( k) k C A C Q C A S t x x x  _ _   _     k=1,2,...,TK (2-3) (1 ) ( k m) bk 0 k ak BE Q p S S t x           k=1,2,...,TK (2-4) 1 ( ) (1 ) ( ) 0 TK b bk k k BZ Q p S t  x       



 (2-5) 上三式中：Ck=某一代表粒徑 k 之懸浮載濃度；p=孔隙率；βk=作用層 內粒徑 k 之百分組成；Em=作用層厚度；Qbk=粒徑 k 之河床載通量；Sk=粒徑 k 之懸浮載源項；Sak=粒徑 k 之作用層源項；Zb=底床高程；B=渠道寬；TK= 非均勻沉滓之代表粒徑數。 至於凝聚性沉滓部分，則因為河床載之影響遠小於懸浮載部分，故在 此可將河床載通量假設為零，並考慮凝聚性沉滓所造成之懸浮載，故凝聚 性沉滓輸砂質量守恆方程式可表示如下: ( ) ( ) ( ) c c c c C A C Q C A S t x x x            c=1, 2,……TC (2-6) ( ) (1 ) c m 0 c ak BE p S S t        c=1, 2,……TC (2-7) 1 ( ) (1 ) 0 TC b c k BZ p S t      



(2-8) 上三式中：Sc =凝聚性沉滓所造成的懸浮載源；TC=凝聚性沉滓之代表 粒徑數；Cc=代表凝聚性沉滓之懸浮載濃度；βc=作用層內凝聚性沉滓之百分 組成。 本研究為簡化計算，假設沉滓為單一粒徑，則僅存懸浮載與河床輸砂 質量守恆方程式各一條，至於作用層質量守恆方程式則成為β=1 之恆等 式。 2.4 輸砂輔助方程式 2.4.1 河床載通量

21

量。其後，Spasojevic & Holly (1990)針對每一粒徑( )，將 van Rijn 公式稍 加修正。假設河床載運移僅發生在作用層內，其中某粒徑之百分組成表示

為_k，且考慮較小粒徑在水體中會形成懸浮載，可引入懸浮載與河床質載

之比值



(van Rijn, 1984b)，對河床載通量作修正。此外，在一般非均勻之

河床質中，較細顆粒可能被隱藏在較粗顆粒之間，而不易被水流帶動，故 Karim et al. (1987)提出一經驗因子，稱之為隱藏因子(hiding factor，ζ)，對河

床載通量予以修正。綜合上述之影響因子可得粒徑為 之河床載通量如 下： 2.1 0.3 * (1 ) t( ) (1 ) (0.053) ( 1) k bk k k b k k j k k k T Q q D B s D D B D            (2-9) 上式中， 13 * 2 ( 1) [ ] k k s g D D    (無因次顆粒粒徑)； 2 2 * * 2 * ( ) ( ) c k k c k u u T u   (輸送參數)； * u g u c  (有效河床剪力速度)； 90 12 18log( ) 3 d c D  (顆粒蔡司係數)； s s    (砂之 比 重 ) ；  = 運 動 滯 度 ； u*c = 臨 界 剪 應 力 ； 0.85 50 ( k ) k D D   ； 上式中， w_k=粒徑k之沉降速度。 2.4.2 作用層厚度

沖刷現象發生時，根據 Bennet and Nordin (1977)之研究， E_m可以下式 表示： 1 1 ( ) n n n m b b E   C Z  Z (2-10) 式中， C為數值參數(模式設定為20)。當河床表面接近護甲條件時 (armored condition)，作用層厚度接近零，在這種情況下，可用 Borah et al.

D_k D_k                 4 . 0 0 10 4 . 0 ) ln( 325 . 0 25 . 0 10 1 * * * * k k k k w u w u w u w u 

22 (1982)所提出護甲層之厚度(armored-layer thickness)，予以修正： 1 1 1 ( ) 1 n n n m m b b k k k m D E C Z Z p         



(2-11) 式中，D_m為不動之最小顆粒粒徑。另外，作用層在淤積期間可定義為： 1 1 ( ) n n n n m m b b E   E  Z  Z _(2-12) 上式中之上標代表n及(n+1)時刻。 2.4.3 非凝聚性懸浮載源 懸浮載源係由懸浮質向下之通量與底床亂流剪力作用產生河床質向上 之通量交互作用之結果。使懸浮質下移到河床表面，主要是受到重力之影 響。懸浮質之沈降速度決定向下沈淤之通量，所以對某一粒徑 k 之懸浮質 而言，其向下之通量可表為： dk k dk Q Bw C (2-13) 式中， * [3.25 0.55ln( k )] dk k w C C u    (Lin,1984)；C_k=顆粒 k 之平均濃度；



=von Karman 係數。另一方面，床面沈滓成為懸浮質，主要受到底床之亂 流作用所造成。Bennet and Nordin (1977)認為對某一粒徑 k 而言，河床質向 上之通量可表為： Qek=BwlkCak (2-14) 式中， * 3.2 4.5 log lk w

u     1.2 (Hu & Hui, 1996)；

=3.1  1.2 ( ) b s gD       1.5 0.3 * 0.015 k k ak k D T C a D  (Van Rijn, 1984b)；

23 lk w =顆粒向上躍起之速度；



b=底床剪應力；a=砂丘高度之一半。 由(2-13)及(2-14)式知，對某一粒徑k之懸浮載源可表為： Sk B w C( lk ak w Ck dk) (2-15) 2.4.4 作用層源 作用層源係表示介於作用地層(active stratum)控制體積間河床粒徑改變 量，由於作用地層頂面之升降而產生，當其下降時， (1 ) [ ( )] ak sk b m S p BZ BE t        (2-16) 其中，_sk為作用地層內某一粒徑之百分組成。如作用地層之厚度增加， 即其頂面上升時，(2-16)式中之_sk則改為_k。 2.4.5 凝聚性沉滓輸砂輔助關係式 輔助機制乃屬懸浮載源部份，懸浮質為向下之通量與底床亂流剪力作 用產生河床質向上之通量交互作用之結果。運移機制包括啟動條件、沖刷 公式(向上懸浮通量)、沉降速度、淤積公式(向下沉淤通量)。目前探討凝聚 性沉滓運移過程中，不計算河床載部分而只針對凝聚性沉滓所造成之懸浮 載源作計算: 凝聚性沉滓懸浮載源S_c B(_cE D ) /100_{s c=1, 2,……TC (2-17)} 上式中， E _{:凝聚性沉滓沖刷造成之向上懸浮通量(} 2 / g cm s )； D _{:凝聚性沉滓淤積造成之向下懸浮通量(} 2 / g cm s )。 在模式中定義粒徑小於等於 0.01mm 屬於凝聚性沉滓，大於 0.01mm 屬 於非凝聚性沉滓。 2.4.6 凝聚性沉滓起動條件 對於凝聚性沉滓起動方面除了頇考慮重力、升力、阻力外，還必頇考 慮黏滯力，Wilbert et al. (2004)根據 Roberts et al. (1998)實驗數據，迴歸出一

24 套半理論半經驗的凝聚性沉滓啟動機制。凝聚性沉滓力矩平衡示意如圖 2-1。

圖 2-1 凝聚性沉滓力平衡示意圖

(資料來源：Wilbert et al. 2004) 利用力矩平衡F b_d (F_g F a_l) ，並考慮黏滯力F_c可以推求出起動剪應力τc: 3 4 2 3 ( g c) (1 )(0.414 10 ) c cn g F F _c d F c d        (2-18) 上式中，c₃   ( _s ) / 6g ；c₄ 4 ( ,5) 10 c b      ，而 1 1 ( , 5) b b c b a e     表示為粒 徑 5 m時根據實驗數據推導出不同濕密度下之迴歸式，其中a₁=7 8 2 10 N m/  ， 1 9.07 / b  L Kg，最後可以推算出如圖 2-2 對應在不同粒徑及濕密度下之啟動 剪應力。

25

圖 2-2

不同粒徑及濕密度下之啟動剪應力

(資料來源：Wilbert et al. 2004)

2.4.7 凝聚性沉滓沖刷公式

Krone (1999)採納 Roberts et al. (1998)實驗結果迴歸出經驗公式，認為沖 刷速率和土壤濕密度存在著線性關係，在土壤濕密度為 1.77g/cm3 時，其土 壤沖刷率大小之分界，公式如下: 3 2 1.84 10 (1.80 b) E     ρ s<1.77g/cm3 (2-19) 4 2 3.65 10 (1.92 b) E     ρ s>1.77g/cm3 (2-20) 2.4.8 凝聚性沉滓沉降速度 沉降速度在凝聚性沉滓沉積佔了重要角色，顆粒在水中可能相互結合 成團絮而影響沉降速度， You (2004)認為凝聚性沉滓沉降非顆粒大小為主 要原因，其紊流、濃度、鹽度亦是重要因素，因此可推得一沉滓沉降速度 與濃度關係之公式：

26





0 0 ( , ) ( , ) ( ) h c c c t t t C z t C z t t dz W C t dt     





(2-21) 迴歸後其經驗公式 ' 2 exp( 6.1952 0.9779 0.108 ) c c c W t C C h  _{   } (2-22) 上式中， ' t  :模擬輸砂時間間距；h :水深。 2.4.9 凝聚性沉滓沉淤公式 當水流剪應力低於沉降剪應力時，其凝聚性顆粒將落淤於河床上；本 研究採用 Krone (1962)之凝聚性沉滓之沉淤公式: /10 c c D W  C P _(2-23) (1 ) cd P     (2-24) 式中 P 為沉淤機率因子，_cd為凝聚性沉滓沉降剪應力，在物理上為沉 滓對於底床的附著力，由於目前乃無學者對此機制有一完整定義，一般以 率定方式給之。Krone (1962)利用沉淤機率因子 P 判斷淤積與否，當 大於_cd 時，P 為零，即表示無沉淤產生。 2.5 水體承載濃度剖面 所謂水體承載濃度剖面，指的是水流中在參考高程 δa以上懸浮沉滓的 分布情形；在垂直二維或三維模式中，此濃度剖面的發展可透過懸浮滓在 水深方向的延散機制來模擬。然而在一維模式中不存在水深方向的維度， 無法模擬濃度剖面的發展過程，因此模式中簡單假設水體承載濃度剖面與 平衡濃度剖面具有相同的無因次濃度剖面(van Rijn, 1984b)，表示如下： c(z)= c〃nce(z) (2-25)

27 z 為水深方向座標，以底床為原點，c是斷面平均泥砂濃度，nce(z)為無 因次平衡濃度剖面；無因次平衡濃度剖面定義如下式： nce(z)= e e c (z) c (2-26) 上式中的 ce(z)為平衡濃度剖面(圖 2-3)，在平衡濃度的狀態下，沉滓向 上擴散與向下沉降的速率相等。ce(z)可以下式表示： ce(z)=ca[ ) ( _a a h    ] Z e-4Z(z/h-0.5) z/h≧0.5 (2-27a) ce(z)=ca[ ) ( ) ( a a h z z h     ]Z z/h<0.5 (2-27b) 式(2-27)中的 ca為參考平衡濃度，以下式表之： ca=0.015 0.3 * 5 . 1 50 D T D a  (2-28) 其中 D50=中值粒徑，h=水深， 參考高程 δa設定為 0.03 倍水深 (van Rijn, 1986) ， T= ₂ * 2 * 2 * ) ( ) ( ) ' ( c c u u u  = 傳 輸 參 數 ，u*' = ' c g u = 顆 粒 底 床 剪 力 速 度 ，c' =18log( 90 3 12 D R_b )=顆粒蔡司係數，Rb=底床水力半徑，D90=90%沉滓小於的粒徑， c u_* =Shields 臨界底床剪力速度， 3 1 2 50 * ] ) 1 ( [  g s D D   =無因次沉滓粒徑，s=  s ₌ 砂比重，ρs=沉滓密度，ρ=清水密度。 式 (2-27) 的 Z=ws/σκu* +ζ= 懸 浮 沉 滓 參 數 ， ws= 沉 滓 沉 降 速 度 ， σ= ) ( ) ( 流體擴散係數 沉滓擴散係數 f s   =比例係數，κ=von Karman 常數，u*為底床剪力速度， ζ=2.5( * u w_s )0.8( 0 c c_a )0.4懸浮沉滓參數修正值，c0=0.65 最高沉滓濃度體積比。應 用描述由參考高程至自由水面平衡濃度剖面之式(2-27)，可得平衡濃度剖面

28 之水深平均濃度c_e，表示如下式： e c = c dz h h z z e a



a    ₎  ( 1 (2-29)

圖 2-3 懸浮載泥砂平衡濃度剖面

(資料來源：van Rijn,1984b) 模式使用 Jacobian 係數矩陣與 Newton-Raphson 法，可疊代求解得斷面 平均泥砂濃度。將此平均濃度代入式(2-25)，可求得水深方向的濃度剖面， 亦即不同高程處的懸浮載泥砂濃度。假設運移至壩址前之泥砂可藉由排水 設施排出庫區，則壩址處不同高程的泥砂濃度即可近似表示為石門水庫不 同高程放水口所排放的泥砂濃度。

29

第三章 數值方法

河道水理計算之控制方程式為非線性聯立方程組，本數值模式沿用葉 等(1996)成果進行河道水理演算。水理模式採用顯式有限解析法(EFA)進行 水理控制方程式之離散化，此數值方法主要係用以求解雙曲線型偏微分方 程式，符合移流項之數學形式。輸砂模式方面，引用許(2002)之成果，採用 與水理分離演算(uncoupled)的計算方式，即在每一計算時段內先求解水理條 件，再以此推估輸砂量與底床沖淤量等，反之水理條件受到輸砂行為的影 響則在累進時間的過程中反應。 3.1 水理方程式 水理連續方程式以控制體積的觀念差分，用以求解通水斷面積變化量。 動量方程式則因其具有雙曲線型方程式之特性，故針對移流項部份採用顯 式有限解析法予以處理。 經離散後之水理連續方程式如下：





0 2 1 2 1 1 1 1 1 1 1 1                                     _{ }     x q Q Q x q Q Q t A A _in _in _lin c n li n i n i c n i n i (3-1) 式中， n1 i A 為未知數，上標為(n+1)者，係先給定 n 時刻之量測值，經反 覆疊代後，再將(n+1)時刻所計算之值帶入；



c 與以下的m分別為連續方

程式與動量方程式之時間加權因子(time weighting factor) ，其範圍在 0~1 之間。qli為支流流量，合流時qli為負，分流時qli為正。

30









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-2) 式中， n1 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）投射，此為 顯式有限解析法求解特色之一。 3.2 輸砂方程式 將懸浮載質量守恆方程式、作用層質量守恆方程式及整體河床質載守 恆方程式，利用 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-3)

31 整體河床質載守恆差分式：









      _       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 b k n i b k Q S S Q    _(3-4) 上式中，S=Sk+Sc。 distance, x time, t t n t  t n t( 1) t  2 t  1 t  LNS t  c D c l=1 l=2

cur r ent computed node

x



char acter istics tr ajector y l=LNS+1 l=LNS l=3

圖 3- 1

懸浮沈滓之移流特性軌跡 圖 3- 1 是以一維空間為例示意沈滓之移流軌跡。就_tn1 計算時刻之計算 點 A 而言，移流軌跡存在於_tn1 與_tn 時刻之間，定義_tn1 時刻之端點 A 為到達 點(arrival point)，_tn 時刻之端點 D 為離開點(departure point)。離開點 D 之懸 浮沈滓濃度可由初始條件求得，但由於模式採用固定格點，故離開點 D 並 不能保證剛好落在格點上，因此該點之濃度必頇藉由鄰近格點濃度以內插 的方式求得。 此外，在一般情況下，河床質移動量遠小於懸浮質移動量，即兩者在

32 時間尺度上相差甚大，在必頇聯立求解的前提下，懸浮載方程式頇使用較 大之可蘭數(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-5) 其中，下標l_{為節點編號，}l_{=LNS+1 代表到達點 A，}l_{=1 代表離開點 D。} 利用式(3-5)推求各節點位置(xl)，必頇要先知道各節點上的移流速度(ul)，但 移流速度又與節點位置有關，可利用疊代收斂的方式來推求一正確的移流 軌跡。 懸浮載質量守恆特性方程式： 1 1 1 ( ) ( ) [( ) ] 2 LNS k k k arr k dep k k c k t t C A C A S _ S S     



   1 1 1 1 1 1 1 2 2 (1 ) [( ) 2( ) ( ) ] [( ) 2( ) ( ) ] ( ) ( ) n n n n n n k i k i k i k i k i k i A A C C C t C C C t x x                     x u u A C i i i _  _  2 1 1 (3-6) 由上式可清楚瞭解到，在斷面一維模式中沈滓交換速率 S 所造成之影 響會直接反應在泥砂濃度的改變上，應被視為一源項，為影響沈滓濃度分 佈的重要物理量。 經由以上數值離散處理後，在非均勻沈滓共區分為 TK 個粒徑區間的情 況下(TK≧2)，任一計算點共可得(2TK+1)條代數關係式，包括 TK 條懸浮 載質量守恆離散式、TK 條作用層質量守恆離散式及 1 條整體河床輸砂質量 守恆離散式。

33 在計算點上之未知量可以向量形式表式如下： 1 1 1 1 ) , ,..., , ,..., , , (   _ n TK TK k k b n c c c z s    (3-7) 或更簡潔地寫成： 1 1 2 2 1 1 ) , , ( _   _ n k k n s s s s k=1,2,…,TK (3-8) 其中，k為粒徑區間之標號。則整體河床輸砂質量守恆離散式、作用層 質量守恆離散式與式(3-6)可分別寫成： 0 ) ( 1 1   n s F  (3-9) 0 ) ( 1 2   n k s F  k=1,2,…,TK (3-10) 0 ) ( 1 1 2    n k s F  k=1,2,…,TK (3-11) 以上三式為非線性代數式，可線性化後利用 Newton-Raphson 法疊代求 解： ) ( ] [ 1 1 1      _{m n} s F s s F    (3-12) ) ( ] [ 2   ₂ 1   _{m n} k k s F s s F    (3-13) ) ( ] [ 2 1 _{2 1} 1   _{ }   m n k k s F s s F    (3-14) 式中，F s為 Jacobian 係數矩陣中之列向量；msn1為前一次疊代所得 之向量，上標為疊代計數；s為疊代修正向量，可表為s(s1,s2k,s2k1)  。 解得修正向量s後，可得新的m1sn1向量： s s sn m n m1 1   1  (3-15) 當 s小於某一收歛容許值時，疊代方可結束，本模式設定收歛值為 0.01。

34 3.3 一維顯式有限解析法應用於水庫模擬之假設 水庫泥砂運移行為乃三維過程，本研究以一維的數值模式模擬其水砂 運移情況，作了以下假設： (1) 洩水設施條件 河道放水口之高程與電廠相近，且其設計流量只有 34 cms，故將河道 放水口併入電廠考慮。石門大圳設計流量 13 cms，相對流量甚小，排砂量 亦少，佔整體排砂比不到 1%，故忽略不計。電廠與河道放水口高程 173 公 尺，合計流量 110 cms；排洪隧道高程 220 公尺，設計流量 1600 cms。小於 1710 cms 的出流量由電廠及河道放水口排洪隧道洩出，而大於 1710 cms 的 出流量就經由高程 235 公尺的溢洪道洩出。 一維模式中，無法表現不同高程洩水口的出流量情況，故以石門水庫 各設施之設計流量及基本的水理運算來假設下游邊界幾何斷面，以近似真 實情況中，壩址處不同高程洩水口對應到不同的流量。以電廠與河道放水 口高程 173 公尺，合計流量 110 cms 為例，石門水庫設計滿水位為 245 公尺， 以 v=√2𝑔ℎ估算高程 173 公尺處之流速為 37.6 公尺/秒，再以 Q=Av 計算流 量通過之斷面積約為 3 平方公尺，則可以假設高程 173 公尺處之斷面寬度 約為 5 公分。以此法可推得排洪隧道處之斷面寬度約為 3 公尺，加上電廠 與河道放水口的 5 公分為 3.05 公尺。由於假設大於 1710 cms 之流量由高程 235 公尺之溢洪道排出，故溢洪道之寬度不再作計算。 依照以上假設條件，調整數值模型出流處之斷面幾何形狀如圖 3-2。

35

圖 3-2 石門水庫下游邊界幾何斷面

(2)水深方向濃度 因為水庫異重流現象，高含砂渾水會在庫區底層運移，使得靠近底床 之泥砂濃度較高。一維模式無法模擬三維的水庫泥砂分布情況，故假設水 深方向濃度可以 van Rijn (1984b)所提出的平衡泥砂濃度剖面推算；由模式 所算出的斷面一維平均泥砂濃度，推估水面下不同高程處的泥砂濃度(可參 考 2.5 節)。 (3)單一粒徑 應用 van Rijn (1984b)的平衡泥砂濃度剖面，則不同泥砂代表粒徑可得 不同的濃度剖面。為簡化排砂量計算，選取單一代表粒徑，使最下游斷面 處可求得單一之泥砂濃度剖面。而藉水庫既有洩水設施所排放的泥砂，亦 假設為單一粒徑。本研究探討目標為水庫水力排砂之整體效率，並不對所 排放泥砂的粒徑百分比作討論，故在模式中假設石門水庫之泥砂為單一粒 徑尚稱合理。 (4)排砂量計算與流量假設 將壩址不同高程處之泥砂濃度乘上對應該高程處的洩水設施流量，即可 170 180 190 200 210 220 230 240 250 0 50 100 150 200 250 300 350 高程 (m ) 累距(m)

水庫下游邊界幾何斷面

36 得該高程設施單位時間所排放的砂量。累加模擬時間內不同洩水設施每一 時刻之排砂量，即可得不同洩水設施之總排砂量，並分別將其除以上游入 砂量，則可得不同洩水設施之排砂效率。在辛樂克颱風、莫拉克颱風與艾 利颱風之洩水設施改善前的案例中，流量採用北區水資源局各洩水設施逐 時流量實測資料；而艾利颱風之洩水設施改善方案一與改善方案二的案例 中，則是根據實測之上游入流量，在壩址處以不同洩水設施放流能力作分 配，得假設之各洩水設施放流量。

37

第四章 模式檢定驗證

4.1 石門水庫概述 石門水庫位於桃園縣龍潭鄉之大漢溪主流上，集水面積 763.4 km2_。水 庫於民國 45 年 7 月開工，民國 52 年 5 月開始蓄水，至民國 53 年 6 月施工 完成。石門水庫為一多標的水庫，具有防洪、灌溉、給水及發電功能。石 門壩頂標高為 252.1 m，最大壩高 132.1 m，壩頂長 360 m，壩身體積 706 萬 m3_{。水庫滿水位標高為 245 m，呆水位標高 195 m。} 石門水庫主要設施包括大壩、溢洪道、排洪隧道、發電進水口、桃園 大圳進水口、石門大圳進水口、河道放水口、後池及後池堰等；相關設施 位置分布及大壩洩水設施高程位置分布如圖 4-1 和圖 4-2 所示。

圖 4-1 石門水庫既有洩水結構平面佈置圖

38

圖 4-2 石門水庫各進水口高程示意圖

民國 93 年 8 月艾利颱風襲台，因極端降雨造成石門水庫集水區大量的 土壤沖蝕與崩塌，隨著洪流進入庫區內，使得水庫大量淤積。民國 93 年 12 月北水局進行庫區斷面測量，與民國 92 年 12 月之紀錄比較，石門水庫總 庫容約減少了 2,788 萬 m3_{；統計至民國 97 年 12 月，總庫容約剩餘 2 億 m}3_， 較原始設計庫容 3 億立方公尺約淤積了 30%。 將模式應用於民國 97 年與 98 年颱洪事件，比較底床變化與排砂濃度 的實測值與模擬值，以驗證一維顯式有限解析法數值模式於石門水庫之適 用性。 4.2 模式檢定 4.2.1 民國 97 年檢定案例模擬條件 (1) 斷面資料：以北水局提供之民國 96 年 12 月石門水庫大斷面為初始 底床，範圍自上游斷面 32(羅浮)至下游斷面 3(石門壩址，斷面 3 為 最靠近壩址處，斷面 1、2 不在庫區主深槽上)。斷面範圍如圖 4-3 所示。