平原河網地區水庫潰壩洪水數值模擬分析
——以江西省柘林水庫為例

吳家陽，劉佳明，單定軍，張飛峰

(1.長江勘測規劃設計研究有限責任公司，湖北 武漢 430010； 2.國網江西省電力有限公司 柘林水電廠，江西 九江 332000)

0 引 言

中國水能資源集中分布的西南地區河道以山區型為主，河道斷面呈深“V”形。分布于上述地區的大壩一旦潰決，潰壩洪水被限制在河道兩側山體之間，可采用一維水動力學模型模擬潰壩洪水的傳播過程。范鴻杰[1]以雅礱江官地、二灘水電站為研究對象，采用一維水動力學模型耦合DAMBRK潰壩模型，計算潰壩洪水對下游沿岸居民點和涉河建筑物的影響；方崇惠等[2]以清江招徠河水庫為研究對象，采用Hec-Ras軟件構建了一維潰壩洪水演進模型，并分析了潰壩洪水對兩岸風險點和隔河巖水庫的影響；羅琳等[3]采用MIKE 11軟件，計算了瀘州白節河玉水水庫不同潰決形式、潰決歷時方案下的潰壩洪水；張大偉等[4]基于自主研發的GIS平臺，開發了一維潰壩洪水分析系統，并成功應用于法國Malpasset拱壩和湖北咸寧市鵝公頸水庫的潰決計算中；周興波等[5]基于MIKE 11軟件，采用一維模型計算了唐家山堰塞湖潰決洪水流量過程及潰口的沖刷過程；趙雪瑩等[6]采用MIKE 11搭建了貴陽松柏山、花溪梯級水庫一維潰壩洪水模型，分析了梯級水庫群不同水位組合潰壩洪水對下游的影響。

對于平原河網地區的大壩潰決，潰壩洪水往往超過下游河道安全泄量，漫溢堤防后進入地勢寬廣和平坦的洪泛區，此時一維模型不再適用，必須采用高維模型計算潰壩洪水在洪泛區的傳播過程。王雯等[7]基于不可壓縮的N-S方程建立了潰壩水流運動的二維數學模型，模擬了漢江石泉水庫的潰壩洪水傳播過程；潘存鴻等[8]從Boltzmann方程出發，應用有限體積法建立了基于無結構三角形網格的二維潰壩波整體數值模型，并應用于浙江省某水庫潰壩后的洪水傳播過程；葛杭建等[9]采用MIKE 21中的二維水動力學模型，模擬了浙江省亭下水庫潰壩洪水演進過程；王曉玲等[10]建立了耦合VOF法的三維紊流數學模型，模擬并繪制了河北省東武仕水庫潰壩洪水風險圖。雖然高維模型在潰壩洪水演進計算中取得了一定進展，但仍存在計算量較大、計算效率偏低，以及不便于設置閘、堰、泵站等水工建筑物等不足。

為了充分發揮一維、高維水動力學模型各自的優點，模擬平原河網地區大壩潰決洪水及其在河道及洪泛區的傳播過程，本文提出了一維-一二維耦合水動力學模型，分別采用一維-二維模型模擬潰壩洪水在河道內和洪泛區的演進過程。為證明模型的有效性，以江西省柘林水利樞紐為例，模擬了不同潰決工況下潰決洪水在修河尾閭地區的傳播過程，并統計了潰壩洪水的淹沒風險。

1 流域及工程概況

1.1 流域概況

1.1.1 自然地理特征

修河又名修水，為江西省五大河流之一，屬于鄱陽湖水系，流域水系詳見圖1。修河干流全長約419 km，流域面積約14 797 km2。修河干流抱子石以上為上游山區，抱子石至柘林大壩為中游丘陵區，柘林水庫以下為下游段，自虬津以下進入濱湖平原地區，圩堤縱橫，河道交錯，河道平均比降僅為0.12‰。

圖1 修河流域水系Fig.1 River network of Xiuhe Basin

1.1.2 修河下游圩垸概況

經過長期建設，修河下游河段現已初步形成以圩垸為基礎，配合水庫的防洪工程體系。在易受洪水威脅的沿河兩岸臺階地和沖積平原，大多已初步形成由圩堤圈成的規模不等、防洪標準不同的防護保護圈。根據2021年江西省堤防工程基本情況統計數據，修河下游河段共分布有4處重點圩和5處一般圩。其中，4處重點圩分別為三角聯圩、九合聯圩、永北聯圩、郭東聯圩；5處一般圩分別為高橋圩、幸福圩、朝陽圩、馬口圩和立新圩。各圩垸基本情況見表1。

表1 修河下游河段圩垸基本情況

1.2 樞紐工程概況

柘林水利樞紐工程位于江西省西北部修河中游的永修縣柘林鎮附近，為修河干流中游大(1) 型水利樞紐工程。壩址以上流域面積9 340 km2，總庫容79.2億m3，為江西省庫容最大的蓄水工程，也是亞洲庫容最大的土壩水庫。樞紐工程按千年一遇洪水設計、可能最大洪水校核。主要建筑物包括主壩、3座副壩、2座溢洪道、泄空洞、電站廠房、船筏道等。主壩為黏土和混凝土心墻土石壩，壩頂高程73.5 m，最大壩高63.5 m，壩頂長590.75 m。3座副壩均為土壩，Ⅰ，Ⅱ，Ⅲ副壩最高壩高分別為20.7，3.0 m和18.4 m。柘林水庫為多年調節水庫，正常蓄水位65.00 m，主汛期(4～6月)防洪限制水位64.00 m，后汛期(7～9月)防洪限制水位65.00 m，死水位50.00 m。大壩千年一遇設計洪水位70.13 m，遭遇可能最大洪水最高調洪水位73.01 m。

2 潰壩洪水計算模型構建

考慮柘林水庫庫容、壩高均較大，潰壩洪峰流量勢必超過修河下游河道安全泄量，潰壩洪水漫溢堤防后進入河道兩側泛洪區。為準確模擬柘林水利樞紐潰壩洪水及其下游演進過程，本文基于丹麥水利研究所開發的MIKE FLOOD平臺，耦合MIKE 11河道洪水演進一維水動力學模塊和MIKE 21地表漫流二維水動力學模塊，搭建柘林水利樞紐潰壩洪水演進一維-二維耦合模型。其中，潰壩洪水在修河、潦河等河道內的演進過程采用一維模型計算，潰壩洪水漫過堤防后在修河下游、鄱陽湖區圩垸內及低洼地區的傳播過程采用二維模型計算。

2.1 模型計算范圍及地形資料

2.1.1 一維模型計算范圍及河道斷面

一維模型計算范圍包括長江干流九江-大通河段、鄱陽湖湖區及鄱陽湖五河尾閭河段。其中，長江干流九江-大通河段，總長179.7 km，斷面總數59個，斷面平均間距3.09 km；修河干流柘林水利樞紐主壩以下河段總長78.3 km，斷面總數34個，斷面平均間距2.4 km；潦河干流萬家埠以下河段總長30.1 km，斷面總數16個，斷面平均間距2.01 km；贛江干流外洲以下河段，包括北支、中支和南支，總長159.5 km，斷面總數37個，斷面平均間距4.43 km；撫河干流李家渡以下河段，包括北支和南支，總長120.5 km，斷面總數23個，斷面平均間距5.48 km；信江干流梅港以下河段，包括東支和西支，總長96.1 km，斷面總數22個，斷面平均間距4.58 km；樂安河干流虎山以下河段，總長38.7 km，斷面總數7個，斷面平均間距6.45 km；昌江干流渡峰坑以下河段，總長70.9 km，斷面總數17個，斷面平均間距4.43 km。一維模型研究范圍見圖2。

圖2 一維模型計算范圍示意Fig.2 Studying area of 1D model

2.1.2 二維模型計算范圍及河道斷面

二維模型計算范圍被修河、潦河分為3塊子區域，詳見圖3。其中，修河左岸區域，南以修河干流白槎-吳城河段為邊界，北以50 m等高線為邊界，東以蚌湖為邊界，總面積870.61 km2；潦河左岸區域，東以潦河干流萬家埠-潦河河口河段為邊界，北以修河干流虬津-潦河河口段為邊界，西、南以50 m 等高線為邊界，總面積285.69 km2；潦河右岸區域，西以潦河干流萬家埠-潦河河口河段為邊界，北以修河干流潦河-吳城河段為邊界，東以贛江北支象山-吳城河段為邊界，南以50 m等高線為邊界，總面積626.16 km2。

圖3 二維模型計算范圍示意Fig.3 Studying area of 2D model

2.2 一維、二維模型耦合方式

本次采用側向耦合方式進行一維-一二維模型的耦合求解，耦合方式詳見圖4。側向耦合允許一維模型從河道兩側連接到二維模型，并采用傳統堰流公式或用戶自定義的泄流能力曲線計算通過側向連接點的流量。該種耦合方式非常適用于洪水漫溢河道堤防進入河道外側洪泛區的場景。當河道內洪水位未超過兩側堤防堤頂高程時，僅一維模型執行計算任務，二維模型不參與計算；當河道內洪水位高出兩岸堤防堤頂漫溢進入洪泛區后，采用二維模型計算洪水在洪泛區的演進過程。

圖4 一維、二維模型側向耦合方式示意Fig.4 Schematics of lateral link for the 1D-2D coupling

當河道內側洪水位高出兩岸堤防堤頂高程后，本文假設堤防漫而不潰，并采用寬頂堰的堰流公式計算漫溢堤防進入洪泛區的單寬流量。

(1)

式中：Q為漫溢堤防的單寬流量，m2/s；σs為淹沒系數；m為自由溢流的流量系數；g為重力加速度，本文取9.81 m/s2；H為堰上水位，m；h為堤頂高程，m。計算所得流量將以點源方式加入二維模型中。

2.3 邊界條件及初始條件

2.3.1 邊界條件

模型邊界條件包括柘林水庫入庫流量邊界條件、長江干流和鄱陽湖五河上游入流邊界條件、鄱陽湖區間入流邊界條件以及長江干流下游水位流量關系邊界條件。根據后續擬定的潰壩洪水計算方案，柘林水庫入庫邊界條件采用柘林水庫不同頻率P入庫設計洪水過程、可能最大洪水過程；長江干流上游邊界采用九江站的多年平均流量或枯水期多年平均流量；鄱陽湖五河分別采用各尾閭控制站的多年平均流量或枯水期多年平均流量；鄱陽湖區間入流采用多年平均流量或枯水期多年平均流量。柘林水庫不同頻率設計洪水過程線和可能最大洪水過程線詳見圖5和圖6，設計洪水過程采用1954年實際洪水按照洪峰和洪量控制進行縮放。模型入流邊界條件詳見表2。

圖5 柘林水庫不同頻率設計洪水過程線Fig.5 Designed flood hydrograph of different return period of Zhelin Reservoir

圖6 柘林水庫可能最大洪水過程線Fig.6 Hydrograph of probable maximum flood of Zhelin Reservoir

表2 模型入流邊界條件取值

2.3.2 初始條件

模型初始條件包括初始流量和初始水位等，是潰壩計算開始時的初始狀態。在柘林水利樞紐潰壩前，將根據長江干流、鄱陽湖五河來水，以及鄱陽湖區間入流情況等運行程序直至達到恒定流條件，并以此作為潰壩模擬的初始條件。

2.4 河道及洪泛區糙率選取

2.4.1 河道糙率選取

長江干流九江-大通河段糙率分別采用1998年和2016年湖口站、大通站實測水文資料進行率定和驗證，糙率取值范圍為0.022～0.032。其中，九江-安慶糙率取值為0.022，安慶-大通糙率取值為0.032。鄱陽湖五河尾閭地區河道糙率取值參考文獻[11]中的經驗糙率表，通過衛星影像資料和現場查勘情況，綜合考慮河床組成及床面特性、河道平面形態、河道岸壁特性等因素，對于順直河段的平灘河槽內糙率取0.020～0.026，彎曲且有茂密植被覆蓋的河段糙率取0.040～0.065。

2.4.2 洪泛區糙率選取

對于河道以外洪泛區糙率，由于缺乏實測洪水及淹沒情況資料，無法對其取值進行率定和驗證，本次參考文獻[11]中河道灘地糙率表，綜合洪泛區地形、下墊面情況等因素，糙率取值范圍為0.030～0.085。

2.5 土石壩逐漸潰決模型

根據相關規范，土石壩的潰決形式一般為逐漸潰決，由漫頂和滲透破壞引起。本文采用中國應用較廣泛的BREACH模型[12-14]模擬柘林水利樞紐主、副壩的潰決過程。BREACH模型由美國國家氣象局的Fread[15]提出，可用于勻質土壩、堆石壩和天然堰塞壩的漫頂破壞和管涌破壞這2種形式的潰壩過程計算。其中，潰口的底部沖刷下切過程采用Meyer-Peter & Müller推移質輸沙公式模擬，潰口兩側的擴展速率取為底部下切速率的80%。對于土石壩的漫頂逐漸潰決，其潰口底部沖刷率可表示為剪應力的函數：

(2)

式中：hd為潰口底部高程，m；Φ(τ)為剪應力τ的函數。Chen等[16]提出雙曲線形式的侵蝕率模型，解決了指數形式侵蝕率模型中泄流渠流速較大時計算不收斂的問題，其表達式為

(3)

式中：k為在剪應力τ范圍內允許dhd/dt接近其極值的單位變換因子；τc為臨界剪應力，Pa；a為τ=τc時，雙曲線斜率的倒數；b為dhd/dt極值的倒數。Chen等[16]在反演唐家山堰塞湖實際潰決過程時得到：k=100，a=1.1，b=0.007。通過曼寧公式，可以得到剪應力與流速之間的關系：

(4)

式中：v為水流流速，m/s；γ為水的容重，N/m3；R為水力半徑，m；J為水力坡降；n為泄流渠糙率系數。

3 潰壩洪水計算結果分析

3.1 計算方案

柘林水利樞紐主、副壩均為土石壩，其潰決形式均按逐漸潰決考慮。根據相關技術規范并參考國內外類似土石壩潰決實例，引發潰決的誘因重點考慮遭遇超標準入庫洪水漫頂潰決、管涌滲透破壞潰決等2種情形，形成2種潰壩洪水計算方案如下。

3.1.1 遭遇超標準入庫洪水漫頂潰決方案

對于遭遇超標準入庫洪水漫頂潰決情形，入庫邊界條件采用柘林水利樞紐遭遇鋒面雨的可能最大洪水過程(洪峰流量25 600 m3/s，詳見圖6)，初始庫水位采用后汛期防洪限制水位(65 m)。對于水庫運用方式，考慮第二溢洪道無法開啟，樞紐第一溢洪道和泄空洞全開，下泄流量按照樞紐下泄能力曲線確定。考慮到Ⅱ副壩最大壩高僅3 m，潰壩洪水量級及淹沒范圍遠小于主壩及其他副壩，因此本方案僅考慮主壩及Ⅰ，Ⅲ副壩的潰決。

3.1.2 管涌滲透破壞潰決方案

考慮到管涌滲透破壞在汛期和枯水期均可能發生，本次從偏不利角度出發，假定在汛期發生，入庫邊界條件采用柘林水庫1 000 a一遇設計洪水過程線，初始庫水位為后汛期防洪限制水位?？紤]到主、副壩同時發生管涌滲透破壞的概率極小，本文僅考慮主壩潰決情形。柘林水利樞紐主壩一旦發生管涌滲透破壞，水庫開啟全部泄洪設施盡快降低庫水位。

3.2 結果分析

潰壩洪水計算結果包括壩址流量過程、壩址下游沿程流量及水位過程，以及壩址下游各典型斷面的洪峰流量、最高水位等。為便于展示及分析潰壩洪水的下游傳播過程，本文沿修河下游河段共選取11個典型斷面，如圖7所示。

圖7 修河下游河段典型斷面分布Fig.7 Typical section distribution in the downstream reaches of Xiuhe River

3.2.1 沿程洪峰流量過程

不同計算方案柘林水利樞紐壩下各典型斷面洪峰流量過程和洪峰流量詳見表3和圖8。相較于管涌滲透破壞潰決方案，漫頂潰決方案入庫洪水更大、壩體潰決歷時更短、起潰水位更高，其壩址洪峰流量更大，洪峰流量過程曲線更為尖瘦。主壩及Ⅰ，Ⅲ副壩同時發生漫頂潰決，壩址洪峰流量87 328 m3/s，而管涌滲透破壞潰決的壩址洪峰流量為46 173 m3/s。

柘林水庫庫容較大，潰壩洪水峰高量大，加之修河干流柘林水利樞紐壩址-吳城鎮段總長僅為78.35 km，潰壩洪水傳播至吳城鎮后削峰效果并不明顯。例如漫頂潰決方案，潰壩洪水傳播至吳城時仍有79 608 m3/s，相較于壩址洪峰流量僅坦化了7 720 m3/s。

表3 不同潰壩方案沿程洪峰流量

3.2.2 沿程水位漲落過程

不同計算方案柘林水利樞紐壩下各典型斷面水位漲落過程和最高水位詳見表4和圖9。對于超標準洪水漫頂潰決方案，壩址洪峰流量更大，壩體沖刷更為強烈，壩下游各斷面水位漲落過程更為陡峭，由于入庫可能最大洪水過程為典型的雙峰形態，且潰決洪水與第一個入庫洪峰遭遇，導致本方案壩下游各斷面水位過程也呈現一定的雙峰形態。對于管涌滲透破壞潰決方案，由于其壩體沖刷強度較漫頂潰決情形小，沖刷過程持續時間長，導致其壩下游各斷面沿程水位漲落過程更為平緩。

圖8 沿程洪峰流量過程Fig.8 Peak discharge profiles along Xiuhe River

表4 不同潰壩方案沿程最高水位與最大水位漲幅

從各方案沿程水位漲幅可見，受修河干流柘林水利樞紐壩址-虬津段地形控制，洪水基本被約束在河道兩側山頭之間，水位漲幅較大，潰壩洪水出虬津后進入較為寬廣的濱湖平原地區，潰壩洪水可充分利用河道兩側平原灘地行洪，加之修河尾閭地區圩垸對潰壩洪水的分蓄作用，水位漲幅有較明顯的下降。例如漫頂潰決方案，柘林水利樞紐壩址-虬津段最大水位漲幅為18.14～18.36 m，虬津以下河段最大水位漲幅為3.73～14.07 m，潰壩洪水出虬津后最大水位漲幅驟降約4.10 m。

4 潰壩洪水風險分析

由于柘林水庫庫容較大，潰壩洪水經修河和河道兩側泛洪區進入鄱陽湖區后，會引起湖區水位出現不同程度的抬升，勢必導致五河尾閭地區泛濫成災。本次僅統計九江市永修縣、廬山市、共青城市，以及南昌市建新區等受災最嚴重地區的淹沒風險數據。

圖9 沿程水位變化過程Fig.9 Water levels along Xiuhe River

本文利用高分辨率遙感衛星影像，采用目視及機器解譯相結合的方式獲取淹沒區域的房屋、耕地、道路和橋梁等地物信息，根據各類地物的頂、底部高程與其所處位置的最高洪水位關系判斷受淹情況。其中，對于房屋，若最高洪水位高于房屋底部高程，則認為房屋被淹沒；對于橋梁，若最高洪水位高于橋面高程，則認為橋梁被淹沒；對于道路，若最高洪水位高于路面高程，則認為道路被淹沒；對于碼頭，若最高洪水位高于碼頭平臺高程，則認為碼頭被淹沒。對于影響人口的統計，由于難以獲取淹沒區域人口密度的準確分布情況，本次采用江西省2020年發布的人均住房建筑面積，結合淹沒房屋面積統計各影響鄉鎮的影響人口數量，結果見圖10～11。

圖10 超標準洪水漫頂潰決方案淹沒水深Fig.10 Flooded area of dam break flood in overtopping breach failure scenario

圖11 管涌滲透破壞潰決方案淹沒水深Fig.11 Flooded area of dam break flood in piping seepage failure scenario

柘林水利樞紐主壩及Ⅰ，Ⅲ副壩遭遇入庫可能最大洪水漫頂潰決方案，淹沒范圍詳見圖10。修河下游地區受淹面積約632.88 km2，共涉及九江市永修縣、南昌市新建區、九江市廬山市和九江市共青城市等4個區縣、23個鄉鎮街道。其中，淹沒房屋總面積6.6 km2，淹沒農田耕地總面積298.97 km2，受淹橋梁共19座，受淹高速公路總長28.91 km，受淹國道與省道總長93.77 km，受淹鐵路總長57.99 km，1座火車站和3座碼頭受影響，受影響總人口31.25萬人。

柘林水利樞紐主壩發生管涌滲透破壞潰決方案，淹沒范圍詳見圖11。修河下游地區總受淹面積為486.19 km2，共涉及4個區縣(九江市永修縣、南昌市新建區、九江市廬山市和九江市共青城市)、20個鄉鎮街道，其中淹沒房屋總面積3.16 km2，淹沒農田耕地總面積213.14 km2，受淹橋梁13座，受淹高速公路總長10.79 km，受淹國道與省道總長73.99 km，受淹鐵路31.22 km，1座碼頭受影響，受影響總人口14.84萬人。

可見，無論遭遇可能最大洪水漫頂潰決還是管涌滲透破壞潰決方案，由于修河下游地區和鄱陽湖區圩垸較多，潰壩洪水漫堤進入三角聯圩、永北聯圩、九合聯圩、郭東聯圩等圩垸后，洪水位迅速降低，加之鄱陽湖對潰壩洪水的蓄納作用，潰壩洪水仍可被限制在修河下游地區，不至影響南昌市、共青城市、德安縣等重要城市的主城區。

5 結 論

本文為模擬柘林水利樞紐工程潰壩洪水及其下游傳播過程，構建了潰壩洪水一維-一二維耦合水動力學模型，分別采用一維-一二維模型模擬潰壩洪水在河道內和洪泛區的演進過程，克服了一維模型不適用于平原河網地區潰壩洪水模擬，以及高維模擬計算效率偏低和水工建筑物設置不便等不足，所得主要結論如下。

(1) 柘林水利樞紐主、副壩均為土石壩，其潰壩形式應為逐漸潰決，致潰誘因重點考慮遭遇超標準入庫洪水漫頂潰決和管涌滲透破壞潰決等2種情形。

(2) 柘林水庫遭遇可能最大洪水，主壩及Ⅰ，Ⅲ副壩同時發生漫頂潰決方案，壩址洪峰流量約87 000 m3/s，壩下柘林鎮最大水位漲幅達18.14 m，洪水傳播至吳城鎮洪峰仍有約79 000 m3/s，最大水位漲幅達3.73 m。柘林水利樞紐主壩發生管涌滲透破壞潰決方案，壩址洪峰流量約46 000 m3/s，壩下柘林鎮最大水位漲幅達13.23 m，洪水傳播至吳城鎮洪峰仍有近41 000 m3/s，最大水位漲幅達1.34 m。

(3) 修河干流柘林水利樞紐主壩壩址-虬津段兩側有山頭控制，潰壩洪水雖然漫溢堤防，但仍被約束在河道兩側山地之間，水位漲幅較大。虬津以下河段為地形開闊的濱湖平原地區，潰壩洪水一旦漫溢河道，兩側堤防可一瀉千里，水位漲幅明顯下降。

(4) 由于修河下游地區和鄱陽湖區圩垸較多，潰壩洪水漫堤進入三角聯圩、永北聯圩、九合聯圩、郭東聯圩等圩垸后，洪水位迅速降低，加之鄱陽湖對潰壩洪水的蓄納作用，潰壩洪水不至影響南昌市、共青城市、德安縣等重要城市的主城區。