江春波,周 琦,申言霞,柳高飛,張 帝
(1.清華大學 水沙科學與水利水電工程國家重點實驗室,北京 100084;2.大連理工大學 建設工程學部水利工程學院,遼寧 大連 116024)
洪澇災害受自然環(huán)境變化和人類高強度活動雙重影響,是全球影響范圍較廣、對人類生存和發(fā)展危害較為顯著的自然災害之一。隨著全球氣候變化引發(fā)的極端暴雨事件增加,洪澇災害呈現(xiàn)出突發(fā)性強、峰高量大和危害面積廣等特征。近年來我國諸多流域尤其是中小流域洪澇災害頻發(fā),給當?shù)厝嗣裆踩拓敭a(chǎn)帶來嚴重的損失。根據(jù)文獻[1],全國山洪災害防治區(qū)面積386 萬km2,涉及受威脅村莊57 萬個,人口3 億人,其中直接受威脅人口近7000 萬人。2010年以來,水利部在全國29個省(自治區(qū)、直轄市)和新疆生產(chǎn)建設兵團的305 個地(市)、2076 個縣(市、區(qū)、旗、團場等)持續(xù)開展山洪災害防治工作,共投入項目建設資金352 億元;非工程措施建設資金269 億元,重點山洪溝防洪治理資金83 億元。我國洪澇災害防治取得了顯著的效果,2011—2019年因山洪災害平均死亡353 人,較項目實施前的2000—2010年平均死亡1179 人大幅減少7 成。其中2018年死亡人口129人,為歷史最低。而在剛剛過去的2020年,我國再次遭受嚴重洪澇災害。據(jù)應急管理部統(tǒng)計[2],2020年7月份全國洪澇災害造成3817.3 萬人次受災,56 人死亡失蹤,299.6 萬人次緊急轉(zhuǎn)移安置;2.7萬間房屋倒塌,24 萬間不同程度損壞;農(nóng)作物受災面積3868.7 千公頃;直接經(jīng)濟損失1097.4 億元。與近5年同期均值相比,7月份洪澇災害因災死亡失蹤人數(shù)、倒塌房屋數(shù)量分別下降74.2%和67.3%,但受災人次、緊急轉(zhuǎn)移安置人次分別上升62.5%、88.6%。因此,盡管我國洪澇風險管理取得了顯著成績,但氣候變化和城市化背景下洪澇風險管理的問題與挑戰(zhàn)依然嚴峻,需要深刻梳理與針對性適應。
因山區(qū)流域洪澇災害突發(fā)性強、預見期短,降雨及洪澇時空變化的復雜性和不確定性日益突出,由此引發(fā)的洪澇災害日益嚴峻,洪澇災害治理任務更加艱巨。及時準確地評估洪澇災害及經(jīng)濟損失,是提高流域洪澇災害監(jiān)測、預報與風險防范能力建設的主要內(nèi)容,為減少洪澇災害損失提供技術(shù)支撐。
目前對于洪澇災害預報,常采用水文模型或水動力模型,兩種模型的研究方法和計算過程側(cè)重點不同。單一的水文或水動力模型,因模型功能的局限性,其洪澇預報結(jié)果可靠性難以滿足流域防洪減災決策的需要。為了提高流域洪澇災害模擬和預報結(jié)果的可靠性,研究流域水文過程與水動力過程的互饋機制,開發(fā)能夠更加合理地反映流域洪澇形成及發(fā)展過程的模型,評價工程措施對減緩洪澇災害損失的效果并優(yōu)化工程設計,是當前洪澇災害治理需要解決的問題。
流域洪澇模擬和預報是基于已知的氣象、水文、地形地貌、水工建筑等條件,利用數(shù)值模型和監(jiān)測設備預報洪澇發(fā)生、發(fā)展的科學方法,其發(fā)展過程也是由多個學科共同決定的[1]。洪澇模擬及預報涉及數(shù)據(jù)的采集、傳輸存儲、后臺計算、前端顯示等多個環(huán)節(jié),每個環(huán)節(jié)都必須有相應的技術(shù)支持才能實現(xiàn)。洪澇預報從早期的經(jīng)驗公式模式逐漸演變成當前囊括遙感和地理信息技術(shù)、計算機科學、水文學、水動力學、數(shù)據(jù)庫管理等多學科、多層次的綜合性系統(tǒng)學科,充分利用水文、水力信息,準確、快速地模擬真實的洪澇演進情景,為洪澇災害治理工作提供科學支持。
洪澇預報數(shù)學模型是洪澇災害治理信息決策支持系統(tǒng)的核心,模型的精度直接影響洪澇計算結(jié)果的可靠性及防洪減災工程布置的成效。研制洪澇預報模型是一項理論性強、應用要求較高的工作[4]。流域洪澇模擬和預報模型可以分為水文模型、水動力模型和水文與水動力耦合模型,下面分別對這三類模型進行介紹。
2.1 水文模型水文模型的發(fā)展過程是伴隨著對自然界水循環(huán)過程的不斷深入認識而逐漸改進和完善的,其目的在于能夠準確模擬大氣、地表、土壤中水的遷移轉(zhuǎn)化過程[4]。水文模型包括產(chǎn)流和匯流過程,根據(jù)降雨、蒸發(fā)、入滲等信息進行產(chǎn)流水量計算,通過坡面流理論或經(jīng)驗公式法進行匯流計算,最終輸出流域或子流域出口處的流量隨時間變化過程。構(gòu)建水文模型需要考慮水文循環(huán)的各個環(huán)節(jié),對表征各環(huán)節(jié)的數(shù)學方程進行合理選擇,計算所需參數(shù)種類較多,包括氣象、土壤條件和土地利用、植被覆蓋等多個方面,需要大量數(shù)據(jù)進行率定驗證后才能進行實際應用[5-6]。早期的簡單水文模型包括經(jīng)驗公式、推理法、馬斯京根法[6-7]。20 世紀中期以后,逐漸形成了“集總式水文模型”和“分布式水文模型”兩種模型并沿用至今,當前應用的大部分模型都可歸于這兩類[8],還有部分模型介于兩者之間被視為半分布式模型,如TOPMODEL[9]。
集總式水文模型不考慮流域內(nèi)部地質(zhì)、地貌、土壤、植被等要素空間分布不均勻性對水文循環(huán)的影響,將流域作為一個整體進行研究,如水箱模型[10]和新安江模型[11]。該類水文模型忽略了流域內(nèi)降雨分布、下墊面空間變化等因素的差異性。隨著社會經(jīng)濟發(fā)展,為了反映流域下墊面變化對徑流的影響,流域水土流失治理等對環(huán)境的影響,研發(fā)了分布式水文模型[12]。分布式水文模型根據(jù)現(xiàn)實流域中降雨、土壤、下墊面等地理信息的空間分布不均勻特征,將流域劃分為若干子流域,在每個子流域上計算不同水文要素變化對水循環(huán)過程的影響。分布式水文模型可分為松散型和耦合型兩類[13]。相對于耦合型模型,松散型水文模型假定每個響應單元對整個流域響應的貢獻互不干擾,通過每個單元的疊加確定整個流域響應。該求解算法簡單,但反映徑流形成機制不夠完善。耦合型水文模型,考慮各個水文子單元的產(chǎn)匯流相互影響,其精度優(yōu)于其它類水文模型。陳海梅等[14]采用擴散波水文模型求解坡面流,分析了擴散波模型的精度和數(shù)值格式穩(wěn)定性,結(jié)果表明擴散波水文模型在坡面流計算中具有較好的精度,格式穩(wěn)定性和計算效率高于求解淺水方程。將水文模型與城市排水管網(wǎng)一維水力流計算耦合,可以將水文模型延伸到城市洪澇治理中。如文獻[15-16]采用了地表-路網(wǎng)-管網(wǎng)-河網(wǎng)的四層結(jié)構(gòu)。文獻[15]的地表徑流模型為二維擴散波方程,缺少淹沒區(qū)的二維水動力模型,水文模型與城市街道流動動態(tài)單向耦合,街道水流與城市排水管網(wǎng)一維流動動態(tài)雙向耦合,管網(wǎng)水流最后流向河道。文獻[17]構(gòu)建了匯流過程精細化城市雨洪模型,采用地表徑流-路網(wǎng)-管網(wǎng)-河網(wǎng)的四層結(jié)構(gòu),通過擴散波水文模型計算地表二維匯流,地表水匯入路網(wǎng)后再通過雨水口匯入管網(wǎng),最終流入排水河道。
綜上,水文模型應用廣泛,計算效率高。但計算結(jié)果是根據(jù)流域水文條件獲得產(chǎn)匯流的流量過程,不能提供淹沒區(qū)洪水演進的動力特征。
2.2 水動力模型水動力模型不僅能模擬水流在河道及河漫灘的演進過程,也可模擬泛濫洪水在受堤防保護的城市與農(nóng)村的演進過程,以及由暴雨形成的洪水內(nèi)澇過程,輸出結(jié)果通常是積水淹沒范圍的時空變化過程及淹沒區(qū)域特定位置的水深和流速。水動力模型的控制方程包括質(zhì)量守恒方程和動量守恒方程。模型涉及的參數(shù)相對較少,如河床糙率系數(shù)等。
水動力模型根據(jù)是否考慮水力要素的橫向和垂向變化,可分一維、二維和三維模型[18]。一維模型只考慮流動要素在河道順流方向的變化,建模所需數(shù)據(jù)量少且計算效率高,但是計算精度不高。常見的一維模型包括求解非恒定流的圣維南方程,求解恒定流的伯努利方程;二維模型,如二維淺水方程(SWE),常用于考慮沿河道橫向水力要素變化的河湖及低洼積水區(qū),適用對江、湖、河口等區(qū)域的水位和流速分布的描述;三維模型可考慮水力要素沿垂向的變化,常用于江河入???、城市大型地下蓄水隧洞進口附近的復雜流態(tài)條件下的水動力特性研究。在流域洪澇模擬和預報中,通常一維和二維水動力模型即可滿足要求,應用三維模型的情況不是很多。對于城市地下大型蓄水隧洞入口流態(tài)、排水管內(nèi)網(wǎng)明滿流等可以采用三維模型,因計算相對復雜費時,在流域洪澇模型中不考慮。
水動力模型的數(shù)值求解方法包括有限差分法、有限元法和有限體積法等[13]。差分法數(shù)值離散格式簡單易懂,通過前差分、后差分和中心差分等算法,可以構(gòu)造出不同精度的差分格式,將偏微分方程轉(zhuǎn)換為差分方程求解。有限差分法多使用結(jié)構(gòu)化網(wǎng)格,不便于網(wǎng)格局部加密,同時矩形構(gòu)造網(wǎng)格劃分對復雜邊界形狀的適用性較低。使用差分法的模型有英國Wallingford 和Halcrow 公司開發(fā)的ISIS 模型[19]、丹麥水力學研究所開發(fā)的MIKE 11 模型[20-21]、美國陸軍工程兵團水文工程中心開發(fā)的HEC-RAS 模型[22]、澳大利亞WBM 公司開發(fā)的TUFLOW 模型[23]等。國內(nèi)采用有限差分法對洪澇模擬的研究也比較成熟,如文獻[24]采用有限差分法求解二維淺水方程,在二維差分網(wǎng)格內(nèi)部“開溝”以模擬山區(qū)小河流的行洪過程、提高模型的空間分辨率;并引入水文學法考慮不同類型土地的入滲飽和過程,將流域降雨轉(zhuǎn)為形成地表徑流的凈雨量。由水力學與山區(qū)河流開發(fā)保護國家重點實驗室(四川大學)研發(fā)的WWL(West Water-Lateral),是立面二維水溫分析模型[25],采用有限差分格式對立面二維方程進行離散。文獻[26]采用有限差分法求解二維水動力方程,模擬近海風暴潮。文獻[27-28]將一維水動力與二維水動力模型進行耦合,模型不但保證耦合節(jié)點的質(zhì)量守恒,還保證了耦合節(jié)點的動量守恒。
有限元法求解問題的基本步驟是將所討論問題的域劃分成若干微小單元,選取基函數(shù)將節(jié)點離散變量轉(zhuǎn)化成連續(xù)變量,先在單元上積分形成單元系數(shù)矩陣,再合成全區(qū)域的整體方程,得到關(guān)于未知量的代數(shù)方程組。常規(guī)的有限元格式是針對結(jié)構(gòu)力學問題的,不適合求解洪水運動的雙曲方程組。對于洪水運動的模擬需要選取適合對流比較強的高精度有限元格式,數(shù)值格式穩(wěn)定且計算精度良好的有限格式也有很多,這里不準備展開介紹。
有限體積法是當前水動力模型中比較常用的數(shù)值計算方法,F(xiàn)luent[29]、MIKE 21[30]、OpenFOAM[31]等均采用有限體積法。Godunov 型格式將水流運動的數(shù)值計算近似為局部黎曼問題,數(shù)值格式精度良好[32]?;诤闈衬M結(jié)果,結(jié)合社會經(jīng)濟信息可進一步開展洪澇災害風險評價,如Mike Flood 模型的應用[33]。張大偉等[34]建立一維、二維潰壩水流耦合數(shù)學模型,采用非構(gòu)造網(wǎng)格對控制方程進行空間離散。文獻[35]采用有限差分法和有限體積法相結(jié)合的數(shù)值格式,對一維水動力方程和二維水動力方程進行求解,并將水文模型與水動力模型進行耦合,建立FRAS(Flood Risk Analysis System)模型。
大連理工大學的Hydroinfo[25]是基于水動力模型的軟件,在連續(xù)方程中考慮降雨、入滲等產(chǎn)流水量,匯流及洪水運動通過水動力模型計算。該模型具有計算簡便、功能齊全等優(yōu)點,包括了各種水文模型、一維水動力與二維水動力的動態(tài)耦合、城市管網(wǎng)模型和泥沙、污染物輸移模型,具有廣泛的應用價值。
近年來,隨著計算機技術(shù)的快速發(fā)展,一些涉及水動力方面的商業(yè)模型,如Mike 系列、EFDC、HEC-RAS 等模型,在工程計算中經(jīng)常應用。為了開發(fā)精確良好/數(shù)值計算穩(wěn)定的洪澇模擬和預報模型,亞利桑那大學于2014年和2017年提出CHRE2D[36-37],該模型屬于流域水動力模型,需在新的時間步計算之前判斷控制體單元內(nèi)實際水深與水深閾值的關(guān)系,被動切換淺水模型和擴散波模型的使用區(qū)域。在水深較淺區(qū)域,因數(shù)值格式穩(wěn)定性差,用運動波或擴散波方程代替淺水方程的求解。該模型采用矩形均勻網(wǎng)格,時空步長均受淺水方程的穩(wěn)定條件控制,在應用于中小流域規(guī)模的洪澇預報時,其計算效率有待提高。文獻[38]提出的洪澇預報模型,采用非構(gòu)造網(wǎng)格在流域上求解淺水方程,通過薄層水假定提高水深較淺區(qū)域的數(shù)值格式的穩(wěn)定性,模型適合城市的洪澇預報。西安理工大學提出的GAST 模型(GPU Accelerated Surface Water Flow and Transport Model)[39],在二維淺水方程中考慮水文產(chǎn)匯流計算功能,使用二階Godunov 格式有限體積法來保證計算精度,并采用GPU 加速技術(shù)提高運算效率。
綜上所述,單純的水動力模型不能滿足流域洪澇模擬與預測的需求,將水文計算公式與全動力模型進行耦合,水動力模型在流域洪澇模擬與預報中發(fā)揮了重要應用。二維水動力模型在連續(xù)方程中考慮凈雨量,降雨、入滲和蒸發(fā)等產(chǎn)流水量,通過水文模型計算凈雨量。對于中小流域及城鎮(zhèn)的洪澇問題,全動力模型具有計算分辨率高等優(yōu)點[24,34-39]。為了提高計算效率和數(shù)值格式穩(wěn)定,可在二維水動力計算中使用粗網(wǎng)格[24,35],將一維河道和管網(wǎng)、行洪街道作為特殊通道,通過特殊通道實現(xiàn)了一維水動力與二維水動力模型的耦合。粗網(wǎng)格上對淺水方程進行空間離散,時間離散采用顯式格式,模型計算節(jié)省時間;全流域使用二維淺水方程,在水深很薄區(qū)域切換到二維運動波或擴散波方程[36-37],保證數(shù)值格式穩(wěn)定;在全流域求解二維淺水方程時,判斷水深較淺的網(wǎng)格,使用薄層水假設[38,62]、保證數(shù)值格式的穩(wěn)定性;根據(jù)二維淺水方程求解產(chǎn)匯流及洪水運動,通過開發(fā)GPU 加速算法[39],保證計算精度與計算效率的配置。
節(jié)省計算時間和保證數(shù)值格式穩(wěn)定的另外一個途徑,是將流域在空間上分成淹沒區(qū)和非淹沒區(qū),在包括河湖、漫灘和低洼積水區(qū)的淹沒區(qū)域采用二維水動力計算[40-45],淹沒區(qū)采用比較精細的計算網(wǎng)格,網(wǎng)格能反映堤防寬度、河道寬度等局部尺寸,細網(wǎng)格為城市排水管網(wǎng)入口位置等處水深確定提供更加合理的水力條件;在淹沒區(qū)域外側(cè)(水文產(chǎn)匯流區(qū)域)采用水文模型,水文計算獲得的產(chǎn)匯流的流量過程,作為水動力計算的邊界條件。水文模型與水動力分區(qū)計算可以達到數(shù)值格式穩(wěn)定和節(jié)省計算時間的目的,適合面積較大的山區(qū)流域洪澇模擬與預報。探討水文與水動力模型的不同耦合方式,可以在提高計算精度的同時提高計算效率與格式穩(wěn)定性。流域分區(qū)計算,水文模型與二維水動力模型動態(tài)耦合,可以根據(jù)需要切換到流域全動力計算,這種分區(qū)耦合模型將具有流域全動力模型的功能。因此水文與水動力動態(tài)耦合模型是具有發(fā)展?jié)摿Φ哪P椭弧?/p>
2.3 水文與水動力模型耦合將水文與水動力模型進行耦合,利用水文模型和水動力模型的各自優(yōu)點,可以分別彌補水文模型和水動力模型的不足,是流域及城市洪澇模擬和預報的發(fā)展趨勢。根據(jù)水文過程與水動力過程的連接關(guān)系及計算時間順序,可以歸納為水文與水動力的串聯(lián)耦合模型、水文與水動力的動態(tài)單向耦合模型、水文與水動力的動態(tài)雙向耦合模型,如圖1 所示。圖中符號T 為總計算時間長度,t 為當前計算時刻。
如圖1(a),串聯(lián)耦合模型從開始時刻到計終了時刻先獨立計算水文過程,然后獨立計算水動力過程,水文計算獲得的在子流域出口的流量作為水動力計算的邊界條件,水文過程影響水動力過程,水動力過程對水文過程沒有影響。如圖1(b),水文與水動力的動態(tài)單向耦合模型,是在時刻t 同時計算水文模型和水動力模型,水文計算影響水動力計算,但水動力計算不影響水文計算。水文與一維水動力動態(tài)單向耦合模型在某種條件下與串聯(lián)耦合模型等效的。如將動態(tài)單向耦合節(jié)點的流量信息儲存記憶,在水文計算結(jié)束后,將耦合點的流量加載到一維水動力計算模型中,動態(tài)單向耦合就相當于串聯(lián)耦合模型。如圖1(c),動態(tài)雙向耦合模型,在時刻t 水文與水動力同時計算,水文計算信息影響水動力計算,水動力計算信息影響水文計算,水文與水動力的影響是雙向的。高分辨率的分布式水文模型與二維淺水方程的動態(tài)雙向耦合,顯然水文模型和水動力模型均具有高分辨率,這樣的動態(tài)雙向耦合模型國內(nèi)外還沒有文獻報道,也是將要建立的耦合模型。
圖1 水文模型與水動力模型耦合
將水文和水動力模型進行更合理的耦合,充分反映水文過程與水動力過程的互饋機制,提高洪澇預報結(jié)果的可靠性,還需要進行深入探討,宋利祥等[46]指出水文與水動力耦合模型的本質(zhì)區(qū)別在于是否合理考慮了水文與水動力過程的相互影響,尤其是水文驅(qū)動水動力模型的邊界形式以及局部淹沒狀態(tài)對匯水單元水文響應過程的影響。因此亟需在探索水文水動力自適應耦合機制的基礎上,建立水文與水動力動態(tài)雙向耦合的洪澇模擬和預報模型,可為流域洪澇水文水動力過程精細化與高效模擬奠定理論基礎。
2.3.1 水文與水動力模型串聯(lián)耦合 水文與水動力耦合模型中應用最成熟的耦合方法是“串聯(lián)耦合”。從洪澇過程的起始時刻到結(jié)束時刻,串聯(lián)耦合模型首先獨立運行水文模型,進行產(chǎn)匯流計算,獲得河道上游斷面及河道支流控制斷面處的流量隨時間變化過程。將這個流量過程作為水動力模型的輸入邊界條件,然后獨立進行水動力計算,獲得河道、河漫灘等積水區(qū)的水流運動信息。水文計算和水動力計算兩個過程在空間上是相互獨立的,屬于前后的串聯(lián)關(guān)系。
由于各種水文模型均可輸出控制斷面的流量過程,容易滿足串聯(lián)耦合模擬的需求,常用于水文與水動力串聯(lián)耦合。McMillan 利用降雨產(chǎn)流模型與一維動力波方程,將降雨產(chǎn)流過程作為動力波方程的上游邊界,分析不同頻率降雨產(chǎn)生的洪水對城市的淹沒情況[47]。采用串聯(lián)耦合模型的有Mon?tanari[48]、Choi[49]和Grimaldi[50]等,研究將 水文模型 的產(chǎn)流過 程 作為一 維、二維水 動力模型 的上游流 量邊界,驅(qū)動水動力模型運行。一些串聯(lián)耦合模型的應用與研究,還可以參考文獻[51-52]。此外,一些成熟的模型如SWAT[53]、HEC-HMS[54]、MIKE系列[55-56]等的應用,多是采取水文模型與水動力模型的串聯(lián)耦合方式。采用基于地貌特征的流域分布式水文模型GBHM與MIKE11一維水動力模型進行串聯(lián)耦合,模擬三峽區(qū)間降雨產(chǎn)匯流過程。其它水文與水動力耦合模型研究與應用,參見參考文獻[57-60]。
串聯(lián)式耦合模型應用方便,對于較小的空間區(qū)域可以獲得合理的精度。但是對于較大范圍的流域洪澇預報問題,水文過程與水動力過程的互饋作用更加敏感。水文與水動力串聯(lián)耦合模型,僅在包括河湖、河漫灘和低洼積水區(qū)進行水動力計算,在積水區(qū)之外的流域進行水文產(chǎn)匯流積水,水文計算結(jié)果為水動力計算提供邊界條件。在應用水文與水動力的串聯(lián)耦合模型時,需要確定淹沒區(qū)周圍流量邊界點的位置和邊界點數(shù)量的合理性,保證每個邊界點的流量隨時間變化與實際匯流過程相符,盡量避免不合理的邊界點位置和流量過程的選取,減小模擬和預報結(jié)果的誤差。事實上,在暴雨過程中洪澇淹沒區(qū)范圍隨時間變化,水動力模型的流量邊界點位置也隨時間變動,水文匯流進入淹沒區(qū)的邊界點在空間上是連續(xù)分布的,不是有限個分散點,而串聯(lián)耦合模型在確定這些流量邊界點位置具有一定難度。
因此需要探討水文與水動力模型的動態(tài)耦合,在時間上同步求解水文模型和水動力模型,兩種模型需在空間上建立合理的信息交換方式,預期可提高流域洪澇模擬及預報結(jié)果的精度。水文與水動力的動態(tài)耦合模型,水文產(chǎn)匯流計算和淹沒區(qū)水動力計算在時間上同時進行。根據(jù)水文與水動力的相互作用關(guān)系,動態(tài)耦合模型又分為“單向”耦合和“雙向”耦合模型,參考文獻中除文獻[13]和[61]中研制的模型外,已有的動態(tài)耦合模型均為動態(tài)單向耦合模型,也常被稱為動態(tài)耦合模型。對于單向耦合模型,水文產(chǎn)匯流計算結(jié)果是水動力計算的驅(qū)動條件,但水動力計算結(jié)果不影響水文計算;雙向耦合模型中水文計算結(jié)果作為水動力計算的驅(qū)動條件,而水動力計算結(jié)果(如水位壅高等)也會對水文匯流計算產(chǎn)生影響。
2.3.2 水文與水動力模型動態(tài)單向耦合 水文產(chǎn)匯流模型與一維水動力模型在時間上同步進行計算,已有的模型可分為水文與一維水動力模型的動態(tài)單向耦合模型和水文與二維/一維水動力的動態(tài)單向耦合模型。
(1)水文與一維水動力模型動態(tài)單向耦合。MIKE 系列軟件包含了水文模型、水動力模型等多個模塊,具有水文與水動力模型動態(tài)耦合功能,如MIKE-SHE 和MIKE11 的耦合屬于動態(tài)單向耦合[56]。水文計算采用分布式水文模型(SHE),水動力計算采用一維圣維南方程。兩個模型在連接通道根據(jù)水位關(guān)系計算兩種模型之間的流量交換,在指定節(jié)點處發(fā)生水量交換(參見圖2)。因一維水動力模型在模擬堤外洪澇積水區(qū)具有局限性,如一維水動力模型將河道外的匯流水量強制性地加載到河道里,但實際上降雨匯流水體是從河道兩側(cè)向河道方向匯集的,因此水文與一維水動力的動態(tài)耦合模型往往不能準確反映堤防外側(cè)低洼區(qū)的洪澇積水。將MIKE11 與MIKE21 耦合,組成MIKE Flood,可以克服一維水動力模型的應用限制。
圖2 MIKE SHE/MIKE 11 動態(tài)耦合示意
(2)水文模型與二維水動力模型的間接動態(tài)耦合。中國水利水電科學研究院提出CUHHM(Cou?pled Urban Hydrological-Hydrodynamic Model)[40-41]是水文模型與一維水動力模型進行動態(tài)單向耦合。對流域進行分區(qū)計算,水文產(chǎn)匯流區(qū)域采用分布式水文模型計算,河道和排水管網(wǎng)進行一維水動力計算,水文計算與一維水動力計算直接動態(tài)耦合。為了考慮河湖、漫灘及低洼積水區(qū)的二維水動力特征,在局部區(qū)域進行二維水動力計算,二維水動力計算與一維水動力計算動態(tài)雙向耦合,屬于水文產(chǎn)匯流模型與二維水動力模型的間接耦合。二維水動力計算可為排水管網(wǎng)入口提供更準確的水位條件,地表二維水流排入管網(wǎng),二維計算影響一維計算,管網(wǎng)水流外溢也影響二維水動力計算。一維水動力與二維水動力的耦合,類似于MIKE Flood,但有許多改進。文獻[41]采用TELEMAC 有限元格式在局部區(qū)域求解二維淺水方程,為保證有限元格式的穩(wěn)定性,推薦時間隱式格式和求解非線性代數(shù)方程組方法,每一個時間步內(nèi)求解聯(lián)立方程組或非線性方程組。對于二維計算區(qū)域面積有限的情況下(不是全部流域求解二維淺水方程),時間隱格式計算花費時間不會過多,計算精度良好。CUHHM 模型功能齊全,不但包括了水文與一維水動力的直接動態(tài)單向耦合,還包括了二維水動力模型與地下管網(wǎng)的動態(tài)雙向耦合計算,模型已經(jīng)應用于北京市等地的洪澇模擬和預報。
華南理工大學提出的IHUM(Integrated Hydrology and Hydrodynamics Urban Flood Model)[42-43]是SWMM 分布式水文模型與一維管網(wǎng)、一維河道水動力模型的動態(tài)單向耦合,這部分為MIKE Urban 類似,但有很多改進?;谧訁R水單元的水文計算結(jié)果為包括管網(wǎng)入口的一維水動力計算提供水位、流量條件。二維水動力計算僅在低洼積水等局部區(qū)域進行,也是水文產(chǎn)匯流模型與二維水動力模型的間接耦點,一維水動力計算與局部區(qū)域的二維水動力計算進行動態(tài)雙向耦合,這類似于MIKE Ur?ban, 但有不少改進之處,如通過正向、側(cè)向等連接方式實現(xiàn)了一維水動力模型與二維水動力模型動態(tài)雙向耦合,保證水量守恒,局部區(qū)域的二維水動力計算結(jié)果可為一維管網(wǎng)入口提供更加合理的水位、流速信息。IHUM 采用時間隱格式在局部區(qū)域求解二維淺水方程,隱格式具有較好的數(shù)值格式穩(wěn)定性,但時間隱格式需要聯(lián)立求解代數(shù)方程組。二維水動力計算僅在局部區(qū)域計算,不會花費過多的計算時間,并能獲得良好的計算精度。
美國密西根大學2012年提出的TRIBS- OFM 模型是水文產(chǎn)流模型與二維水動力模型動態(tài)單向耦合[62],由于水文模型僅計算產(chǎn)流(不計算匯流),坡面匯流及河道洪水演進均采用二維淺水模型進行計算,導致全流域運行二維水動力模型,對于大范圍的流域洪澇預報,模型的運行時間相對較長。這樣的全動力模型精度較高,但計算效率及格式穩(wěn)定性均有待提高。
中國水利水電科學研究院建立的FRAS 模型[35,63]是水文產(chǎn)流模型與二維水動力模型的動態(tài)單向耦合,水動力計算結(jié)果不影響水文計算。通過特殊通道法實現(xiàn)包括排水管網(wǎng)的一維水動力與二維水動力的動態(tài)雙向耦合,其不同維數(shù)的水動力模型耦合方式具有獨創(chuàng)性。該模型中水文計算采用SCS 模型,將水文產(chǎn)流水量作為源項加載到水動力模型中,通過水動力模型計算坡面匯流及河道中洪水演進過程。該模型功能齊全,水動力計算包括一維河道及管網(wǎng)、二維淺水模型等,在城市洪澇模擬和預報方面,模型具有良好的精度。FRAS 模型中的水文模型僅計算產(chǎn)流,匯流過程和河道淹沒區(qū)的水動力計算均需采用二維淺水方程,模型已經(jīng)成功地應用于城市洪澇模擬和預報。對于二維水動力計算區(qū)域面積有限的情況下,通過二維稀疏網(wǎng)格與一維精細網(wǎng)格的耦合,F(xiàn)RAS 模型的計算效率良好。FRAS[35,63]模型提出時間要早于TRIBS- OFM[62]模型,排水管網(wǎng)的計算功能優(yōu)于TRIBS- OFM 模型。
2.3.3 水文產(chǎn)匯流與二維水動力模型直接動態(tài)雙向耦合 如圖3 所示,淹沒區(qū)的空間范圍隨著降雨條件和時間在不斷變化,水文產(chǎn)匯流的流量匯入淹沒區(qū)的位置也是隨時間變化的,水文與水動力過程的影響是雙向的,即水文徑流影響水動力過程,水文產(chǎn)匯流流量是水動力計算的邊界條件;反過來淹沒區(qū)的水動力條件變化,如淹沒區(qū)水位變化會影響水文匯流計算。將陸地水文產(chǎn)匯流區(qū)稱為非淹沒區(qū),河道等有水區(qū)域稱為淹沒區(qū),連接兩個區(qū)域的邊界位置是動態(tài)變化的,在動邊界上不但有水量的傳遞,還有動量的傳遞。清華大學提出的水文與二維水動力的動態(tài)雙向耦合模型(Dynamic Bilat?eral Coupling Model,DBCM)[13,61],尚處于研制階段,模型能夠更真實反映流域水文過程與水動力過程的互饋機制,除文獻[61]的模擬結(jié)果外,目前國內(nèi)尚無水文產(chǎn)匯流模型與二維水動力模型直接動態(tài)耦合的研究成果報道。
圖3 流量邊界點位置變化示意
圖4 為動態(tài)雙向耦合模型框圖,將計算流域在空間上劃分成洪澇淹沒區(qū)和非淹沒區(qū),隨著洪水的漲落,淹沒區(qū)和非淹沒區(qū)的范圍隨時間變化,連接這兩個區(qū)域的動邊界位置也隨著時間變化。將兩個區(qū)域動邊界位置變化范圍作為過渡區(qū)。非淹沒區(qū)的水文計算采用擴散波、運動波或其它精度良好的分布式水文模型;淹沒區(qū)采用一維或二維動力波,僅對河道、河漫灘等低洼積水等局部區(qū)域采用二維淺水方程求解。在不同區(qū)域采用不同的數(shù)學模型,在保證模擬精度的同時,預期能提高計算效率。淹沒區(qū)通常需要反映堤防橫向斷面等微觀尺寸,網(wǎng)格劃分較細,網(wǎng)格的空間尺寸較小,受庫朗條件限制,時間步長也較小。而非淹沒區(qū)的空間范圍遠遠大于河道等淹沒區(qū)的面積,水文單元網(wǎng)格尺寸可以遠遠大于水動力計算的網(wǎng)格尺寸,時間步長也可相對較大。水文區(qū)域的較大網(wǎng)格尺寸和時間步長,與非淹沒區(qū)采用二維水動力模型相比可以明顯節(jié)省模型的計算時間。二維擴散波水文計算和二維水動力計算均采用時間顯式格式,避免聯(lián)立求解代數(shù)方程組,是提高計算效率的一個因素。由于水文計算和水動力計算采用不同的時間步長,在一個水文計算時間步長內(nèi),需要進行多次水動力計算。兩個區(qū)域通過動邊界進行空間連接,需要根據(jù)動邊界兩側(cè)的流態(tài)、特征波方向等確定動邊界移動方向及邊界界面的水量、動量通量。
圖4 DBCM 模型構(gòu)成框圖
2.3.4 水文與水動力過程回饋機制分析及耦合模型發(fā)展歷程 圖5 為水文與水動力耦合模型的發(fā)展歷程示意圖,水文模型由簡單的經(jīng)驗公式發(fā)展到分布式水文模型,在流域產(chǎn)匯流計算方面不可缺少,但是單一的水文模型不考慮水動力過程,不能描述洪水演進的動力特征;水動力模型由謝才公式發(fā)展到流域全動力模型,在城市洪澇模擬和預報方面發(fā)揮了重要作用。單一的水動力模型,包括二維淺水方程,在城市洪澇模擬和預報方面具有精度高的優(yōu)勢,但是推廣到流域范圍的洪澇模擬和預報,還需深入討論數(shù)值格式穩(wěn)定性和計算效率問題。圖5 中右側(cè)第二列,水文與2D 水動力動態(tài)單向耦合, 包含兩種類型:(1)水文產(chǎn)流模型與二維(包括一維)水動力動態(tài)單向耦合,因匯流計算依賴二維淺水方程,屬于流域二維全動力洪澇模型。(2)水文產(chǎn)匯流模型首先與一維水動力模型直接動態(tài)單向耦合,然后是一維水動力模型與二維水動力模型動態(tài)雙向耦合,可以視為水文產(chǎn)匯流模型與二維水動力模型的間接動態(tài)耦合。圖5 中右側(cè)第一列,水文與2D 水動力動態(tài)雙向耦合模型,是分布式水文產(chǎn)匯流模型與局部二維(包括一維)水動力模型直接耦合,也可稱為DBCM。BDCM 是由清華大學于2019年提出的,詳細參考文獻[13,61]。DBCM 的構(gòu)建有如下兩個特點:(1)水文產(chǎn)匯流模型與二維水動力的數(shù)值離散采用時間顯格式,避免求解聯(lián)立方程所花費的過多的計算時間;(2)提出的水文產(chǎn)匯流模型與二維水動力模型之間的動態(tài)雙向耦合方法,避免在兩個模型連接處設定內(nèi)邊界條件或通過固定內(nèi)邊界位置來實現(xiàn)水文產(chǎn)匯流模型與二維水動力模型的耦合。對于山區(qū)流域的洪澇模擬和預報問題,在揭示水文過程與水動力過程的動態(tài)互饋機制的基礎上,充分利用水文模型和水動力模型的特點,建立水文與水動力的動態(tài)雙向耦合模型值得深入開展研究,在提高流域洪澇模擬和預報精度的同時,達到模擬精度與計算效率的優(yōu)化配置是今后有待深入探討的方向。
圖5 水文與水動力耦合模型發(fā)展過程示意
流域洪澇發(fā)生和發(fā)展的水文過程與水動力過程在時間是同時進行的,在空間上是緊密聯(lián)系的,水文信息與水動力信息存在動態(tài)互饋關(guān)系。如圖6(a)中所示,常用的水文與水動力的串聯(lián)耦合模型將水文過程與水動力過程在時間上分開進行;在空間上,水文計算的匯流水量結(jié)果作為水動力過程的邊界點,邊界點流量數(shù)值集中了河道兩岸的匯流水量,常常會夸大河道上游入口的流量,流量邊界點位置選定具有人為因素,如果邊界條件加載的不合理,串聯(lián)耦合模型的模擬結(jié)果將會存在一定的誤差。如圖6(b)所示,水文與一維水動力模型的動態(tài)單向耦合(如MIKE SHE 與MIKE 11 的動態(tài)單向耦合),全流域采用分布式水文模型計算產(chǎn)匯流流量,產(chǎn)匯流水量作為邊界條件或源項加載到一維河道上游或側(cè)向有限個節(jié)點上,計算網(wǎng)格一旦確定,邊界點位置不再變化,水動力過程不對水文過程產(chǎn)生影響,一維水動力模型功能限制了模型的應用范圍。如圖6(c)所示,已有的水文產(chǎn)流模型與二維水動力模型動態(tài)單向耦合,如TRIBS-OFM、FRAS、Hydroinfo,克服了串聯(lián)耦合及水文與一維水動力動態(tài)單向耦合模型的不足,有效地提高了模型精度。這類水文與二維水動力的動態(tài)單向耦合模型中,水動力過程不對水文過程產(chǎn)生影響,如果推廣到流域范圍的洪澇模擬和預報,需要解決數(shù)值格式穩(wěn)定性和計算效率問題。圖6(d)為水文與二維水動力的直接動態(tài)雙向耦合模型(DBCM),流域水文計算與水動力計算在時間上同時進行;在空間上,水文計算產(chǎn)匯流水量結(jié)果作為水動力計算的邊界條件,僅在積水等局部區(qū)域進行二維水動力計算,水文計算結(jié)果按實際發(fā)生位置動態(tài)加載到淹沒區(qū)的周圍,淹沒區(qū)與水文產(chǎn)匯流區(qū)以動邊界動態(tài)連接,動邊界位置隨時間變化。圖6(d)為水文產(chǎn)匯流模型與二維水動力的動態(tài)雙向耦合模型,即DBCM 模型示意圖,該模型僅在局部區(qū)間進行二維水動力計算,實現(xiàn)水文產(chǎn)匯流模型與二維(包括一維)水動力模型的直接動態(tài)雙向耦合。如果在全流域采用二維淺水模型進行計算,此時DBCM 等價于圖6(c)的水文產(chǎn)流模型與二維水動力模型的動態(tài)單向耦合。反過來,圖6(c)的水文與水動力耦合模型不能包含DBCM 的功能,因為圖6(c)中的水文模型,僅進行產(chǎn)流計算,匯流計算需要通過淺水方程實現(xiàn),不能開展局部二維水動力計算,需要在全流域進行二維水動力計算,相當于流域全動力洪澇模型。DBCM 模型能更合理地反映流域水文過程與水動力過程的動態(tài)雙向互饋機制,動態(tài)雙向耦合模型在流域范圍的洪澇模擬和預報中將具有發(fā)展?jié)摿Α?/p>
圖6 水文與水動力的不同耦合示意圖
對于洪澇災害模擬和預報,水文與水動力耦合模型具有優(yōu)勢。水文與水動力耦合模型理論需要反映水文過程與水動力過程的互饋機制,水文過程驅(qū)動水動力過程,反過來水動力過程影響水文的匯流條件??紤]水文與水動力這兩個物理過程的相互影響,國內(nèi)外學者提出了多種的耦合方法。
水文與水動力的串聯(lián)耦合模型、動態(tài)耦合模型各有其優(yōu)勢,應該根據(jù)實際要解決的問題選取合理的耦合模型。分布式水文模型與二維水動力模型的動態(tài)雙向耦合更反映水文產(chǎn)匯流過程與水動力洪水演進過程的互饋機制,模型預期具有良好的精度,在流域洪澇災害預報中具有應用前景。
水文模型與二維、一維水動力模型的耦合方式可以分成如下三種模式:耦合模式1,基于水動力模型的動態(tài)耦合方式;耦合模式2,水文模型與一維水動力模型直接耦合及與局部區(qū)域二維水動力的間接耦合;耦合模式3,水文模型直接與二維和一維水動力的動態(tài)耦合。三種耦合模式如圖7 所示。
圖7 耦合模式
耦合模式1:在連續(xù)方程中考慮降雨、蒸發(fā)、入滲等水文產(chǎn)流過程,通過求解淺水方程的動量方程計算匯流及水動力過程。耦合模式1 在全流域劃分差分或控制體網(wǎng)格,進行高分辨率的二維水動力計算,可獲得更加準確的排水管網(wǎng)入口處水位和流速。耦合模式1 的代表性軟件,包括中國水利水電科學研究院的FRAS[35,63]、西安理工學研發(fā)的GAST[39]、大連理工大學研發(fā)的Hydroinfo[25]、密西根大學研發(fā)的TRIBS-OFM[62]、亞利桑那大學研發(fā)的CHRE2D[36-37]。中國水利水電科學研究院的FRAS 在二維網(wǎng)格內(nèi)通過特殊通道法,實現(xiàn)二維水動力與一維水動力的動態(tài)連接,實現(xiàn)了河道、排水口附近的精細水動力模擬;西安理工大學的GAST模型,求解二維淺水方程,通過GPU技術(shù)提高計算效率[39];大連理工大學的Hydroinfo[25]通過非構(gòu)造網(wǎng)格等手段,提高局部區(qū)域的計算精度,并提高整體計算效率。耦合模式的很多模型,均在城市及流域洪澇防治工程中得到廣泛應用。
耦合模式2:由兩部分組成,第一部分是分布式水文模型直接與河道、管網(wǎng)的一維水動力的動態(tài)耦合,與SWMM 或Mike Urban有類似之處,也有不少改進;第二部分是淹沒區(qū)域的局部二維水動力模型與一維水動力的動態(tài)雙向耦合,與MIKE Flood有類似之處,但有所改進。耦合模式2的二維水動力計算僅在局部區(qū)域進行,在局部區(qū)域進行網(wǎng)格加密,可以獲得排水口入口處的更加準確的水位和流速信息。局部區(qū)域精細網(wǎng)格的二維水動力計算,可以節(jié)省計算時間,提高計算效率。中國水利水電科學研究院的CUHHM[40-41]、華南理工大學的IHUM[42-43]等屬于耦合模式2,該模式已經(jīng)在城市洪澇防治工程應用。
耦合模式3:是分布式水文模型與二維及一維水動力模型的直接耦合,屬于動態(tài)雙向耦合。將流域空間分成非淹沒區(qū)和淹沒區(qū),僅在非淹沒區(qū)采用粗網(wǎng)格求解二維擴散波方程,計算流域的產(chǎn)流和匯流;在淹沒區(qū)等局部區(qū)域采用精細網(wǎng)格求解二維淺水方程,淹沒區(qū)包括河湖、漫灘、低洼積水區(qū)及排水管網(wǎng)密集區(qū)域;模型預期具有良好的數(shù)值精度。由于僅在局部區(qū)域(淹沒區(qū))進行二維水動力計算,耦合模式3 預期具有良好的計算效率。清華大學提出的DBCM[61]模式屬于耦合模式3,該耦合模式提出的時間較短,模型尚在完善之中。
流域洪澇災害發(fā)生頻發(fā),研制并應用洪澇災害模擬及預報模型,對于減緩洪澇災害影響具有現(xiàn)實意義。在評價工程措施對減緩洪澇災害風險的效果時,國內(nèi)自主開發(fā)的模型和國外的商業(yè)模型都具有其優(yōu)點。國外的Mike 系列、SWMM、SWAT、HEC-RAS 等模型,在高校和工程設計院等部門經(jīng)常使用,在國內(nèi)有較好的應用和銷售市場。由于國外的商業(yè)模型模塊化比較強,不同功能在同一模型中不能同時兼顧,對國內(nèi)自主研發(fā)模型不利。國內(nèi)自主研發(fā)的模型在模型功能、計算精度、計算效 率和格式穩(wěn)定等方面都具有獨到之處和可持續(xù)開發(fā)的 潛力,如FRAS[35,63],GAST[39]、Hydroinfo[25]、CUHHM[40-41]、IHUM[42-43]、FRAS[35,63]、DBCM[61]等模型,各種耦合均勻有應用案例,模型的先進性已經(jīng)得到國內(nèi)外同行的良好評價。揭示水文過程與水動力過程的真實相互影響機制,在保證數(shù)值格式穩(wěn)定性的同時,實現(xiàn)水文與水動力耦合模型的計算精度和計算效率的優(yōu)化配置,增加模型的功能,將模型應用于流域范圍的洪澇災害及水環(huán)境綜合整治,是今后需要深入研究的方向之一。
將水文模型和水動力學模型作為一個整體統(tǒng)籌考慮,這種耦合模式稱為緊密耦合或雙向耦合[61]。緊密耦合在機理上最為完善,通過緊密耦合實現(xiàn)水文、水動力、水質(zhì)等多過程的交互,進而實現(xiàn)模擬復雜的洪澇過程和水環(huán)境,是洪澇模型的重要發(fā)展趨勢[64]。
本文介紹了流域洪澇模擬和預報模型的研究進展,包括水文模型、水動力模型、水文與水動力耦合模型。在實際應用時應該根據(jù)流域或城市洪澇問題的實際,選取合適的模型。單一的水文模型具有計算簡便、模擬范圍大、計算效率高等優(yōu)點,在流域洪澇模擬與預報中已經(jīng)發(fā)揮了重要作用。在流域范圍較小或在城市洪澇的模擬與預報中,單一的水動力模型具有較高的空間分辨率,可以反映洪水與建筑物耦合作用的動力過程,能為城市排水管網(wǎng)進口位置提供更加合理的水深和流速等信息。通過加大二維區(qū)域的網(wǎng)格尺寸、一維水動力計算與二維水動力計算耦合、在水深較淺區(qū)域使用“薄層水”假設、二維水動力與二維擴散波模型切換等方法,可以改進數(shù)值格式的穩(wěn)定性,提高計算效率。
水文與水動力的耦合模型可分為串聯(lián)耦合、動態(tài)單向耦合、動態(tài)雙向耦合。水文與水動力的串聯(lián)耦合模型因使用簡便,應用靈活等特點,在洪澇災害治理中工程中經(jīng)常使用。串聯(lián)耦合模型適合淹沒區(qū)范圍較小的洪水演進分析,在水動力計算的邊界點位置選取和流量邊界條件選取上應該注意與物理過程相符合等問題,盡量避免邊界點位置固定和流量邊界條件帶來的誤差。水文與一維水動力模型的動態(tài)單向耦合模型,因水動力模型是一維的,模型的應用范圍受限。水文與二維、一維水動力的動態(tài)耦合模型應用范圍較廣,已有的三種耦合模式如圖7 所示,各種耦合模式均具有各自的優(yōu)勢,可以根據(jù)流域或城市洪澇的實際情況選取和使用合理的耦合模式。由于耦合模式3 的開發(fā)時間相對較短,模型的功能還處于完善過程中。流域洪澇模型的計算精度與計算效率是一對矛盾,往往不能同時兼顧。常規(guī)的洪澇預測模型常常是計算精度的提高,以花費更多的計算時間及格式穩(wěn)定性為代價,采用并行計算技術(shù)也是減少模型運行時間的常用方法。在取得良好的計算精度的同時,能夠提高計算效率,需要在模型理論上進行創(chuàng)新。CUHHM、IHUM、IFMS Urban[65]等正是朝著這個方向開展工作,已經(jīng)研制出具有獨立知識產(chǎn)權(quán)的水文產(chǎn)匯流模型與二維水動力的間接動態(tài)耦合模型。于2019年提出的DBCM,是水文產(chǎn)匯流模型與二維(包括一維)水動力模型的直接動態(tài)雙向耦合,在模型理論上能更加真實地反映水文過程與水動力過程的互饋機制,同時避免水文產(chǎn)匯流模型與二維水動力模型之間互為邊界條件的算法、盡量減少采用二維隱式格式(聯(lián)立求解大型方程組)來達到保證格式穩(wěn)定的目的。通過改進水文與水動力模型的耦合方式,研制計算精度高、數(shù)值格式穩(wěn)定、計算效率高的流域洪澇模擬及預測模型,為流域洪澇模擬和預報、智慧水務建設提供技術(shù)支撐。