山丘區(qū)小流域洪水風(fēng)險模型研究及應(yīng)用

李曉英，談亞琦，蘇志偉，田佳樂

(1.河海大學(xué) 水利水電學(xué)院，南京 210098； 2.杭州誠禹水利科技有限公司，杭州 310008)

為預(yù)防和減輕洪水帶來的危害，抵御洪水侵襲，研究洪水演進的規(guī)律顯得十分重要。隨著計算機科學(xué)的不斷進步，越來越多的洪水演進模型得到開發(fā)研究，成為防洪減災(zāi)理論體系的組成部分和分析工具。

洪水演進數(shù)值模擬是研究洪水演進規(guī)律的有效方法之一。Caleffi等人[1]采用二維淺水方程對意大利托切河進行洪水演進模擬。Bladé E等人[2]基于有限體積法分別建立一維、二維水動力學(xué)模型，并提出基于數(shù)值通量的耦合方法，加強一維～二維模型之間的動量傳遞，在位于西班牙的埃布羅河進行驗證；李傳奇等人[3]利用構(gòu)建的一維二維水動力耦合模型，針對不同典型降雨過程，對模型進行驗證，并最終用于模擬濟南市不同重現(xiàn)期下洪水淹沒情況；苑希民等人[4]建立了一、二維水動力耦合模型，其中一維模型采用Preissmann格式離散，二維利用Roe格式離散，應(yīng)用模型于黃河青銅峽河西灌區(qū)潰堤洪水的模擬，較為真實地反映了洪水在計算區(qū)域的演進過程與淹沒范圍；吳天蛟等人[5]以MIKE11為基礎(chǔ)，區(qū)間入流采用分布式模型，對三峽庫區(qū)進行洪水演進模擬；謝作濤等人[6]從求解Holly-Preissmann格式出發(fā)，結(jié)合有限差分法建立一維洪水演進模型；槐文信等[7]建立一、二維水動力學(xué)模型對渭河下游河道及洪泛區(qū)洪水進行數(shù)值仿真模擬；付成威等人[8]建立了一、二維實時動態(tài)耦合模型，并利用濕水深和干水深理論，改進傳統(tǒng)洪水演進模型，將模型在谷堆圩蓄滯洪區(qū)進行了驗證，結(jié)果良好。

基于上述研究，本文應(yīng)用圣維南方程組作為一維模型的基本方程；應(yīng)用水流運動的能量守恒原理和質(zhì)量守恒推導(dǎo)得到平面二維非恒定流數(shù)值模型的基本方程組，進而得到二維水動力學(xué)模型。通過側(cè)向連接的方式建立一維、二維耦合水動力學(xué)模型，利用該模型于太平溪山丘區(qū)小流域，對該流域作P=10%設(shè)計暴雨條件下的洪水演進模擬分析，為該區(qū)域防洪風(fēng)險管理提供數(shù)據(jù)參考。

1 原理介紹[9-12]

1.1 一維洪水演進模型原理

以圣維南方程組為一維水動力學(xué)模型采用的基本方程。方程表示如下：

式中：Q為流量，m3/s；A為斷面面積，m2；B為河道水面寬度，m；α為動量校正系數(shù)；g為重力加速度,m/s2。

1.2 二維洪水演進模型原理

根據(jù)水流運動的質(zhì)量守恒以及能量守恒原理，可以推導(dǎo)出平面二維非恒定流數(shù)值模型的基本方程組，該方程組由連續(xù)性方程以及X、Y方向的動量方程組成：

式中：h為實際水深,m；η為水面高程,m；u、v分別為X方向和Y方向的流速,m/s；S為源匯項,m2/s；ρ0為水的密度,kg/m3；f為科式力系數(shù)；τsx、τsy分別為X方向和Y方向風(fēng)切應(yīng)力；τbx、τby分別為X方向和Y方向底部切應(yīng)力；Tih為黏滯項、紊動項和擴散項的總和。

1.3 一維、二維模型耦合方式及模型求解過程

本文以側(cè)向連解的方式完成一維洪水演進模型與二維洪水演進模型的耦合，耦合模型求解過程見圖1。

圖1 耦合模型求解過程Fig.1 Coupling model solving process

2 模擬實例[13-14]

2.1 區(qū)域概況

太平溪流域?qū)儆谏酱ㄠl(xiāng)在安吉縣的南部，該流域是苕溪流域的源頭支流，安吉縣境內(nèi)區(qū)域受饅頭山水庫控制，控制范圍內(nèi)的集雨面積約為27.6 km2，河流長約為11.9 km。該流域所屬區(qū)域包括安吉縣山川鄉(xiāng)大里村、九畝村、船村3個行政村，該區(qū)域在汛期來臨時，天然來水量十分巨大且流速迅急，給該流域及其下游防護區(qū)安全造成極大的威脅。

2.2 太平溪流域流水網(wǎng)生成

在特定的流域，徑流匯聚的多少與水流累積值成正比。因此，可用水流累積個數(shù)來定義河流的起始點。在河網(wǎng)生成的第一步中，必須先設(shè)置水流累積矩陣累積值的閾值，認為在水流累積個數(shù)小于該值的網(wǎng)格上不能產(chǎn)生足夠的徑流而形成水道，而大于該值的網(wǎng)格上能形成水道。由此，水流累積個數(shù)大于該值的網(wǎng)格所組成的流水網(wǎng)，即為一個柵格形式的河網(wǎng)，可稱之為模擬河網(wǎng)。

當采用閾值來確定河網(wǎng)起始點時，由于流域地形、下墊面產(chǎn)流機制以及流域所處區(qū)域氣候條件等因素的不同，其閾值也不同。本文基于DEM數(shù)據(jù)，生成了太平溪流域的模擬河網(wǎng)；采用自然流域分塊的方法生成流域子單元，將太平溪流域分為6個子單元(圖2)，其中子單元3、4、5位于太平溪的干流；子單元6位于上游；子單元1和2分別為區(qū)間入流。

圖2 太平溪流域模擬河網(wǎng)和劃分子單元示意圖Fig.2 A schematic diagram of simulating river network and dividing sub - unit in Taipingxi basin

2.3 太平溪流域網(wǎng)格劃分

本文基于無結(jié)構(gòu)網(wǎng)格來生成太平溪流域坡面二維洪水演進的計算網(wǎng)格。

二維洪水演進模型的邊界可以分為開邊界與陸地邊界，開邊界就是指流量過程、水位過程等人為設(shè)定一部分水體所形成的有界計算域；而陸地邊界是指水流運動的邊界，一般使用無滑移邊界來設(shè)定，認為水流不能越過陸地邊界。陸地邊界是實際存在的邊界。由于流域邊界為山脊分水嶺，一般水流不能越過分水嶺，故本文以太平溪流域邊界作為生成網(wǎng)格的陸地邊界。此外，考慮到二維水動力模型需要與一維水動力模型進行耦合，故本文在生成網(wǎng)格時將太平溪干流區(qū)域也當作單獨的邊界進行處理。生成網(wǎng)格時，每個網(wǎng)格的最大面積限制在2 000 km2以內(nèi)。生成網(wǎng)格后將太平溪內(nèi)的高程信息插值至各計算網(wǎng)格，結(jié)果見圖3。

圖3 太平溪流域二維洪水演進網(wǎng)格示意圖Fig.3 Schematic diagram of two-dimensional flood evolution grid in Taipingxi basin

3 耦合模型構(gòu)建

3.1 邊界條件

二維洪水演進模型的邊界條件包括兩類：①水陸邊界條件：此類邊界條件為圖4中二維洪水演進模擬范圍的邊界，認為水流無法越過該邊界，在水陸邊界法向處的流量及流速均為零。②流量/水位邊界條件：此類邊界條件為太平溪干流向流域坡面漫灘的過程。本文采用側(cè)向連接的方式(圖4)將一維太平溪干流與二維太平溪流域坡面進行耦合，利用堰流計算通過側(cè)向連接的水流:

式中：q為交換水量，m3/s；W為寬度，m，一般取單元格和河道相連的邊長；C為堰流系數(shù)；k為堰指數(shù)；Hus為堰上游水位,m；Hds為堰下游水位,m；Hw為堰頂高程,m。

圖4 側(cè)向連接示意圖Fig.4 Side connection diagram

3.2 動邊界處理

由于洪水在流域坡面上運動的特殊性，計算區(qū)域中會存在隨水位漲落而變化的動邊界, 為保證模型計算的連續(xù)性, 采用干/濕處理技術(shù)來對動邊界進行處理。其中，干水深設(shè)置為0.005 m，濕水深設(shè)置為0.07 m，淹沒水深設(shè)置為0.05 m。即當計算區(qū)域水深小于0.005 m時，該計算區(qū)域記為“干”，不參加計算；當水深大于0.07 m時,該計算區(qū)域記為“濕”,重新參加計算。

3.3 初始條件及模型參數(shù)

在太平溪干流未發(fā)生漫灘時流域坡面上的積水較淺，本文假定此時區(qū)域內(nèi)是無水的，認為初始水位為流域坡面的地面高程值。

由于缺少實測資料，各類地形、區(qū)域取值因地制宜，但大致保持在一定范圍內(nèi)。本文依據(jù)相關(guān)經(jīng)驗及相關(guān)文獻各糙率值，確定為如下：一、二維的河道取值為0.03，村莊取值為0.07，水田取值為0.05，樹叢取值為0.05，道路以及空地取值為0.035。

4 太平溪流域洪水淹沒風(fēng)險分析

選取P=10%頻率下的50 h的設(shè)計暴雨，基于本文構(gòu)建的耦合洪水演進模型對太平溪流域的洪水淹沒風(fēng)險進行分析。圖5為P=10%設(shè)計暴雨條件下不同計算單元的出流過程。

圖5 太平溪流域P=10%設(shè)計洪水模擬過程 (子單元1、2、6)Fig.5 P=10% design flood simulation process in Taipingxi basin(subunit 1、2、6)

根據(jù)以上邊界條件模擬了洪水太平溪干流及坡面上的演進過程。圖6為太平溪干流在P=10%時最高水位的沿程分布情況，最高水位出現(xiàn)在第20 h，也即邊界入流達到峰值的時刻。從圖6可以看出，太平溪干流大多數(shù)斷面所面臨的風(fēng)險情況較小，水位均在干流主槽內(nèi)，未發(fā)生漫灘的情況。但仍有3個斷面發(fā)生了漫灘的風(fēng)險事件，這3個斷面位置見圖7，其計算起點距分別為1 913.4、3 044.9和3 676.7 m。

圖6 太平溪干流P=10%最高水位分布Fig.6 P=10% maximum water level distribution in Taipingxi mainstream

圖7 太平溪干流P=10%發(fā)生風(fēng)險事件斷面示意圖Fig.7 P=10% risk profile schematic in the Taipingxi basin mainstream

圖8～圖11分別為設(shè)計暴雨時間內(nèi)35 h的太平溪流域坡面淹沒的范圍及深度示意圖。從這些圖可以看出，由于漫灘事件的發(fā)生，從設(shè)計暴雨的20 h起洪水由太平溪干流進入流域坡面進行傳播，發(fā)生淹沒的區(qū)域位于風(fēng)險斷面2，并造成太平溪流域船村的風(fēng)險事件。但在P=10%設(shè)計暴雨條件下，洪水淹沒的范圍相對有限，主要集中在干流附近，并且主要向下游方向進行演進，未對船村造成過大的影響。從淹沒發(fā)生的時間過程來看，發(fā)生漫灘事件后2 h(圖9)的淹沒深度達到最大值，隨后隨著洪水的擴散與歸槽，淹沒深度逐漸減小，發(fā)生漫灘事件后10 h后的淹沒高程已顯著降低。

圖8 太平溪流域淹沒范圍及深度示意圖(19 h)Fig.8 Flooding scope and depth schematic in Taipingxi basin(19 h)

圖9 太平溪流域淹沒范圍及深度示意圖(22 h)Fig.9 Flooding scope and depth schematic in Taipingxi basin(22 h)

圖10 太平溪流域淹沒范圍及深度示意圖(25 h)Fig.10 Flooding scope and depth schematic in Taipingxi basin(25 h)

圖11 太平溪流域淹沒范圍及深度示意圖(35 h)Fig.11 Flooding scope and depth schematic in Taipingxi basin(35 h)

通過對比設(shè)計暴雨過程、徑流過程以及淹沒過程發(fā)現(xiàn)，由于太平溪流域的匯流時間較短，故暴雨中心到洪峰的時間以及漫灘風(fēng)險發(fā)生的時間均較短，在P=10%設(shè)計暴雨條件下的時間僅為1 h。為此，需要建立洪水預(yù)報警報系統(tǒng)。通過將實時降雨、水文等數(shù)據(jù)輸入模型得到預(yù)報信息，在必要的情況下，啟用警報系統(tǒng)，為抗洪救災(zāi)、居民撤離等提供寶貴的準備時間，減少經(jīng)濟損失。

從工程措施來看，通過提升3個風(fēng)險斷面的堤防可以防止漫灘事件的發(fā)生。具體來說，風(fēng)險斷面1加高5 cm、風(fēng)險斷面2加高15 cm、風(fēng)險斷面3加高20 cm。通過數(shù)值實驗分析發(fā)現(xiàn)，在加高以上3個斷面的高程后，太平溪流域在P=10%設(shè)計暴雨條件下無風(fēng)險事件發(fā)生。

5 結(jié) 語

1) 基于DEM數(shù)據(jù)生成太平溪流域的模擬河網(wǎng)，基于無結(jié)構(gòu)網(wǎng)格生成流域坡面二維洪水演進的計算網(wǎng)格，以側(cè)向連解方式完成一維洪水演進模型與二維洪水演進模型的耦合。

2) 利用建立的耦合模型對太平溪流域P=10%設(shè)計暴雨條件下的洪水淹沒進行分析，并根據(jù)分析結(jié)果，為山丘區(qū)太平溪流域提出風(fēng)險應(yīng)對對策及工程和非工程措施。

3) 采用基于數(shù)學(xué)模型技術(shù)的洪水過程研究，可以準確地分析流域的洪水過程特征和洪水風(fēng)險因子，利用耦合模型進行洪水風(fēng)險分析研究將成為洪水風(fēng)險管理的一項重要的技術(shù)支撐。