賁 鵬,胡 勇
(1.安徽省(水利部淮河水利委員會)水利科學(xué)研究院,合肥 230000;2.水利水資源安徽省重點實驗室,合肥 230000)
沂沭泗流域地處我國南北氣候過渡帶,年降雨分布不均,多集中在汛期。由于歷史上黃河長期奪泗奪淮,沂沭泗河下游原有水系被破壞,洪水出路不暢,洪澇災(zāi)害頻繁發(fā)生。流域上游大多為山區(qū)性河道,源短、坡陡,洪水來勢迅猛,峰高流急;下游河道則地勢平坦,河床淤積嚴(yán)重,行洪緩慢且持續(xù)時間長,極易發(fā)生洪澇災(zāi)害[1,2]。
我國處于從控制洪水向洪水管理的轉(zhuǎn)變時期,區(qū)域洪水淹沒特征涉及洪澇模擬、洪水避難、災(zāi)害預(yù)警、災(zāi)情評估和公眾洪水風(fēng)險意識等方面,是進行洪水管理的科學(xué)依據(jù)之一[3,4],而水動力數(shù)值模擬是分析區(qū)域洪水淹沒特征的基礎(chǔ)支撐平臺。目前,利用各種水力學(xué)模型對洪水的數(shù)值模擬技術(shù)做了大量研究,苑希民[5]等對潰口、多源洪水耦合模型進行了研究,并運用在防洪保護區(qū)洪水風(fēng)險計算中;李義天等[6]開展了淮河中游行蓄洪區(qū)開啟方式、口門位置及其分洪風(fēng)險研究;虞邦義等[7,8]建立淮河中游水動力數(shù)學(xué)模型,研究了行洪區(qū)閘門開啟過程的不同階段對淮河干流上下游洪水演進的影響和行蓄洪區(qū)內(nèi)洪水演進規(guī)律以及過流能力;劉衛(wèi)林、金玲等對中小河流潰堤進行了模擬和風(fēng)險分析[9,10]。這些研究廣泛運用于洪水演進、防洪排澇、仿真模擬、洪水調(diào)度等方面,而對防洪保護區(qū)的山區(qū)性洪水模擬、遭遇組合以及河道與保護區(qū)洪水互饋等綜合分析研究較少。
以沂沭泗河流域沭西片防洪保護區(qū)研究對象,建立河道一維、保護區(qū)二維以及耦合水動力數(shù)學(xué)模型,模擬了河道及保護區(qū)水流運動情況,以沂河馬頭潰口為典型案例,分析了潰堤洪水淹沒要素、風(fēng)險點分布等信息,并對成果合理性進行了分析,模型可以較好地反映區(qū)域洪水演進情況,作為洪水淹沒要素以及洪水風(fēng)險分析的計算平臺,提高了區(qū)域防洪減災(zāi)能力。
耦合模型以潰口或者堤防為媒介,補充一維模型和二維模型的水位或者流量關(guān)系,實現(xiàn)模型實時交互計算[11],一維、二維模型連接斷面的水位和流量的關(guān)系如下:
Z1=Z2
(1)
(2)
式中:Z1、Z2分別為一、二維模型在連接斷面處的水位;Q1、Q2分別為一、二維模型在連接斷面處的流量;Uk為連接斷面法向流速;hk為水深;dk為距離;k為連接斷面的編號。
一維水動力模型的控制方程為Saint Venant方程組。
連續(xù)方程:
(3)
動量方程:
(4)
式中:Q為流量;Z為水位;A為過水?dāng)嗝娴拿娣e;B為水面寬度;q為旁側(cè)入流流量;K為流量模數(shù)。對上式采用Abbott六點隱格式進行離散求解。
對用Navier-Stokes方程沿水深進行積分,可得平面二維淺水水流控制方程。
連續(xù)性方程:
(5)
動量方程:
(6)
(7)
式中:h,ξ分別為水深和水位;t為時間;u,v分別為x,y方向的垂向平均流速;Ex,Ey分別為x方向和y方向的水流紊動黏性系數(shù);τbx,τby為x方向和y方向的底部摩阻;τsx,τsy分別為風(fēng)對自由表面x方向和y方向的剪切力;f=2ωsinφ為科氏力,φ為計算水域的地理緯度,q為源或匯的流量。為更好地適應(yīng)復(fù)雜邊界,采用基于非結(jié)構(gòu)網(wǎng)格的有限體積法對控制方程進行離散求解。
潰口口門概化為寬頂堰,按照Villemonte公式計算其過流能力,主要影響因素為潰口寬度、頂部高程及其兩側(cè)的水位等。
(8)
式中:Q為流量;C為堰流系數(shù);b為寬度;k為堰流指數(shù);Hus為堰上游水位;Hds為堰下游水位;Hw為堰頂高程。
沭西片防洪保護區(qū)屬于沂沭泗河流域,由沂河左堤、沭河右堤、分沂入沭水道右堤、駱馬湖東堤、新沂河北堤等保護區(qū)域,涉及江蘇省徐州市新沂市、邳州市及山東省臨沂市臨沭縣、郯城縣,面積約2 380 km2。保護區(qū)主要外部河流水系有沂河、沭河、分沂入沭水道、新沂河和駱馬湖以及保護區(qū)內(nèi)的白馬河、墨河等河流,區(qū)域地形北高南低,被多條河流和道路分隔[12,13],保護區(qū)水系和主要控制工程見圖1。
圖1 區(qū)域概化圖Fig.1 Regional generalization
劉家道口至駱馬湖段河道,分沂入沭水道、沭河人民勝利堰至新沂河段河道、新沂河嶂山閘至沭河閘段河道等建立了一維水動力數(shù)學(xué)模型,并考慮沿程閘壩、橋梁等跨河建筑物。為了滿足計算時間和精度的要求,一維模型空間步長取500~1 000 m,時間步長為30 s。河道主槽糙率為0.021~0.03,河道灘地糙率為0.032~0.040。
2.3.1 模擬范圍
二維模擬計算范圍是沂河左堤、沭河右堤、分沂入沐右堤、新沂河左堤和駱馬湖東堤組合成的封閉區(qū)域以及駱馬湖湖區(qū)。
2.3.2 結(jié)構(gòu)物概化及網(wǎng)格優(yōu)化
考慮保護區(qū)內(nèi)道路、堤防等高于地面0.5 m的線性阻水建筑物以及內(nèi)河對于面上洪水演進過程的影響,將道路、堤防、水系河道作為網(wǎng)格剖分的內(nèi)部約束條件處理。在ARCGIS軟件中,從基礎(chǔ)地理數(shù)據(jù)圖層中,分別提取保護區(qū)內(nèi)阻水、導(dǎo)水地物要素。阻水要素主要包括:鐵路、高速、省國道、縣道等道路,白馬河、墨河、新戴運河等河道堤防;導(dǎo)水要素主要包括:橋梁和涵洞。通過現(xiàn)場調(diào)研對道路、堤防、橋梁、涵洞等信息進行復(fù)核,最終確定道路和堤防等阻水建筑物249條,模型中采用堤防概化;橋梁涵洞等導(dǎo)水建筑物136個,模型中概化為涵洞。
為了反映地形變化趨勢,模型網(wǎng)格邊長不超過300 m,網(wǎng)格面積不超過0.05 km2,區(qū)域高程范圍15~160 m,網(wǎng)格總數(shù)114 439 個,網(wǎng)格節(jié)點總數(shù)57 715 個。重要地區(qū)、地形變化較大部分的計算網(wǎng)格要適當(dāng)加密。
2.3.3 分布式糙率參數(shù)
從基礎(chǔ)地理數(shù)據(jù)中分別提取道路、水系河流、居民地、植被等SHP圖層要素,利用ArcGIS軟件“相交”工具,分別與網(wǎng)格圖層“相交”計算,統(tǒng)計各網(wǎng)格中道路、水系河流、居民地、植被等不同下墊面元素面積,通過計算各網(wǎng)格中各圖層元素所占面積比,按照加權(quán)平均計算出綜合糙率值。
沂河、沭河、新沂河和駱馬湖采用潰堤或漫溢形式與保護區(qū)連接,區(qū)域內(nèi)的白馬河、墨河等重要河流采用二維網(wǎng)格或者一維河道漫溢模式,模型概化見圖2。
圖2 耦合模型概化示意圖Fig.2 Sketch of coupling model generalization
采用1993年、2003年和2012年實測洪水過程對一維河網(wǎng)模型進行率定與驗證。計算結(jié)果表明,閘壩和橋梁計算水位落差和設(shè)計落差相當(dāng),各測站實測水位與計算水位之差得絕對值≤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
根據(jù)模型運算結(jié)果統(tǒng)計計算區(qū)域進出水量、河道槽蓄量與淹沒區(qū)積水量,基本滿足“保護區(qū)積水量=上游來水量-下游出水量-河道槽蓄量”關(guān)系,所建模型洪水計算滿足水量平衡要求。
洪水淹沒水深、洪水流速、流場分布均與區(qū)域地形起伏分布特征相匹配且存在一定規(guī)律性,即洪水由高向低流動、低洼地帶積水且水深較大、地形較高區(qū)域水深較小或無淹沒。二維模型能夠反映面上洪水演進過程和特征。
保護區(qū)洪水威脅主要來源于外河潰堤,外河主要包括沂河、沭河、新沂河和駱馬湖。結(jié)合河勢地形、地質(zhì)狀況、工程狀況、歷史出險等情況,綜合考慮潰口對保護區(qū)影響較大和各種不利情況的組合等,保護區(qū)周邊共設(shè)置了9處堤防潰口以及部分河段堤防漫溢。以沂河50年一遇洪水為例,對馬頭堤防潰決進行洪水風(fēng)險分析。
沂河發(fā)生50年一遇洪水,左堤馬頭險工發(fā)生潰決,潰口寬度100 m,馬頭斷面達(dá)到最高水位時開始潰堤。結(jié)果表明,馬頭潰口進洪歷時共計13.2 h,洪峰流量537 m3/s,進洪量833 萬m3,少量洪水漫過白馬河,淹沒區(qū)總面積107.7 km2,沂河與白馬河交匯處西北處附近洼地淹沒水深達(dá)到1.9 m,其他區(qū)域淹沒水深較淺,一般在1.0 m以內(nèi)。潰堤洪水起于馬頭鎮(zhèn),自北向南演進,受縣道郯勝線、白馬河右堤的阻隔影響,洪水主要將沿白馬河右堤向下游演進,部分洪水在省道S232與白馬河交匯處漫過白馬河堤防向南演進。洪水主要淹沒區(qū)域包括郯城縣的郯城鎮(zhèn)、馬頭鎮(zhèn)、港上鎮(zhèn)、楊集鎮(zhèn)、花園鄉(xiāng)、歸昌鄉(xiāng)、銀杏產(chǎn)業(yè)開發(fā)區(qū);邳州市的港上鎮(zhèn);新沂市的合溝鎮(zhèn)、瓦窯鎮(zhèn)、港頭鎮(zhèn)等11個鄉(xiāng)鎮(zhèn),淹沒水深見圖7。馬頭潰堤后,3 h洪水淹沒至白馬河堤防,24 h淹沒至花園鄉(xiāng),48 h淹沒至白馬河入沂河處。洪水到達(dá)時間見圖8。
圖7 馬頭潰口洪水淹沒水深Fig.7 The Submergence depth of Matou breach
沂河馬頭潰口設(shè)置了100年一遇、50年一遇、20年一遇、1957年型、1974年型等5個年型洪水方案。按照洪水量級的排序為100年一遇>1974年型>1957年型>50年一遇>20年一遇。以淹沒水深差值作為疊加分析的依據(jù),統(tǒng)計了0.05~0.5,0.5~1.0,1.0~1.5,1.5~2.5 m和>2.5 m等分級水深的淹沒面積,各方案的淹沒面積為76~210 km2,約為保護區(qū)面積的4.6%~12.8%,各分析方案淹沒面積變幅較大,詳見表1。淹沒面積排序為100年一遇>1974年型>1957年型>50年一遇>20年一遇,與洪量一致。通過合理性分析可知,模型
表1 沂河馬頭潰口淹沒水深和面積分級統(tǒng)計 km2Tab.1 lassification statistics of inundation depth and area of Matou breach in Yihe River
可以反映河道和保護區(qū)洪水演進和洪澇分布情況,作為區(qū)域洪澇計算平臺。
(1)基于沭西片防洪保護區(qū)地形地貌及其外河洪水特征,構(gòu)建了實時動態(tài)耦合模型水動力數(shù)學(xué)模型。河道一維模型采用典型年實測洪水資料驗證,水位和流量誤差較小;保護區(qū)二維模型通過現(xiàn)場調(diào)研查勘設(shè)置了阻水導(dǎo)水構(gòu)筑物、分布式糙率等,可以反映保護區(qū)的水流特性;模型參數(shù)設(shè)置合理,具有較高的精度,可以用于區(qū)域洪水淹沒特征分析計算。
(2)運用耦合模型,開展了沭西片防洪保護區(qū)洪澇模擬計算,潰口設(shè)置了100年一遇、50年一遇、20年一遇、1957年型、1974年型等5個洪水方案。以沂河50年一遇洪水碼頭潰口為例,潰口進洪歷時13.2 h,洪峰流量537 m3/s,進洪量833 萬m3,淹沒區(qū)總面積107.7 km2。洪水主要沿白馬河右堤向下游演進,少量洪水漫過白馬河,淹沒范圍涉及郯城鎮(zhèn)、馬頭鎮(zhèn)等11個鄉(xiāng)鎮(zhèn)。
(3)沂河馬頭潰口進洪歷時短、峰值流量大,防洪保護區(qū)內(nèi)洪水演進速度快,大部分區(qū)域淹沒水深淺,洪水將沿著白馬河向下游演進,結(jié)合淹沒水深、洪水路徑和到達(dá)時間對人員和財產(chǎn)轉(zhuǎn)移進行轉(zhuǎn)移和安置,為搶險救災(zāi)提供依據(jù)。
□