DOI: 10.29417/JCSWC.201912_50(4).0001
簡易一維土石流型堰塞湖潰決演算及運動波洪水傳遞模型：
以荖濃溪為例
陳慈愔[1] 洪啟耀[2]*
摘 要 2009 莫拉克風災造成大量崩塌土砂堆積於台灣山區支流集水區，形成後續土石流事件之料源。大 規模之土石流將擴大主支流匯流口處之沖積扇，並可能進一步淤塞主流形成堰塞湖，進而對下游村落產生威 脅。本研究以土砂連續方程式與水流連續方程式為基礎的潰壩模式模擬堰塞湖之形成與潰決過程，再以一維 非均勻河道運動波模型進行洪水演算，探討堰塞湖潰壩對下游保全對象之影響。本研究所提出之演算模型僅 需少量可由現地資料率定之參數即可進行高效率且可靠之模擬，可同時將主流上游流量與潰壩所產生的洪峰 流量一併納入計算，亦可在精簡的計算中模擬出河道寬度與坡度變化對洪水傳遞速率與洪峰削減之影響。本 研究以荖濃溪與其支流布唐布納斯溪交會口之堰塞湖做為案例，採用蒙地卡羅法配合現地歷史降雨資料和率 定完成之降雨逕流模式生成大量洪水事件作為堰塞湖潰決模式之上游邊界條件，再以系列模式模擬之結果進 行分析，得到不同情境下此區域兩百年重現期之堰塞湖潰壩洪水水位預測值。最後，根據分析結果，提出針 對台灣山區洪水事件之危害性評估建議。
關鍵詞：堰塞湖、蒙地卡羅法、潰壩模式、運動波輸送理論
Simple One-Dimensional Dam Breach and Flood Route Models for Tributary-Dammed Lake: A Case Study of the Laonong River in
Southern Taiwan
Tzu-Yin Chen[1] Chi-Yao Hung[2]*
ABSTRACT Due to numerous landslides caused by Typhoon Morakot in 2009, much loose sediment was deposited in the upper reach of the tributary of Laonong River. Loose material was transported into the main stream during the heavy rainfall of the typhoon and formed a tributary-dammed lake. If the dam were to fail, this would pose a significant threat to downstream villages. In this research, we present a diffusion equation–based one-dimensional semianalytical model of a tributary-dammed lake breach. To estimate the most dangerous scenario, the model incorporates the intense discharge caused by torrential precipitation during typhoons. To assess the security of downstream villages, the floods are routed by a one-dimensional kinematic wave solution for a nonuniform valley. The solution can simulate areas impacted by flooding with much reduced computational cost and topographical data requirements. Finally, the research applies combined methods to the tributary-dammed lake at the confluence of the Putunpunas tributary river with the Laonong River. The maximum flood elevation of the dam was estimated for breach floods under different scenarios.
Key Words: tributary-dammed lake, breaching model, kinematic wave model
一、前 言
台灣平均每年有三至四個颱風襲擊，強降雨所引發的土 砂災害事件往往造成人民財產大量的損失。而在莫拉克風災 時，台灣南部受到重創，以超過 3000 mm 的總雨量與高達 1500 mm/day 的降雨強度造成許多支流集水區發生大規模的 崩塌。其所產生的土砂材料大量的輸送進入主流，除了造成 主流河床大幅度的變動，也在主流上形成堰塞湖，而其潰決
的巨大洪峰也對下游村落造成極大的損失。然而，在風災過 後，由於支流上游仍有因崩塌所堆積的大量土砂材料，每當 颱風豪雨事件侵襲時，仍有相當高的機率會透過支流土石流 匯入主流形成堰塞湖，對下游村落、公路都產生相當大的威 脅。
針對堰塞湖過去的研究，多以崩塌所形成的堰塞湖為主 要討論的對象，相對於土石流型堰塞湖主要在壩體形狀與潰 決方式有所區別。土石流型堰塞湖之壩體由於受到主流河水
〔1〕國立台灣大學土木工程學系
Department of Civil Engineering, National Taiwan University, Taipei, Taiwan, R.O.C.
〔2〕 國立中興大學水土保持學系
Department of Soil and Water Conservation, National Chung Hsing University, Taichung, Taiwan, R.O.C.
* Corresponding Author. E-mail: cyhung@nchu.edu.tw
的影響其壩體長度較長，而崩塌所產生的堰塞湖壩體由於崩 塌瞬間進入河道，使其壩體長度較短。兩者在潰決的行為上 也因為這樣的特性而有所區別，在土石流型堰塞湖潰決歷程 中，堰塞湖體的水位將隨著壩體升高而升高，因此潰決的時 間將與支流土石流停止供應材料的時間相同，同時由於壩體 的長度較長，產生滲流破壞的可能性較低，多半以溢頂破壞 的行為為主。而在崩塌形成的壩體上，由於壩體形成的速度 較快，導致水位會需要慢慢蓄滿壩後空間在潰決，同時因崩 塌所產生的壩體因為長度較短，湖水滲流出的時間較短，因 此在蓄滿的過程中極有可能由於滲流的關係導致壩體產生 滲流破壞，在不同的破壞歷程中所產生的洪水歷線也將大不 相同。在過去有許多研究針對堰塞湖潰決進行探討，不論是 以實驗進行探討 (Cao et al., 2011; Chang et al., 2011; Zhou et al., 2013; )，或是以數值模式進行模擬，包含 HEC-RAS(Xiong, 2011; Butt et al., 2013)、BREACH(Fread,1988,1989)或是以偏 微分方程式計算潰壩的模式 (Liu et al, 2012; Wu et al.,2012)，
這些模式可以進行相當精細的模擬，但其參數的選擇都相當 複雜，並且在參數的決定上多以逆推法進行，並且在計算上 也需要耗費許多時間。因此本研究將以Capart(2013) 所提出 的潰壩模式為基礎，考量到台灣河川坡陡水急的特性，進行 堰塞會潰決的演算。
為了能有效地評估出堰塞湖潰決對傳遞至下游的影響，
本研究採用蒙地卡羅法進行未來100 年的事件推估，進行多 次堰塞湖潰決的模擬，在這樣大量的計算考慮下，本研究採 用一維非均勻河道運動波模型來進行演算。和一般常使用之 二維動力波模式 (如 HEC-ras) 相比，一維運動波模式不需 要詳細的地形資料，並可在小量的運算下得到良好之大尺度 洪水傳遞模擬結果 (Hunt, 1995; Kazezyilmaz-Alhan and Me- dina, 2007; Bohorquez, 2010; Capart 2013; Chen and Capart, 2017; Stilmant et al., 2018)。根據 Chen and Capart (2017)，一
維非均勻河道運動波模型除無法應用於回水之情況，以及無 法反映彎道內外側水位差異，但能有效可以反映河道寬度、
坡度變化對洪水傳遞速度及洪峰削減情形之影響，計算出準 確之洪水水位，符合本研究需求。
二、研究區域
荖濃溪為高屏溪之主流，1999-2017 年年平均雨量 3320 mm(高中、梅山、天池、六龜雨量站平均)，於 2009 年莫拉 克颱風更有超過3000 mm 的總降雨量。其中莫拉克颱風期 間，高達1500 mm/day 的高降雨強度使荖濃溪眾多支流上游 產生大量的崩塌 (圖 1)。其中崩塌面積最大之支流集水區為 布唐布納斯溪流域，崩塌面積高達500 公頃，佔整個支流集 水區面積之75 %，為全台崩塌比例最大之支流集水區之一。
布唐布納斯溪上游大量土砂在莫拉克颱風時以土石流 之形式傳遞進入主流，大幅擴大與荖濃溪匯流口處之土石流 沖積扇的規模。此沖積扇於2007 年曾經發生過阻塞主流形 成堰塞湖的案例，在莫拉克颱風時也有很高之堰塞湖發生可 能性。此外在莫拉克風災過後，大量的土砂材料並沒有完全 被傳輸進入主流中，而是在支流布唐布納斯溪上游堆積 (圖 1(b))，其後每當有颱風侵襲台灣或是大豪雨發生時，土砂材 料就被以土石流的形式傳輸進入荖濃溪主流堆積。
根據本研究現地調查及衛星影像得知，布唐布納斯溪於 2010、2011、2012、2013、2017、2018 年皆有大型土石流事 件發生，有新增之沖積扇堆積。而在2010、2011、2016、2017、
2018、2019 年皆有衛星影像或現地觀測記錄到此沖積扇阻 塞主流河道形成堰塞湖。若在有堰塞湖發生的時期發生暴雨 洪水事件，暴雨洪水疊加堰塞壩之潰決洪水，將會提高洪水 流量，可能對主流下游聚落產生嚴重的威脅。
圖1 研究區域與對象。(a) 研究區域位置與衛星影像；(b) 主支流匯流口堰塞壩與堰塞湖示意圖
Fig.1 Study area and target. (a) Location and satellite image of the study area; (b) Tributarydammed lake
布唐布納斯溪下游20 公里內荖濃溪沿岸有勤和、桃源、
寶來等主要聚落 (圖 1(a))，其中以勤和聚落受洪水災害之風 險最大。勤和聚落隸屬高雄市桃源區，位於布唐布納斯溪與 荖濃溪匯流口下游3.2 公里處的左岸河階地上，有省道台 20 線（南橫公路）通過。勤和聚落原距離荖濃溪河床有45 公 尺，然2009 年莫拉克颱風重創勤和聚落，並大幅提升荖濃 溪河道高程達30 公尺。近年來本團隊持續進行長期河道測 量，可見荖濃溪河道仍有逐漸淤高之趨勢，根據2018 年地 形測量之結果，勤和主要聚落與河床的高差僅餘8-12 公尺。
而連接勤和聚落與桃源聚落之撒拉阿塢橋位於台20 線 94k+720 處，橫越荖濃溪主流河道，有極高之受災風險。舊 橋於莫拉克颱風時毀損，並於2011 年重建。新撒拉阿塢橋 為一上坡型橋梁，橋梁中點之橋面底部標高560 公尺。其下 方河道高程於2011 年為 545 公尺，於 2019 年 1 月量測高程 為550 公尺，於 2019 年 6 月則淤高至 555 公尺，距離橋面 底部僅於5 公尺。在 2019 年 5 月 18-20 日大豪雨事件，荖 濃溪河水暴漲，公路局在此期間進行了兩次預警性封路，封 閉台20 線勤和至上游復興聚落間 5 公里路段。由以上資料 可見，撒拉阿塢橋為荖濃溪流域十分關鍵之高危區域，若有 堰塞湖潰壩產生之大洪水事件，此連通聚落之間的交通要道 之安全性將受到極大威脅，故本研究將著重探討各情境下此 區域之洪水水位。
三、研究方法
1. 降雨逕流模式
在過去的研究中，已發展大量降雨逕流模式 (USACE, 2004; Chang and Huang, 2013; Wu et al., 2014 )。在本研究中 考量到若採用過度複雜的模式，其參數的選擇較多，且須進 行較為大量的運算，同時，本研究希望能將結合降雨逕流模 式、堰塞湖潰壩模式及運動波模式進行演算，因此本研究採 用Criss and Wiston (2003) 針對山區河川所提出的擴散方程 式單位歷線模型進行分析。該方程式以達西定理為基礎，將 集水區簡化為一土壤模型，假設於上游降雨後透過擴散作用，
最後進入河道，其單位歷線模型可寫為
1 32 b
Q b e t
Vtotal b t
(1)
其中b 為一待率定的時間常數，e 為自然指數，而 Vtotal 為 總逕流量，在本研究中使用合理化公式計算其總逕流量進行 分析。
Vtotal CAw (2) P
其中 為累積總雨量， AP w 為上游集水區面積， C 為逕流 係數。透過歷史資料進行參數率定，可建立本研究區域之降 雨逕流模式。
2. 堰塞湖潰壩模式
在過去的研究中，山區河川泥砂運動的所造成的地形地 貌變化可將以單純一維的質量守恆方程式 (Exner Equation) 進行考慮 (Paola and Voller, 2005; Capart, 2013)，透過流動方 向上均勻寬度之底床高程隨時間的變化率與輸砂率在空間 上的變化率相同進行考慮，而本研究將堰塞湖壩體生成至潰 決以相同方程式進行模擬 (圖 2)
( ) zc J 0 B t t x
(3)
其中，x是指沿河道方向，t 是指時間， ( )B t 代表河道有效 寬度，其值將隨時間變化，在本研究中假設在壩體生成的過 程中有效寬度為河道的完整寬度BV 而在壩體下切潰決時有 效寬度僅為下切寬度BC ， zc 代表堰塞湖流動渠道底部之高 程，J 代表主流之輸砂率。在輸砂率的部分，過去的研究中 指出輸砂率與流量、坡度成正比指數關係 (Lane, 1955; Ca- part et al, 2010; Capart, 2013)，在本研究中採用最簡化版的輸 砂公式
( ) zc J KQ Sd d KQ td x
(4)
其中，K 為河川的擴散係數，其值將隨不同河川而改變，
Qd為往下游的流量，而Sd為坡度，其坡度可將其寫為壩體 高程對空間方向的微分 (zc/x)，將式 (3) 與式 (4) 合併 後，令D t( )KQd( )/ ( )t B t 將其簡化，可得
( ) 2 0 2 zc D t zc
t x
(5)
圖2 土石流型堰塞湖潰決模型與參數
Fig.2 The breaching model and parameters of the tributarydammed lake
同時，可以發現此類擴散方程式的形式較為複雜，由於係數 ( )
D t 也會隨時間而改變，求解的過程將會變得相當複雜，因 此為了讓方程式簡化，可使用變數轉換得到較為簡化的方程 式。將累積流量的概念導入方程式中，可令
( ) ( ) d KQ td
dt B t
(6)
其中， 代表的物理意義近似於單寬的累積流量，透過變數 轉換後，其方程式將由變係數的的擴散方程式可簡化為
2 2
( ) 2 2 0 zc D t zc zc zc
t x x
(7)
倘若重新思考該方程式所代表的物理意義，可以發現除了透 過變數轉換得到比較簡化版本的方程式外，此方程式在時間 軸上的定義也更具有物理意義，相較於先前的方程式，可想 像目前的方程式是一個隨流量變化而有不同變化程度的方 程式，不再是僅考量時間。而這樣的方程式也更為合理，由 於潰壩的過程本應隨流量變化而變化，不應是僅隨時間變化 而變化。
若需求解此方程式，邊界條件與初始條件是必須的。在本 研究中，可依據累積流量的變化將邊界條件寫為
( )dzc(0, ) KQ td dx t ID
(8)
其中，整體邊界條件所代表的物理意義為在主支流匯流口處 (x ) 之輸砂率完全來自於支流所提供的輸砂率0 ID。在本 研究中假設主流的輸砂量在堰塞湖開始形成時會淤積在堰 塞湖的上游端，而到匯流口時主要的砂量均完全來自於支流。
在支流的土砂供給上，本研究假設其土砂量體與主流上游流 量 (Qu) 成正比關係，寫為
D u
I Q (10)
其中為一無因次待率定參數。而在初始條件的部分，可寫 為
( ,0) 0
z xc z S xd (11)
其中，z0為匯流口的初始高程，而整體河床的高程應與坡 度相關。因此透過初始條件與邊界條件的設立，使用擴散方 程式之基礎解與邊界條件，可得解
2 ( )
0 2
zc z S xd Sd ierfc x
(12)
其中erfc 為誤差函數。此方程式可描述壩體整體隨時間空間
之變化，但其中累積流量變化仍尚未解出，因此簡化其條件，
僅針對潰口部份 (x ) 進行討論，並將其對時間進行微分0 寫為
22 2 ( )
dhc d KSd Q td
dt Sd dt B hc (13) 為求此方程式之解，本研究納入堰塞湖水體質量守恆方程式
dVL Q Qdt u- d (14)
其中，VL 為堰塞湖中水體的體積。又湖水的體積可以簡單 的三角形進行近似，因此可將湖水體積寫為
( + )2 2Bv
VL Su hc (15)
其中，本研究將 定義為水深即為水面高程與渠底之差值 (zwzc)，將湖水體積帶入湖水質量守恆方程式後，可寫成
2 2 ( ) ( ) Q Q
d Su u d KSd Q td dt Bv hc B t hc
(16)
假設通過潰口時之流量滿足寬頂堰之堰流公式 278
Qd BC g (17)
最後，將式(13)與式(17)同時進行聯立數值解，可得到壩體隨 時間、空間變化之歷程。
3. 一維非均勻河道運動波洪水傳遞
本研究採用一維非均勻河道運動波模型來模擬堰塞湖 潰壩洪水向下游之傳遞情形。運動波模式考慮連續方程式
A Q 0 t x
(18)
其中A x t 為通水斷面積， ( , )( , ) Q x t 為流量。並假設流況可 視為擬均勻流 (Quasi-uniform)，各河道區段內，重力水頭與 摩 擦 水 頭 達 成 平 衡 (Lighthill and Whitham, 1955;
Whitham,1974; Hunt, 1984; Kazezyilmaz-Alhan and Medina, 2007; Capart, 2013)。因此動量方程式可簡化為均勻流公式，
每個河段流量Q 與洪水斷面積 A 有單一之關係。
( , )
Q Q A x (19)
此關係可由該區段之河道斷面形狀、糙度、坡度，配合曼寧 公式進行計算。若河道為寬矩形斷面，則
0( ) 3/ 2
( , ) ( ) ( , )
f
Q x t gS x B x h x t
C (20)
其中h x t( , )A B/ v為水深，而C x 為糙度係數，f( ) B x V( ) 為河道寬度，S x0( ) 為底床坡度，此三者可隨河道位置 x 變 化。以有限體積法進行運動波洪水傳遞之數值計算，將第 (18) 式改寫為離散之形式
1 t * 1/2 * 1/2
k k k k
Ai Ai x Qi Q i
(21)
其中， 為空間網格大小， tx 為計算之時間步階， kAi 為 第i 個網格在第 k 個時間點時之通水斷面積， * 1/2k
Q i 為此 時間下第i 和 i+1 個網格之間之平均通量。而根據 Caleffi et al. (2003), Valiani et al. (2002), Zhang et al. (2016), and Lai and Khan (2018)，通量 * 1/2k
Q i 可以由 HLL 法 (Harten, Lax and van Leer, 1983)計算。
( )
1 1
* 1/2
k k
R R L k k
S Q S QL S S A A R
k i i i i
Q i S SL
(22)
其中 kQi 為第 i 個網格在第 k 個時間點時之平均流量，可由 對應之通水斷面積 kAi 和該網格所屬的河道斷面之流量率定 關係 (式 (19)) 進行計算。不同網格可以代入不同之流量率 定關係，以呈現非均勻之河道斷面形狀、寬度、糙度和坡度。
而SR 和 SL 分別代表此位置向下游及上游傳遞之最大波 速
max( ,i i 1,0), min( , i 1,0)
SR SL i (23) 而 i 為第 i 個網格之波速。由於運動波只向下游傳遞，波
速恆為正值，故而SL恆為0。因此，通量 * 1/2k
Q i 之計算可 以簡化為
* 1/2k k
Q i Qi (24)
而上游邊界之第一個通量 * 1/2k
Q 則可代入 t=tk時之潰壩流 量，作為邊界條件。
*
1/ 2k d( )k
Q Q t (25)
其中Q 為堰塞湖潰壩流量。此運動波數值計算模型對於上d
游邊界之輸入流量歷線沒有任何形式之限制，適用於由堰塞 湖潰壩而產生之多洪峰之流量歷線。配合河道各位置之初始 流量假設作為初始條件，即可進行洪水傳遞之演算模擬。為 了確保模擬過程之數值穩定，模擬所使用之時間步階t必 須滿足 Courant-Friedrichs-Lewy (CFT) 條件 (Harten et al., 1983)。
max 1 Cr t
x
(26)
其中Cr 為 Courant 數，max為所有網格中最大之波速。
max( ~ ) max
max 0
Q
i n A
i i Qik
(27)
這樣的一維洪水演算模型非常單純，雖然無法反應出回水效 應等較複雜的動力反應，但對於計算在大尺度的情況下洪水 向下游的傳遞相當有效 (Hunt, 1995; Singh, 2001)。透過不同 位置採用不同之流量率定關係，此模型可以應用於非均勻河 道，直接反應出不同河段河道形狀、寬度、糙度、坡度變化 導致之洪水傳遞速度差異和水深差異。
四、模擬參數率定
在過去不同的研究中，根據不同的模式所需率定的參數 及其數量均不相同，在本研究使用三種不同的模式進行堰塞 湖潰決的模擬推估，考量到實際運用的可行性，本研究採用 參數較少之模式進行模擬。在本研究所使用的模型中，除了 參數較少的優點外，參數的選擇也盡可能使用從觀測資料所 得到的結果進行率定分析，本節將針對各模式所需使用之參 數進行說明，並於本節最後列出在本模擬中之總參數表。
1. 降雨逕流模式
在本模式中有兩項參數需要率定，分別為 (1) 集流時間 b 及 (2) 集水區之逕流係數，針對集水區之逕流係數，本研 究中採用「水土保持技術規範」(2016) 之建議，考量本次模 擬區位為山區的情境下進行考量。而集流時間的推估，本研 究則使用水利署於 2004~2007 於興輝大橋所實測之流量資 料與同時段之高中氣象站之雨量資料進行分析，透過最小二 乘法的計算，進行集水區集流時間之單一參數率定，其結果 顯示於圖3，圖中可見透過本模式所計算出之流量資料與實 測之流量資料相當吻合。
圖3 降雨逕流模式驗證。(a) 高中氣象站之降雨資料；(b) 計算之流量資料與實測之流量資料比較 (Capart et al., 2010) Fig.3 Validation of rainfall-runoff model. (a) Precipitation data; (b) Theoretical predicted and measured runoff data
2. 堰塞湖潰決模式
在堰塞湖潰決模式中以擴散係數K 較難率定，根據過 去研究顯示該值之分布範圍甚廣，因應不同的河川特性會有 不同對應之值，但仍可透過現地的觀測進行參數的率定。本 研究參考Capart et al. (2010) 中根據現地資料所率定出擴散 係數，該研究同以荖濃溪布唐布納斯溪為研究區域，並且透 過過去堰塞湖實際觀察之案例進行分析，所得之參數供本研 究使用。除了擴散係數外，支流的強度仍需要進行計算，本 研究假設土石流之量體與其區間內之累積雨量成正比關係，
透過多年的地形觀測結果與光達掃描之資料估算其土石流 量體 (表 1)，而累積雨量的部分，根據水土保持局所發布的 土石流警戒值之最小值作為門檻，以高中氣象站中在該時間 區段內超過200mm 事件之雨量進行累加，將其土石量體與 累積雨量之關係繪製於圖4。
表1 土砂量體估算與累積雨量估算結果
Table 1 Estimation results of debris flow volume and accumulated precipitation
時間區間 累積雨量
[mm] 降雨時間
[day] 土石流堆積體積 [104 m3] 2007.08.07-2008.02 2636 13 98a 2010.07.17-2010.08.01 399 3 80b 2010.02.02-2010.08.01 957 11 180a 2010.08.01-2011.03 424 2 11a 2010.08.01-2014.01 5974 39 1880b
2014.01-2017.04 4365 28 790b 2017.04-2018.08 1506 14 980b
a以進入主流之土砂量 (由兩期主流河道縱剖面計算) 計算土石流堆積體積
b以支流集水區減少土砂量 (兩期 Lidar DEM 相減) 計算土石流堆積體積
由於在莫拉克風災後，大量土砂堆積於支流峽谷中，相 較於莫拉克風災前土石流之輸送趨勢資料，在圖中可以發現 風災前的土砂量體與累積雨量比值較風災後的比值為小。因 此在後續的計算上將不納入莫拉克風災前之資料。為建立土 石流量與累積雨量之關係，本研究採用 Takahashi(1991) 所 提出來的平衡濃度公式來估算土石流之量體V D
=(1 )( 1)(tan ) ST
VD nD s ST VR (28)
其中，S 為支流坡度，T n 為土石流之孔隙率，D s 為土石 流材料之比重，tan為土石流材料的內摩擦係數，最後累積 雨量所帶來的體積V ，本研究以合理化公式推估在支流集R 水區內之量體，相關參數均列於表2 中。Takahashi(1991) 所 提出的公式建立了土石流量體與水量之線性關係，參考水土 保持手冊使用之參數將其關係繪製於圖4 中。從圖中可以發 現資料分布趨勢呈現線性關係，但使用之參數非現地材料參 數，導致低估其量體大小。因此，本研究基於Takahashi(1991) 公式的線性關係基礎上，以線性的方式將實測資料進行迴歸，
推估土石流量體與水量之比例關係，最後將其參數納入潰壩 模式中計算使用。
圖4 累積雨量與土石流量體關係
Fig.4 Relationship between accumulated precipitation and debris flow volume
表2 降雨逕流模式與堰塞湖潰決模式參數總表
Table 2 Parameters of the rainfall-runoff model and the dammed lake breaching model
參數 物理意義 值 單位 資料來源
降雨逕流模式
b 集流時間 27 hr 依據2004-2007實測水文資料進行率定計算。
C 集水區逕流係數 0.75 [] 參考「水土保持技術規範」(2016) AW 主流上游集水區面積 542 km2 根據數值地形資料進行計算判定。
堰塞湖潰決模式(主流)
Su 初始上游坡度 0.0095 m/m
根據數值地形資料進行計算判定。
Sd 初始下游坡度 0.021 m/m
BV 溪谷寬度 300 m
BC 出流河道寬度 60 m
K 河川擴散係數 0.04 [] 參考Capart et al., 2010模式設定使用
支流土石流比例係數 6.9 10 3 [] 根據過去現地地形資料與水文資料進行計算 支流材料相關參數
s 土石流材料比重 2.6 []
參考「水土保持手冊:土石流篇」(2017)案例說明使用 材料參數進行設定
tan 材料內摩擦係數 0.7 []
nD 土石流材料孔隙率 0.35 m3/m3
ST 支流坡度 0.10 m/m 根據數值地形資料進行計算判定。
3. 一維非均勻河道運動波洪水傳遞模式
本 模 式 所 需 之 參 數 僅 有 河 道 寬 度B xV( )、 底 床 坡 度
0( )
S x 和糙度係數C 三者。本研究採用圖 5 中 Cross-section f 01- Cross-section 10 等 10 個斷面來簡化河道，將布唐布納斯 沖積扇致勤和聚落之間4 公里之河道分為 10 段。根據 Chen and Capart (2017)，在有限長度的河道中，假設河道參數區段 線性變化與假設河道參數區段均勻所得之洪水傳遞模擬結 果會相近。且在此案例中每段河道之河道參數變異性小，故 可假設每段之河段寬度與底床坡度均為均勻之情況，由對應 之斷面的河道參數來代表進行計算。本研究由Lidar 數值地 形求得各斷面之平均寬度和平均坡度，如表3 所列。而在糙 度係數C x 之部分，本研究假設此區域河道糙度係數為常f( ) 數，採用Capart et al. (2010)中根據現地興輝大橋 2005-2007
年流量─水深資料所率定出之糙度係數值 (Cf = 0.261)。
表3 河道斷面參數表
Table 3 Parameters of the river cross-sections 河道斷面 河道寬度BV[m] 底床坡度S0[-]
01 161.7 0.0162 02 213.8 0.0174 03 219.1 0.0176 04 239.9 0.0212 05 241.6 0.0224 06 283.6 0.0236 07 350.7 0.0249 08 356.8 0.0261 09 497.1 0.0287 10 661.2 0.0299
圖5 河道與斷面位置
Fig.5 Valley and cross-sections
五、研究情境
本研究以蒙地卡羅法進行未來一百年的事件模擬，將降 雨事件區分為 (1) 何時發生 (2) 發生事件之強度，針對事 件發生的時間本研究假設其機率分布符合帕松分布 (Pois- son distribution)。在過去的研究中 (Rodrı́guez-Iturbe 1988;
Dunn, 2004;)，帕松分布適合於描述單位時間內隨機事件發 生的次數的機率分布。透過收集氣象站1970~2017 年的雨量 資料 (高中、梅山、天池、六龜雨量站平均)，分析每個月共 會發生幾次事件 (pulse)，其分布狀況如圖 6(a) 所示，其中 發生事件數之定義如氣象局所公告之降雨事件定義相同進 行分析。而在事件的強度方面，根據過去的研究 (Kurothe et al., 1997; Goel et al., 2000)，本研究假設其強度分布符合指數 分布(exponential distribution)，同樣的利用氣象站的資料，本 研究找出其每月各事件的降雨量之平均值 (圖 6(b))，可透過 將事件之平均降雨量與發生事件之機率相乘進行後，可得到 平均每年的降雨量進行簡單的驗證。根據本研究所收集之資 料計算結果值約為3000 mm/year，與目前所觀測的結果大致 符合。透過發生事件的次數與其降雨量之抽選，可進行未來
100 年的降雨資料模擬。
同時，在本研究也將氣候變遷的因子納入考慮，氣候變 遷主要影響氣溫的變化，而氣溫的變化也影響了大氣中的水 氣含量，因此本研究透過雨量的變化考量氣候變遷對未來模 擬的影響。根據「臺灣氣候變遷推估與資訊平台建置報告」
(2010)，該報告指出透過過去 30 年的降雨資料分析，可計算 出每10 年平均總雨量呈現 85.25 mm 之增長幅度，因此在本 研究中假設在事件發生的分布不變的前提下，將增加每月的 每事件之平均降雨量進行計算，當模擬年達100 年時，其每 月的每事件平均降雨量將增加如圖6(b) 中所示。
待建立出完整之雨量變化資料後，將其雨量資料透過以 建置完成之降雨逕流模式進行流量之模擬。本研究在未考量 氣候變遷與考量氣候變遷的情境下分別進行20 次 100 年的 模擬，其模擬出之上游累積流量如圖6(c) 所示。在圖中可以 明顯發現在考量氣候變遷的情境下，累積流量在模擬初期並 不會有太大的改變，但由於隨時間增加每事件之平均降雨量 也會隨之增加，因此反應在累積流量相對未考慮氣候變遷的 案例下而隨時間增加而有所改變。最後，本研究將以兩個情 境搭配潰壩反應所產生的結果進行討論。
圖6 蒙地卡羅取樣降雨事件。(a)平均各月事件發生次數； (b) 平均各月事件發生強度；(c)累積流量歷線
Fig.6 Precipitation simulated by Monte Carlo method. (a) Average pulse number in each month; (b) Average pulse inten- sity in each month; (c) Cumulative discharge along the 100 years
六、結果與討論
本節將先針對堰塞湖在單一事件潰決之結果進行討論，
透過本模式可得到壩體自生成至潰決的完整歷程(圖 7(a))，
在考慮土石材料輸送主支流交會處的歷程上，相較於崩塌型 堰塞湖的計算中一次性阻斷河道的行為，本模式中土砂材料 不完全一次性阻斷主流，而是保有渠道寬度 (BC) 供主流流 動，同時也帶走部分土砂材料，但壩體仍會隨支流土砂輸入 而抬升，直到支流停止提供材料後即開始潰決行為。在形貌 上與崩塌型堰塞湖不同，在下游壩址處並不會有顯著的坡度 變化，而是與原河道地形相切。本模式中下游流量的計算採 用堰流公式進行，其流量大小將取決於潰口水深的變化，圖
7(b) 中記錄在模式計算中潰口高程與水面高程的變化情況，
整體歷程水深於潰決時達最大值，對應的流量變化繪於圖 7(c)。考慮台灣河川的特性，由於山區河川的坡度較陡，導 致土石流型堰塞湖的湖體較小，若僅有考慮堰塞湖體的水體 潰決所造成的流量，相較於主流的流量其量體甚小，因此若 是單純考慮堰塞湖體之潰決，並不會對下游村落造成太大的 威脅，但若是考慮在颱風期間，主流流量甚大時堰塞湖潰決，
此時潰決的流量將與主流流量以非線性的方式疊加，此時流 量將會對下游村落等造成相當大的威脅。在本研究中以最危 險的狀況進行計算，即假設堰塞湖開始潰決時間與主流之洪 峰抵達時間相同進行計算。
圖7 堰塞湖潰決過程。(a)堰塞壩與渠道剖面高程變化； (b)堰塞壩體與外流渠道底部高度變化；(c)堰塞壩上游入流流量歷線與 下游出流流流歷線
Fig.7 Breaching process of tributary-dammed lake. (a) Evolution of debris dam and outflow channel elevations; (b) Evo- lution of debris dam and channel heights; (c) Upstream and downstream hydrographs of tributary-dammed lake
本研究將未考量氣候變遷與考量氣候變遷的情境下分 別進行20 次 100 年之洪水模擬，得到 20 組 100 年間多次 (每組 100 年約 1000 場，合計約 20000 場) 洪水事件之流量 歷線。考量洪水事件會伴隨堰塞湖之生成與潰決，將這些流 量歷線作為堰塞壩壩體上游流量歷線，代入本研究之堰塞湖 潰決模式，模擬支流與主流匯流口堰塞壩形成與潰決之過程。
本研究篩選出考量氣候變遷與否之兩種情境每次模擬中100 年間堰塞壩高度超越 7.5 m 之事件 (每組 100 年模擬間約 350 場事件，合計約 7000 場事件)，將其堰塞壩下游潰壩流 量歷線帶入一維非均勻河道運動波模型，進行潰壩洪水向下
游傳遞之洪水演算。為了對比未考量堰塞湖之情境，本研究 亦將未考慮堰塞湖的洪水事件流量歷線代入一維非均勻河 道運動波模型進行洪水演算，總結模擬情境與條件列於表4。
根據運動波模型模擬之結果，可以得到每種情境各場事件中 目標斷面之最大洪水水位。本研究採用極端值序列法計算每 次模擬之兩百年重現期之洪水水位，再將每種情境的20 次 模擬進行平均，並求得20 次模擬結果之間的標準差。本研 究著重於勤和聚落周邊三個斷面 (Cross-section 07-09) 及撒 拉阿塢橋斷面 (Cross-section 10)，勤和聚落周邊三個斷面模 擬結果如圖8 所示，撒拉阿塢橋斷面模擬結果如圖 9 所示。
圖8 勤和聚落洪水模擬結果。(a)(c)(e) 斷面形狀與洪水水位；(b)(d)(f) 平均最大洪水深累積機率分布；(a)(b) 斷面 7，(c)(d) 斷 面8，(e)(f) 斷面 9
Fig.8 Maximum flooding elevations at Ching-Ho village. (a)(c)(e) Cross-section shape and maximum flooding elevations;
(b)(d)(f) Accumulated probability of the average maximum flooding depths; (a)(b) Cross-section 7; (c)(d) Cross-sec- tion 8; (e)(f) Cross-section 9
圖9 撒拉阿塢橋洪水模擬結果。(a)斷面形狀、橋梁位置與洪水水位；(b)平均最大洪水深累積機率分布
Fig.9 Maximum flooding elevations at Sala-a-wu Bridge. (a) Cross-section shape and maximum flooding elevations; (b) Accumulated probability of the average maximum flooding depths
在勤和部落 (圖 8) 與撒拉阿塢橋 (圖 9) 的斷面上，都 可以發現在先不考慮潰壩的情境下，氣候變遷的僅產生約10
%的水位變化，若是同在氣候變遷的情境下，潰壩對水位的 影響則會高達30 %；整體來說，以一百年的模擬情境下，氣 候變遷造成的水位變化較少，但本研究是以一簡化的情境，
單純考慮平均雨量的增加，在事件的發生機率上並沒有進行 調整，若是機率也隨之上升，其結果之影響將會更大。反觀 在同一情境 (Scenario1 &2) 下，潰壩事件造成水位的變化更 大，考量到該區域由於布唐布納斯溪的活躍，該地區發生堰 塞湖之機率甚大，在進行該地區的防災規劃上應該也將堰塞 湖的影響一併納入考慮，避免錯估水位的變化。
針對勤和部落，由於在莫拉克風災後，大量土砂抬升了 河道，導致河道寬度增加，反應在模擬的情境下水位都不甚 高，不會有直接淹沒到聚落之危險。而在2017 年 6 月 3-4 日 豪雨事件中，洪水水位亦在遠低於勤和聚落路面高程的情況 下，掏空了該聚落11 戶民宅地基和周邊道路地基，使路面 和屋舍陷落於河中被沖走。由此可見，雖然河水水位沒有直 接淹沒路面，但地下水位之上漲和洪水水位對河岸之衝擊仍 會威脅路基和住宅。
而在撒拉阿烏橋的部分 (圖 9)，由於該橋橫跨荖濃溪，
因此其受到的影響相較於勤和部落更大，在兩百年的重現期 水位影響下，即使是未考量氣候變遷與潰壩事件的情境下 (Scenario 4)，水位也會超過目前設計的橋面高度造成危險。
除了洪水的影響外，根據現地調查和地形Lidar 測量之結果 可知，荖濃溪主流河道底床高程變動劇烈，且近年多呈持續 淤高之情況，然本研究之洪水水位僅以 2018 年河道底床高 程為基準進行計算，並未考慮河道淤高。若將河道淤高之情 況考量進入，本區域之洪水水位亦會有明顯之提高，將對道 路橋樑及聚落產生更高之危害。整體來說，不論是河床變動 的情況或是堰塞湖潰決的水位影響，都應該在未來的防災規 劃中納入考慮。
七、結 論
本研究採用了較簡易的模式計算土石流型堰塞湖潰決 時對下游的影響。考量到參數率定與現地資料收集的困難度，
本研究中的模式將參數的數量減到最小，並且多以較容易取 得之資料進行參數率定，也將洪水傳遞模型中對地形資料的 需求降到最低。在簡易模型的計算下，可使用蒙地卡羅法進 行大量的計算推估未來的趨勢，本研究也將氣候變遷的因素 納入模型中進行計算。經過四種情境的模擬後，本研究提出 在山區支流活躍的情況下，在工程構造物與防災規劃上，相 較於氣候變遷的因素，更應將支流活躍所造成的堰塞湖潰決 情境納入考量中。此外，透過本研究所收集的現地測量資料，
也可得知目前由於支流堆積的大量土砂材料導致河床快速 的抬升，因此除了洪水高度的影響外，在未來也應將河床的 高度變動性一併納入考慮，可進行較為完整的防災評估。
參考文獻
[1] Bohorquez, P. (2010). “Competition between kinematic and dynamic waves in floods on steep slopes.” Journal of Fluid Mechanics. 645, 375-409. <https://doi.org/10.1017/S0022 11200999276X>
[2] Butt, M.J., Umar, M., and Qamar, R. (2013). “Landslide dam and subsequent dam-break flood estimation using HEC-RAS model in Northern Pakistan.” Natural Hazards 65, 241. <https://doi.org/10.1007/s11069-012-0361-8>
[3] Caleffi, V., Valiani, A., and Zanni, A. (2003). “Finite vol- ume method for simulating extreme flood events in natural channels.” Journal of Hydraulic Research. 41, 167-177.
<https://doi.org /10.1080/00221680309499959>
[4] Cao, Z., Yue, Z., and Pender, G. (2011). “Landslide dam failure and flood hydraulics. Part I: experimental investiga- tion.” Natural Hazards, 59, 1003. <https://doi.org/10.1007/
s11069-011-9814-8>
[5] Capart, H., Hsu, J.P.C., Lai, S.Y.J., and Hsieh, M. (2010).
“Formation and decay of a tributarydammed lake, Laonong
River, Taiwan.” Water Resources Research, 46. <https://doi.
org/10. 1029/2010WR009159>
[6] Capart, H. (2013). “Analytical solutions for gradual dam breaching and downstream river flooding”. Water Re- sources Research, 49. 1968-1987. <https://doi.org/10.1002/
wrcr.20167>
[7] Chang, D. S., Zhang, L. M., Xu, Y., Huang, R.Q., (2011).
“Field testing of erodibility of two landslide dams triggered by the 12 May Wenchuan earthquake.” Landslides 8(3), 321-332.
[8] Chang, C., and Huang, W. (2013). “Hydrological modeling of typhoon-induced extreme storm runoffs from Shihmen watershed to reservoir, Taiwan.” Natural Hazards, 67, 747.
<https://doi.org/10.1007/s11069-013-0600-7>
[9] Chen T.Y., and Capart H. (2017). “Kinematic wave propa- gation in variable-width channels applied to the Tangjiashan dam breach flood.” Cross-strait Symposium on Soil and Water Conservation, 2017 August 21-23th, Cheng- du, China.
[10] Criss, R.E., and Winston, W.E. (2003). “Hydrograph for small basins following intense storms.” Geophysical Re- search Letters, 30(6), 1314-1318.
[11] Dunn P.K. (2004). “Occurrence and quantity of precipita- tion can be modeled simultaneously”. International Journal of Climatology, 24, 1231-1239
[12] Fread, D.L. (1988). “BREACH: A Breach Erosion Model for Earthen Dams.” Report, National Weather Service, NOAA, Silver Spring, Md.
[13] Fread, D.L. (1998). Dam-Breach Modeling and Flood Rout- ing: A Perspective on Present Capabilities and Future Di- rections, the U.S. Department of Agriculture, Stillwater, Okla. the U.S. Department of Agriculture, Stillwater, Okla.
[14] Goel, N.K., Kurothe, R.S., Mathur, B.S., and Vogel, R.M.
(2000). “A derived flood frequency distribution for corre- lated rainfallintensity and duration.” Journal of Hydrology, 228, 56-67.
[15] Harten, A., Lax, P.D., and Leer, B.V. (1983). “On upstream differencing and Godunov-type schemes for hyperbolic conservation laws.” SIAM Review, 25, 35-61.
[16] Hunt, B. (1984). “Dam-break solution.” Journal of Hydrau- lic Engineering, 110, 675-686.
[17] Hunt, B. (1995). Fluid mechanics for civil engineers.
Christchurch, New Zealand: Department of Civil Engineer- ing, University of Canterbury.
[18] Kazezyılmaz-Alhan, C.M., and Medina, M.A. (2007).
“Kinematic and Diffusion Waves: Analytical and Numeri- cal Solutions to Overland and Channel Flow.” Journal of Hydraulic Engineering, 133, 217-228. <https://doi.org /10.1061/(asce)0733-9429(2007)133:2(217)>
[19] Kedem, B., Pavlopoulos, H., Guan, X., and Short, D.A.
(1994). “Probability distribution model for rain rate.” Jour- nal of Applied Meteorology and Climatology, 33, 1486- 1493, doi:10.1175/1520-0450 (1994)033,1486:APDMFR.
2.0.CO;2.
[20] Kurothe, R.S., Goel, N.K., and Mathur, B.S. (1997). “De- rived flood frequency distribution for correlated rainfall in- tensity and duration.” Water Resources Research, 33 (9), 2103-2107.
[21] Lai, W., and Khan, A.A. (2018). “Numerical solution of the Saint-Venant equations by an efficient hybrid finite-vol- ume/finite-difference method.” International Journal of Hydrogen Energy, 30, 189-202. <https://doi.org/10.1007 /s42241-018-0020-y>
[22] Lane, E. W. (1955). “The importance of fluvial morphology in hydraulic engineering.” Proceedings of the American So- ciety of Civil Engineers, 1955, 81(7), 1-17.
[23] Lighthill, M.J., and Whitham, G., (1955). “On kinematic waves I. Flood movement in long rivers.” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 229, 281-316. <https://doi.org/10.1098 /rspa.1955.0088>
[24] Liu, F., Fu, X., Wang, G., and Duan, J. (2012). “Physically based simulation of dam breach development for Tangjiashan Quake Dam, China” Environmental Earth Sci- ences. 65, 1081-1094.
[25] Paola, C., and Voller, V.R. (2005). “A generalized Exner equation for sediment mass balance.” Journal of Geophysi- cal Research, 110, F04014, doi: 10.1029/2004JF000274.
[26] Rodrı́guez-Iturbe, I., Cox, D.R, and Isham, V. (1988). “A point process model for rainfall: further developments.”
Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 417, 283-298.
[27] Singh, V.P. (2001). “Kinematic wave modelling in water resources: A historical perspective.” Hydrological Pro- cesses, 15, 671–706. <https://doi.org/10.1002/hyp.99>
[28] Stilmant, F., Pirotton, M., Archambeau, P., Erpicum, S., and Dewals, B. (2017). “Hydraulic Determination of Dam Re- leases to Generate Warning Waves in a Mountain Stream:
Performance of an Analytical Kinematic Wave Model.”
Journal of Hydraulic Engineering, 144, 05017006. <https://
doi.org/10.1061/(asce) hy.1943-7900.0001425>
[29] Takahashi T. (1991). Debris flow, IAHR monograph.
Balkema Publication, Rotterdam.
[30] USACE (2004). Upper Mississippi River System Flow Fre- quency Study:Final report
[31] Valiani, A., Caleffi, V., and Zanni, A. (2002). “Case Study:
Malpasset Dam-Break Simulation using a Two-Dimen- sional Finite Volume Method.” Journal of Hydraulic Engi- neering,128, 460-472. <https://doi.org/10.1061/(asce)0733 -9429(2002)128:5 (460)>
[32] Whitham, G.B. (1974). Linear and Nonlinear Waves.
Wiley-Interscience, New York.
[33] Wu, W., Marsooli, R., and He Z.G. (2012). “Depth-aver- aged two-dimensional model of unsteady flow and sediment transport due to noncohesive embankment break/breaching”
Journal of Hydraulic Engineering, 138, 503-516.
[34] Wu, M.C., Lin, G.F., and Lin, H.Y. (2014) “Improving the forecasts of extreme streamflow by support vector regres- sion with the data extracted by self-organizing map” Hydro- logical Processes, 28(2), 386-397.
[35] Xiong Y. (2011). “A dam break analysis using HEC-RAS”.
Journal of Water Resource and Protection, 3, 370-379.
[36] Zhang, M.L., Xu, Y.Y., Qiao, Y., Jiang, H.Z., Zhang, Z.Z., and Zhang, G.S. (2016). “Numerical simulation of flow and bed morphology in the case of dam break floods with vege- tation effect.” Journal of Hydrodynamics, 28, 23-32.
<https://doi.org/10. 1016/S1001-6058 (16)60604-2>
[37] Zhou, G.G.D., Cui, P., Chen, H.Y., Zhu, X.H., Tang, J.B., and Sun, Q.C. (2013). “Experimental study on cascading landslide dam failures by upstream flows” Landslides, 10, 633. <https://doi.org/10.1007/s10346-012-0352-6>
[38] 科技部 (2010) 「臺灣氣候變遷推估與資訊平台建置」， 科技部補助專題研究計畫成果報告。Ministry of Science
and Technology (2010) “Taiwan Climate Change Projec- tion and Information Platform Project” Published by Minis- try of Science and Technology, Taiwan (In Chinese).
[39] 行政院農業委員會 (2016) 「水土保持技術規範」。
Council of Agriculture (2016) “Technical Regulations for Soil and Water Conservation” Published by Council of Ag- riculture, Taiwan (In Chinese).
2019 年 09 月 09 日 收稿 2019 年 10 月 01 日 修正 2019 年 10 月 30 日 接受 (本文開放討論至 2020 年 03 月 31 日)
![Table 3 Parameters of the river cross-sections 河道斷面 河道寬度B V [m] 底床坡度S 0 [-] 01 161.7 0.0162 02 213.8 0.0174 03 219.1 0.0176 04 239.9 0.0212 05 241.6 0.0224 06 283.6 0.0236 07 350.7 0.0249 08 356.8 0.0261 09 497.1 0.0287 10 661.2 0.02](https://thumb-ap.123doks.com/thumbv2/9libinfo/9122226.406877/7.892.472.820.174.388/Table3Parametersoftherivercrosssections河道斷面河道寬度BVm底床坡度S1161716.webp)
