基于水文水動力模型的平原地表水資源研究

高 慶，劉蕓辰，李 根，衣 鵬

(1.河海大學(xué)，江蘇 南京 210024；2.泰州市水利局海陵分局，江蘇 泰州 225300)

1 研究背景與意義

水資源作為自然環(huán)境至關(guān)重要的組成部分，是人類生存發(fā)展必不可少的自然資源，更是衡量一個國家綜合實力的戰(zhàn)略性經(jīng)濟資源[1]。地表水資源作為水資源的重要組成部分，和人類關(guān)系最為密切，受人類影響最顯著；同時，作為研究區(qū)水資源合理開發(fā)利用、優(yōu)化規(guī)劃、最優(yōu)調(diào)度的依據(jù)，開展地表水資源研究對于優(yōu)化水資源管理具有極重要的現(xiàn)實價值[2]。因此，地表水資源研究一直受到廣泛的關(guān)注。

平原河網(wǎng)區(qū)地平坡緩，河流縱橫，大都是密度高且城市化率較高的地區(qū)，水資源對經(jīng)濟發(fā)展的制約作用極其突出[3]，開展地表水資源研究，科學(xué)制定水資源規(guī)劃迫在眉睫。然而，平原河網(wǎng)區(qū)水流流速緩慢，流向、流態(tài)隨機性大，水利工程密布，水量交換關(guān)系錯雜，且大多缺乏水文站點，無法依據(jù)實測徑流資料研究地表水資源，多數(shù)采用還原法計算。然而，該方法任意性較大，不能較好反映真實的地表徑流情況[4- 6]，加之受氣候變化和人類活動對下墊面的巨大影響，水循環(huán)過程發(fā)生變化，進而影響了地表水的產(chǎn)匯流機制，造成還原計算法的困難更大[7]。目前，水文水動力耦合模型在國內(nèi)外已得到成功應(yīng)用，且應(yīng)用范圍越來越廣泛。李致家等[8]將新安江模型與一維、二維水力學(xué)模型進行耦合，模擬研究南四湖上級湖的流量演進和二級壩水利樞紐的調(diào)度。從水位模擬值的精度來看，該耦合模型能夠較準確地模擬出湖泊水位演進過程。Thompso等[9]通過耦合MIKESHE模型和水力學(xué)模型MIKE11，深入研究了英國東南部某一濕地系統(tǒng)的水文過程。模型的模擬結(jié)果和實際觀測值擬合度較高，可見該耦合模型完善了MIKESHE模型在濕地系統(tǒng)進行水文模擬應(yīng)用的不足。

本次研究區(qū)域選取江蘇省泰州市海陵區(qū)，該區(qū)域地處長江、淮河流域交匯處，屬于典型的平原河網(wǎng)區(qū)，是南水北調(diào)重要輸水通道之一，具有豐富的過境水資源。本研究基于河海大學(xué)開發(fā)的數(shù)字流域模型，建立水文水動力耦合模型。該模型可以在缺乏徑流資料的情況下，充分利用研究區(qū)域的水文氣象資料、水系資料和下墊面資料概化模型，較準確地模擬出區(qū)域的徑流狀況，可用于研究本地地表水資源和過境水資源，并為之后的平原河網(wǎng)區(qū)水資源量研究提供參考。

2 研究區(qū)域概況

海陵區(qū)位于江蘇省中部、長江中下游北岸，西臨揚州，東靠南通，與蘇州、無錫、常州隔江相望，是泰州市的經(jīng)濟政治文化中心，區(qū)域范圍為32.45°N～32.57°N，119.80°E～119.98°E(詳見圖1)。該區(qū)域總面積為218.02 km2，其中水面面積38.08 km2，約占總面積的17.5%，區(qū)域內(nèi)主要有鹵汀河、南官河、泰東河、引江河、新通揚運河和老通揚運河6條骨干河道。海陵區(qū)地形上為平原，屬于江淮沖積平原，地勢平坦，西高東低、南高北低，南部多平地，北部河網(wǎng)密布。海陵區(qū)地處長江水系與淮河水系的交匯處，河流縱橫，水網(wǎng)密布，境內(nèi)水域廣泛。因地形、水系不同，以橫貫中部的老328國道為界，劃分為兩個水資源分區(qū)：北部為淮河流域的里下河腹部區(qū)(俗稱“下河地區(qū)”)，南部為長江流域的通南高沙土區(qū)(俗稱“上河地區(qū)”)(見圖1)。正常情況下，上河水位高于下河，上河歷年平均水位2.20 m，下河歷年平均水位1.30 m，平均水位差為0.90 m。

3 研究方法和模型

海陵區(qū)水文水動力耦合模型包括產(chǎn)匯流模型和河網(wǎng)水動力模型。首先，通過平原區(qū)產(chǎn)匯流模型，模擬海陵區(qū)4種下墊面的降雨徑流關(guān)系和凈雨的坡面匯流過程；然后，通過河網(wǎng)水系概化、水利工程概化和下墊面產(chǎn)流分配，建立河網(wǎng)水動力學(xué)模型；接著，通過耦合這兩個模型，構(gòu)建海陵區(qū)水文水動力耦合模型，來模擬水流在河網(wǎng)中的運動過程，并在此基礎(chǔ)上，創(chuàng)建來水組成模塊，為分析過境水資源來水組成提供依據(jù)；最后，在海陵區(qū)耦合模型構(gòu)建及參數(shù)設(shè)定的基礎(chǔ)上率定產(chǎn)匯流、河道糙率等模型參數(shù)，并基于率定好的耦合模型，開展模擬計算，驗證該耦合模型的在海陵的有效性和適用性。

3.1 產(chǎn)匯流模型

不同下墊面所對應(yīng)的產(chǎn)流規(guī)律差異較大。為此，本研究基于劃分的下墊面類型：水面、旱地、水田和建設(shè)用地這四類分別進行產(chǎn)流模擬，計算對應(yīng)的產(chǎn)流量[10-11]。

匯流過程從空間上分為坡面匯流和河網(wǎng)匯流兩個階段。其中，坡面匯流模型通過由周邊河道圍成的河網(wǎng)多邊形進行構(gòu)建，并采用平原河網(wǎng)區(qū)特有的平原區(qū)坡面匯流單位線來進行模擬，而河網(wǎng)匯流由河網(wǎng)一維水動力模型來模擬。

平原區(qū)的產(chǎn)流經(jīng)過地面匯流才能進入周邊河道，關(guān)于平原區(qū)坡面匯流單位線的確定目前尚無成熟的理論和計算方法。同時，考慮到水文模型的計算結(jié)果存在多種誤差的積累，包括降雨產(chǎn)流模擬誤差、匯流誤差、河網(wǎng)概化誤差、河網(wǎng)糙率和水利工程調(diào)度模擬誤差等，因此，無法通過率定來選擇合適的匯流過程。故，本次研究移用應(yīng)用成熟的里下河水文水動力模型采用的匯流過程，根據(jù)不同計算分區(qū)需要，按不同區(qū)域的匯流路徑及網(wǎng)格，采用坡面匯流流速計算生成匯流單位線匯入河網(wǎng)[12-15]。

3.2 河網(wǎng)水動力模型

3.2.1 模型原理

根據(jù)對平原河網(wǎng)區(qū)水流運動規(guī)律的分析，可將其概化為河道水流運動的一維模擬和堰閘泵站單元模擬兩塊。針對各要素的特點，分別采用合適的數(shù)值求解方法，將各部分的求解有機結(jié)合，形成全耦合求解的區(qū)域動力模擬模型。

(1)一維水動力模型。描述河道一維水流運動的圣維南方程組為

(1)

(2)

式中，Z、Q分別為河道的斷面水位和流量；B、A分別為河寬和河道過水面積；q為旁側(cè)入流；a為動量校正系數(shù)，表示河道斷面流速分布均勻性的系數(shù)；K為流量模數(shù)，表示河道的實際過流能力；Vx為旁側(cè)入流在水流方向上的流速分量，一般可以近似為零。本模型采用有限差分隱式格式的簡化格式，即四點線性隱式格式進行數(shù)值離散[16]，求解圣維南方程組。

(2)堰閘泵站單元模擬。假定水工建筑物的上下游各有一個節(jié)點，建筑物與概化河道經(jīng)該節(jié)點進行連接(見圖2)，而節(jié)點間的水位和流量，主要是根據(jù)堰流公式、水工建筑物運行方式進行計算確定；且在模型中，水工建筑物的運行方式根據(jù)相對應(yīng)的邏輯控制條件進行概化。

圖2 水利工程概化示意

以水閘為例，閘門開啟時的過閘流量以寬頂堰公式為依據(jù)開展計算，寬頂堰上的水流流態(tài)主要包括自由出流和淹沒出流，分別采用對應(yīng)公式計算。即

H0=Z1-Zd，hS=Z2-Zd

(3)

當流態(tài)為自由出流時，為

(4)

當流態(tài)為淹沒出流時，為

(5)

式中，Zd為堰頂高程；Z1、Z2分別為堰上、下水位；m為自由出流系數(shù)，一般在0.325～0.385之間；B為堰寬；φm為淹沒出流系數(shù)，本次取1.0。

3.2.2 模型概化

根據(jù)模型中河網(wǎng)概化的原則和要求，結(jié)合海陵區(qū)實際情況和水循環(huán)模擬需要，海陵區(qū)水文水動力模型概化主要包括河網(wǎng)水系概化和水利工程概化。

河網(wǎng)水系的概化實行物理河網(wǎng)概化的原則。即，不采用小河道合并的概化模式，概化河道全部對應(yīng)實際河道，原則是概化河網(wǎng)的過流能力能夠反映實際情況；最后，對所有河道的連接性、連通性及拓撲關(guān)系進行檢查，形成海陵區(qū)河網(wǎng)水系圖最后成果(見圖3)。概化河道在模型中采用一維水動力模型計算，區(qū)域內(nèi)其余小支流按照模型模擬需要及海陵區(qū)水系現(xiàn)狀確定，均作為陸域面上的調(diào)蓄水面處理。

圖3 海陵區(qū)模型概化后河網(wǎng)示意

水利工程(閘、泵站)概化為聯(lián)系模型要素，其過流流量采用水動力學(xué)方法模擬，建筑物的上下游設(shè)節(jié)點，通過節(jié)點與概化河道連接；節(jié)點間的水位及流量取決于堰流公式及建筑物運行方式，運行方式用一系列的邏輯控制條件進行模擬概化。

3.2.3 下墊面產(chǎn)流分配

本次采用的水文水動力耦合模型，為了更準確地反映下墊面的空間分布，采用柵格化方法，把區(qū)域網(wǎng)格和下墊面圖層進行疊加，將下墊面信息數(shù)字化，再借助河網(wǎng)多邊形將產(chǎn)水量分配到河網(wǎng)。

河網(wǎng)多邊形是指概化河道構(gòu)成的封閉區(qū)域，可通過概化的河道自動生成。模型運行時，為了進行數(shù)值計算，將計算區(qū)域劃分成固定大小的網(wǎng)格，通過劃分好的網(wǎng)格把計算分區(qū)和下墊面信息柵格化。選擇河網(wǎng)多邊形中某網(wǎng)格，用如圖4所示的圖案填充的網(wǎng)格，能夠求出該網(wǎng)格中心點和周邊各條概化河道之間的距離s。倘若沒有詳細地形，一般假設(shè)選定網(wǎng)格的產(chǎn)流將匯入距離最近的一條概化河道，也就是采用距離最短的方法來判斷匯流的方向。但是只根據(jù)距離來判斷網(wǎng)格流向，經(jīng)常會導(dǎo)致一些不合理的現(xiàn)象。譬如說，過水能力較小的河流在洪水時會引發(fā)來不及排泄造成水位太高的情況。

圖4 網(wǎng)格和河網(wǎng)多邊形

因此，本模型在網(wǎng)格產(chǎn)流結(jié)果分配時考慮到河道過水能力，選取綜合系數(shù)進行計算，來避免出現(xiàn)這種不合理現(xiàn)象。綜合系數(shù)計算公式為

f=s/AR0.67

(6)

式中，s為計算網(wǎng)格中心點到概化河道的最小距離；A為河道的過水面積；R為河道的水力半徑。

通過計算每一網(wǎng)格距周邊所有河道的綜合系數(shù)，將網(wǎng)格分配到河網(wǎng)多邊形綜合系數(shù)最小的河道，網(wǎng)格屬于的下墊面也隨之分配到對應(yīng)的概化河道。

3.3 模型耦合

平原河網(wǎng)區(qū)的降雨產(chǎn)匯流不僅無法全部直接匯入河道，還反過來影響以后的產(chǎn)匯流。尤其在產(chǎn)流量計算中，經(jīng)常將其與河道匯流獨立開來計算，未考慮河道頂托的影響，而這種影響在平原地區(qū)是不可忽視的[17-19]。

為了全面、真實地模擬區(qū)域產(chǎn)匯流過程，本研究采用水動力學(xué)與水文學(xué)相結(jié)合的方法，建立水文水動力耦合模型[20-23]。將水文模型計算出的出流量輸入到對應(yīng)的河道斷面中，以此來驅(qū)動水動力學(xué)模型，可以取長補短，構(gòu)建精度較高的平原河網(wǎng)區(qū)水文水動力耦合模型，對區(qū)域地表水資源進行精細化研究[24-28]。

總體而言，平原區(qū)產(chǎn)流模型解決區(qū)間產(chǎn)流部分，將匯流模擬計算結(jié)果分別作為圩區(qū)及圩外調(diào)蓄水面處理區(qū)域的入流，通過構(gòu)建的河網(wǎng)多邊形，將其精確匯入河網(wǎng)水動力模型中的每條河道斷面，實現(xiàn)水文模型和水動力模型的耦合，具體模型結(jié)構(gòu)見圖5。

圖5 水文水動力耦合模型結(jié)構(gòu)

水文模型模擬時對下墊面進行網(wǎng)格化處理，利用網(wǎng)格信息確定概化河道的陸域調(diào)蓄寬度，并通過產(chǎn)流模塊計算網(wǎng)格產(chǎn)水量。通過和該網(wǎng)格聯(lián)系的河道，能夠把每個網(wǎng)格對應(yīng)的產(chǎn)水量匯集到概化河道，求和可以得到所有河道的圩內(nèi)、外和陸域面上交換的水量。

然而，該處理方法工作量較大，對每個網(wǎng)格都進行產(chǎn)流計算，會極大地影響模型的運算速度和效率。因此，本次模擬將每個網(wǎng)格對應(yīng)的下墊面，疊加到與之相聯(lián)系的河道上。在此設(shè)置下，模型計算時先疊加各下墊面的網(wǎng)格數(shù)；再以下墊面為基本單元計算產(chǎn)流，同一下墊面類型在產(chǎn)流模塊中只進行一次計算。各網(wǎng)格的產(chǎn)流量可以根據(jù)對應(yīng)的下墊面類型確定，獲得地理分布下的下墊面產(chǎn)流計算結(jié)果，最終得到海陵區(qū)多年平均地表水資源量為13 601萬m3(見表1)。將模型計算結(jié)果和傳統(tǒng)地面分類法計算出的結(jié)果進行比較，模型計算出來的地表水資源量偏大688.3萬m3，誤差為5.3%，尚在合理范圍之內(nèi)。

表1 海陵區(qū)多年平均地表水資源量計算結(jié)果

4 模型率定與驗證

本次主要對海陵區(qū)水文水動力耦合模型進行率定和驗證，使該模型能夠更好地適用于海陵區(qū)水流過程模擬，為進一步開展海陵區(qū)水資源、水環(huán)境方面的相關(guān)研究提供技術(shù)支撐[23，29-30]。結(jié)合海陵區(qū)近年下墊面特征變化，分別選取2009年～2012年和2013年～2016年作為模型率定和驗證年份。

4.1 模型率定

本次主要對2009年～2012年的水位進行率定，率定資料包括實測降雨、蒸發(fā)、水位、流量等資料，在分析海陵區(qū)實際雨水情和站點分布情況的基礎(chǔ)上，選擇泰州(通)站2009年～2012年全年的日平均實測水位資料，作為模型率定參照值。

由于目前無法對產(chǎn)流模型進行率定，本次主要對照《江蘇省水文手冊》和1984年《江蘇省暴雨洪水圖集》中的次降雨～徑流關(guān)系，對產(chǎn)流進行復(fù)核。同時參照泰州市水資源調(diào)查評價有關(guān)成果，調(diào)整產(chǎn)流模型參數(shù)，使計算結(jié)果基本合理，徑流系數(shù)符合平原區(qū)降雨徑流關(guān)系。

從泰州(通)站水位率定成果看，2009年～2012年模型計算的代表站日均水位過程與實測過程總體擬合較好，主要是5月份、6月份的擬合結(jié)果較差。由于此時為水稻田生長的返青和分蘗期，區(qū)域的實際用水情況與灌溉制度可能存在差異；同時，區(qū)域養(yǎng)殖用水情況復(fù)雜，受水田灌溉用水過程和河網(wǎng)實際取水過程對水位模擬造成影響，導(dǎo)致擬合效果較差。選取水位代表站泰州(通)站見圖6。

圖6 2009年～2012年泰州(通)站計算和實測水位過程對比

對率定結(jié)果進行統(tǒng)計分析，得到率定期模型模擬的納什效率為0.65，相對誤差為15.21%。前者大于0.6，后者小于20%，則表示率定期模型計算值與實測值總體擬合度較好。

4.2 模型驗證

利用率定好的海陵區(qū)水文水動力耦合模型，對2013年～2016年泰州(通)站的實測水位為對比開展驗證，采用的資料主要涵蓋2013年～2016年的實測降雨、蒸發(fā)、水位、流量資料等。

根據(jù)研究區(qū)域海陵區(qū)水雨情特點，本次研究在2013年～2016年中各選取一場有代表性的洪水進行模擬計算，從而驗證所搭建的水文水動力耦合模型是否可以較為準確地體現(xiàn)流域產(chǎn)匯流特性、能否真實描述河道的水流情況及下墊面概化、河道斷面概化和產(chǎn)匯流參數(shù)的設(shè)置是否合理。

本次驗證選擇了泰州(通)站的日水位資料，繪制了泰州(通)站近兩年的觀測水位和計算水位的過程(見圖7)。從泰州(通)站的水位驗證成果來看，水位過程與實測過程總體擬合較好。對驗證結(jié)果進行分析，得到驗證期模型模擬的納什效率為0.77，相對誤差為8.63%，前者大于0.6，后者小于20%，進一步反映了驗證期模型計算值與實測值總體擬合度較佳。

圖7 泰州(通)站計算和實測水位過程對比

5 結(jié) 語

平原河網(wǎng)區(qū)由于水量關(guān)系交換錯雜、缺乏站點資料等原因，地表水資源的研究一直是個難題。傳統(tǒng)水文多采用還原法計算地表水資源，然而該方法依舊不能很好地反映真實情況。本文基于成熟的數(shù)字流域模型，將水文模型和水力學(xué)模型二者相結(jié)合，構(gòu)建了海陵區(qū)水文水動力耦合模型，對平原河網(wǎng)區(qū)的地表水資源進行模擬研究。經(jīng)過率定驗證，得到率定期模型模擬的納什效率為0.65，相對誤差為15.21%，驗證期模型模擬的納什效率為0.77，相對誤差為8.63%；表明擬合較好。最終，得到的海陵區(qū)多年平均地表水資源量為13 601萬m3；同時，為海陵區(qū)科學(xué)制定水資源規(guī)劃，確保區(qū)域經(jīng)濟可持續(xù)發(fā)展提供基礎(chǔ)數(shù)據(jù)支持。

本次研究基于海陵區(qū)典型平原河網(wǎng)區(qū)，改善了傳統(tǒng)研究方法受有限徑流資料的限制狀況，創(chuàng)新性地將平原區(qū)產(chǎn)匯流模型和河網(wǎng)水動力模型二者結(jié)合，綜合考慮區(qū)域降雨、產(chǎn)流及河網(wǎng)匯流運動特性，為進行地表水資源量模擬研究提供了較為準確合理的模型支撐。