沭西片防洪保護區洪水淹沒風險特征

賁 鵬，胡 勇

(1.安徽省(水利部淮河水利委員會)水利科學研究院，合肥 230000；2.水利水資源安徽省重點實驗室，合肥 230000)

0 引 言

沂沭泗流域地處我國南北氣候過渡帶，年降雨分布不均，多集中在汛期。由于歷史上黃河長期奪泗奪淮，沂沭泗河下游原有水系被破壞，洪水出路不暢，洪澇災害頻繁發生。流域上游大多為山區性河道，源短、坡陡，洪水來勢迅猛，峰高流急；下游河道則地勢平坦，河床淤積嚴重，行洪緩慢且持續時間長，極易發生洪澇災害[1,2]。

我國處于從控制洪水向洪水管理的轉變時期，區域洪水淹沒特征涉及洪澇模擬、洪水避難、災害預警、災情評估和公眾洪水風險意識等方面，是進行洪水管理的科學依據之一[3,4]，而水動力數值模擬是分析區域洪水淹沒特征的基礎支撐平臺。目前，利用各種水力學模型對洪水的數值模擬技術做了大量研究，苑希民[5]等對潰口、多源洪水耦合模型進行了研究，并運用在防洪保護區洪水風險計算中；李義天等[6]開展了淮河中游行蓄洪區開啟方式、口門位置及其分洪風險研究；虞邦義等[7,8]建立淮河中游水動力數學模型，研究了行洪區閘門開啟過程的不同階段對淮河干流上下游洪水演進的影響和行蓄洪區內洪水演進規律以及過流能力；劉衛林、金玲等對中小河流潰堤進行了模擬和風險分析[9,10]。這些研究廣泛運用于洪水演進、防洪排澇、仿真模擬、洪水調度等方面，而對防洪保護區的山區性洪水模擬、遭遇組合以及河道與保護區洪水互饋等綜合分析研究較少。

以沂沭泗河流域沭西片防洪保護區研究對象，建立河道一維、保護區二維以及耦合水動力數學模型，模擬了河道及保護區水流運動情況，以沂河馬頭潰口為典型案例，分析了潰堤洪水淹沒要素、風險點分布等信息，并對成果合理性進行了分析，模型可以較好地反映區域洪水演進情況，作為洪水淹沒要素以及洪水風險分析的計算平臺，提高了區域防洪減災能力。

1 研究方法

1.1 多源洪水耦合模擬

耦合模型以潰口或者堤防為媒介，補充一維模型和二維模型的水位或者流量關系，實現模型實時交互計算[11]，一維、二維模型連接斷面的水位和流量的關系如下：

Z1=Z2

(1)

(2)

式中：Z1、Z2分別為一、二維模型在連接斷面處的水位；Q1、Q2分別為一、二維模型在連接斷面處的流量；Uk為連接斷面法向流速；hk為水深；dk為距離；k為連接斷面的編號。

1.2 河道洪水演進模擬

一維水動力模型的控制方程為Saint Venant方程組。

連續方程：

(3)

動量方程：

(4)

式中：Q為流量；Z為水位；A為過水斷面的面積；B為水面寬度；q為旁側入流流量；K為流量模數。對上式采用Abbott六點隱格式進行離散求解。

1.3 保護區洪水演進模擬

對用Navier-Stokes方程沿水深進行積分，可得平面二維淺水水流控制方程。

連續性方程：

(5)

動量方程：

(6)

(7)

式中：h，ξ分別為水深和水位；t為時間；u，v分別為x，y方向的垂向平均流速；Ex，Ey分別為x方向和y方向的水流紊動黏性系數；τbx,τby為x方向和y方向的底部摩阻；τsx,τsy分別為風對自由表面x方向和y方向的剪切力；f=2ωsinφ為科氏力，φ為計算水域的地理緯度，q為源或匯的流量。為更好地適應復雜邊界，采用基于非結構網格的有限體積法對控制方程進行離散求解。

1.4 潰堤洪水模擬

潰口口門概化為寬頂堰，按照Villemonte公式計算其過流能力，主要影響因素為潰口寬度、頂部高程及其兩側的水位等。

(8)

式中：Q為流量；C為堰流系數；b為寬度；k為堰流指數；Hus為堰上游水位；Hds為堰下游水位；Hw為堰頂高程。

2 模型構建及驗證

2.1 研究區概況

沭西片防洪保護區屬于沂沭泗河流域，由沂河左堤、沭河右堤、分沂入沭水道右堤、駱馬湖東堤、新沂河北堤等保護區域，涉及江蘇省徐州市新沂市、邳州市及山東省臨沂市臨沭縣、郯城縣，面積約2 380 km2。保護區主要外部河流水系有沂河、沭河、分沂入沭水道、新沂河和駱馬湖以及保護區內的白馬河、墨河等河流，區域地形北高南低，被多條河流和道路分隔[12,13]，保護區水系和主要控制工程見圖1。

圖1 區域概化圖Fig.1 Regional generalization

2.2 河道一維模型

劉家道口至駱馬湖段河道，分沂入沭水道、沭河人民勝利堰至新沂河段河道、新沂河嶂山閘至沭河閘段河道等建立了一維水動力數學模型，并考慮沿程閘壩、橋梁等跨河建筑物。為了滿足計算時間和精度的要求，一維模型空間步長取500～1 000 m，時間步長為30 s。河道主槽糙率為0.021～0.03，河道灘地糙率為0.032～0.040。

2.3 保護區二維模型

2.3.1 模擬范圍

二維模擬計算范圍是沂河左堤、沭河右堤、分沂入沐右堤、新沂河左堤和駱馬湖東堤組合成的封閉區域以及駱馬湖湖區。

2.3.2 結構物概化及網格優化

考慮保護區內道路、堤防等高于地面0.5 m的線性阻水建筑物以及內河對于面上洪水演進過程的影響，將道路、堤防、水系河道作為網格剖分的內部約束條件處理。在ARCGIS軟件中，從基礎地理數據圖層中，分別提取保護區內阻水、導水地物要素。阻水要素主要包括：鐵路、高速、省國道、縣道等道路，白馬河、墨河、新戴運河等河道堤防；導水要素主要包括：橋梁和涵洞。通過現場調研對道路、堤防、橋梁、涵洞等信息進行復核，最終確定道路和堤防等阻水建筑物249條，模型中采用堤防概化；橋梁涵洞等導水建筑物136個，模型中概化為涵洞。

為了反映地形變化趨勢，模型網格邊長不超過300 m，網格面積不超過0.05 km2，區域高程范圍15～160 m，網格總數114 439 個，網格節點總數57 715 個。重要地區、地形變化較大部分的計算網格要適當加密。

2.3.3 分布式糙率參數

從基礎地理數據中分別提取道路、水系河流、居民地、植被等SHP圖層要素，利用ArcGIS軟件“相交”工具，分別與網格圖層“相交”計算，統計各網格中道路、水系河流、居民地、植被等不同下墊面元素面積，通過計算各網格中各圖層元素所占面積比，按照加權平均計算出綜合糙率值。

2.4 多源耦合模擬計算

沂河、沭河、新沂河和駱馬湖采用潰堤或漫溢形式與保護區連接，區域內的白馬河、墨河等重要河流采用二維網格或者一維河道漫溢模式，模型概化見圖2。

圖2 耦合模型概化示意圖Fig.2 Sketch of coupling model generalization

3 模型合理性分析

3.1 典型洪水過程復演

采用1993年、2003年和2012年實測洪水過程對一維河網模型進行率定與驗證。計算結果表明，閘壩和橋梁計算水位落差和設計落差相當，各測站實測水位與計算水位之差得絕對值≤20 cm，計算流量與實測流量的相對誤差≤10%，具有較高精度。由于模型率定與驗證的測站較多，僅以2003年洪水主要測站為例，說明驗證情況，沂河港上站水位和流量、沭河人民勝利堰壩下水位、新沂河嶂山閘下水位計算過程與實測過程比較見圖3至圖6。

圖3 港上站計算與實測水位過程Fig.3 Calculated and measured water level process of the ShangGang station

圖4 港上站計算與實測流量過程Fig.4 Calculated and measured discharge process of the ShangGang station

圖5 人民勝利堰壩下計算與實測水位過程Fig.5 Calculated and measured water level process of the Renminshengliyan station

圖6 嶂山閘下計算與實測流量過程Fig.6 Calculated and measured discharge process of the Renminshengliyan station

3.2 水量平衡分析

根據模型運算結果統計計算區域進出水量、河道槽蓄量與淹沒區積水量，基本滿足“保護區積水量=上游來水量-下游出水量-河道槽蓄量”關系，所建模型洪水計算滿足水量平衡要求。

3.3 流場及地形分析

洪水淹沒水深、洪水流速、流場分布均與區域地形起伏分布特征相匹配且存在一定規律性，即洪水由高向低流動、低洼地帶積水且水深較大、地形較高區域水深較小或無淹沒。二維模型能夠反映面上洪水演進過程和特征。

4 洪澇風險分析

4.1 典型潰堤洪水淹沒分析

保護區洪水威脅主要來源于外河潰堤，外河主要包括沂河、沭河、新沂河和駱馬湖。結合河勢地形、地質狀況、工程狀況、歷史出險等情況，綜合考慮潰口對保護區影響較大和各種不利情況的組合等，保護區周邊共設置了9處堤防潰口以及部分河段堤防漫溢。以沂河50年一遇洪水為例，對馬頭堤防潰決進行洪水風險分析。

沂河發生50年一遇洪水，左堤馬頭險工發生潰決，潰口寬度100 m，馬頭斷面達到最高水位時開始潰堤。結果表明，馬頭潰口進洪歷時共計13.2 h，洪峰流量537 m3/s，進洪量833 萬m3，少量洪水漫過白馬河，淹沒區總面積107.7 km2，沂河與白馬河交匯處西北處附近洼地淹沒水深達到1.9 m，其他區域淹沒水深較淺，一般在1.0 m以內。潰堤洪水起于馬頭鎮，自北向南演進，受縣道郯勝線、白馬河右堤的阻隔影響，洪水主要將沿白馬河右堤向下游演進，部分洪水在省道S232與白馬河交匯處漫過白馬河堤防向南演進。洪水主要淹沒區域包括郯城縣的郯城鎮、馬頭鎮、港上鎮、楊集鎮、花園鄉、歸昌鄉、銀杏產業開發區；邳州市的港上鎮；新沂市的合溝鎮、瓦窯鎮、港頭鎮等11個鄉鎮，淹沒水深見圖7。馬頭潰堤后，3 h洪水淹沒至白馬河堤防，24 h淹沒至花園鄉，48 h淹沒至白馬河入沂河處。洪水到達時間見圖8。

圖7 馬頭潰口洪水淹沒水深Fig.7 The Submergence depth of Matou breach

4.2 不同量級洪水淹沒分析

沂河馬頭潰口設置了100年一遇、50年一遇、20年一遇、1957年型、1974年型等5個年型洪水方案。按照洪水量級的排序為100年一遇>1974年型>1957年型>50年一遇>20年一遇。以淹沒水深差值作為疊加分析的依據，統計了0.05～0.5,0.5～1.0,1.0～1.5,1.5～2.5 m和>2.5 m等分級水深的淹沒面積，各方案的淹沒面積為76～210 km2，約為保護區面積的4.6%～12.8%，各分析方案淹沒面積變幅較大，詳見表1。淹沒面積排序為100年一遇>1974年型>1957年型>50年一遇>20年一遇，與洪量一致。通過合理性分析可知，模型

表1 沂河馬頭潰口淹沒水深和面積分級統計 km2Tab.1 lassification statistics of inundation depth and area of Matou breach in Yihe River

可以反映河道和保護區洪水演進和洪澇分布情況，作為區域洪澇計算平臺。

5 結 論

(1)基于沭西片防洪保護區地形地貌及其外河洪水特征，構建了實時動態耦合模型水動力數學模型。河道一維模型采用典型年實測洪水資料驗證，水位和流量誤差較小；保護區二維模型通過現場調研查勘設置了阻水導水構筑物、分布式糙率等，可以反映保護區的水流特性；模型參數設置合理，具有較高的精度，可以用于區域洪水淹沒特征分析計算。

(2)運用耦合模型，開展了沭西片防洪保護區洪澇模擬計算，潰口設置了100年一遇、50年一遇、20年一遇、1957年型、1974年型等5個洪水方案。以沂河50年一遇洪水碼頭潰口為例，潰口進洪歷時13.2 h，洪峰流量537 m3/s，進洪量833 萬m3，淹沒區總面積107.7 km2。洪水主要沿白馬河右堤向下游演進，少量洪水漫過白馬河，淹沒范圍涉及郯城鎮、馬頭鎮等11個鄉鎮。

(3)沂河馬頭潰口進洪歷時短、峰值流量大，防洪保護區內洪水演進速度快，大部分區域淹沒水深淺，洪水將沿著白馬河向下游演進，結合淹沒水深、洪水路徑和到達時間對人員和財產轉移進行轉移和安置，為搶險救災提供依據。

□