水稻灌區(qū)基于多過程耦合的分布式水轉(zhuǎn)化模型

王 璞，盧 炤，霍再林

（中國農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，北京 100083）

1 研究背景

水稻種植面積占我國糧食作物種植總面積的27%，稻田灌溉用水量占我國農(nóng)業(yè)總用水量65%以上，但其灌溉水利用系數(shù)僅在0.5 到0.6 之間［1-3］。由于灌溉定額大，水稻灌區(qū)灌溉水垂向及水平向轉(zhuǎn)化頻繁，灌區(qū)水轉(zhuǎn)化模型是灌區(qū)水資源科學(xué)管理的基礎(chǔ)［4-6］。

目前針對灌區(qū)水轉(zhuǎn)化模擬的方法主要包括：（1）基于陸地水文模型的灌區(qū)水轉(zhuǎn)化過程模擬。SWAT 模型已被廣泛應(yīng)用于灌區(qū)水轉(zhuǎn)化、面源污染、農(nóng)業(yè)生產(chǎn)力等方面的模擬［7-9］。特別是近年來，SWAT 模型與MODFLOW 模型結(jié)合，實現(xiàn)了灌區(qū)地表水-土壤水-地下水的耦合模擬［7］。（2）基于農(nóng)田水文模型的分布式模擬。一般是對農(nóng)田水文模型進行率定驗證后，將其拓展到灌區(qū)尺度進行模擬［10-11］。該方法可對農(nóng)田垂向水轉(zhuǎn)化及作物生長過程進行詳細(xì)描述，但由于缺乏對地下水水平運動的模擬以及灌溉-排水-地下水等水轉(zhuǎn)化過程的表征，無法應(yīng)用于不同情景下灌區(qū)水轉(zhuǎn)化的模擬。（3）基于灌區(qū)分布式水平衡的模擬。一般是對灌域水平衡進行分析計算，進而將其拓展至灌區(qū)尺度。如岳衛(wèi)峰等［12］針對內(nèi)蒙古河套灌區(qū)義長灌域基于四水轉(zhuǎn)化關(guān)系進行了灌域水均衡分析，該方法盡管可獲得灌區(qū)整體水平衡，但其時空分辨率低，無法實現(xiàn)灌區(qū)尺度詳細(xì)的水轉(zhuǎn)化過程模擬［12-13］。

綜上所述，現(xiàn)有的灌區(qū)水轉(zhuǎn)化模型尚未充分體現(xiàn)灌區(qū)在灌溉與排水驅(qū)動下，灌溉水、土壤水及地下水在垂向與水平向頻繁轉(zhuǎn)化的特點。基于此，本研究旨在發(fā)展灌區(qū)尺度基于灌排-農(nóng)田水文-地下水過程耦合的水轉(zhuǎn)化模型，并將其應(yīng)用于黑龍江省綏化市慶安縣和平灌區(qū)，為灌區(qū)尺度水平衡分析提供有效方法。

2 研究區(qū)域及水文監(jiān)測

2.1 研究區(qū)概況本研究以黑龍江省和平灌區(qū)為研究區(qū)，該灌區(qū)地處松嫩平原和小興安嶺余脈的交匯地帶，屬呼蘭河流域中上游（圖1）。和平灌區(qū)地理坐標(biāo)東經(jīng)127°21′至127°45′，北緯46°5′至46°92′；屬寒溫帶大陸性季風(fēng)氣候，年均日照時數(shù)2599 h，年均氣溫1.69 ℃，平均年降雨量577 mm。和平灌區(qū)控制面積約170 km2，其中耕地占比90.5%，耕地內(nèi)97%種植作物為水稻。灌區(qū)分為井灌區(qū)與渠灌區(qū)，其中渠灌區(qū)地表引水來源主要包括呼蘭河、拉林清河與安邦河。灌區(qū)內(nèi)現(xiàn)有干渠1 條，主要支渠20 條，主要排水溝道9 條。灌區(qū)和平渠首灌溉期多年平均灌溉引水量5.74 億m3。灌區(qū)地下水埋深由上游的2 ～ 3 m 向下游井灌區(qū)13 ～ 18 m 逐漸增大。

圖1 研究區(qū)位置

2.2 灌區(qū)水文監(jiān)測與數(shù)據(jù)收集灌區(qū)水文監(jiān)測數(shù)據(jù)包括2018—2019年河流水位、干支渠流量、溝道流量及地下水位，用于模型的率定與驗證。河流水位監(jiān)測斷面包括李山屯渠首、鄭文舉渠首、安邦河渠首，于河流水位監(jiān)測斷面布設(shè)水尺并標(biāo)定黃海高程；渠道及排水溝道數(shù)據(jù)除干渠流量外監(jiān)測灌區(qū)上、中、下游3 條典型支渠（二支渠、六支渠、慶勝支渠）渠首流量，2 條典型排水溝道（一支排、七支排）流量，每日監(jiān)測一次。渠道流量數(shù)據(jù)使用水尺監(jiān)測數(shù)據(jù)結(jié)合水位流量關(guān)系曲線計算，溝道流量數(shù)據(jù)采用多普勒超聲流量計測量。由于灌區(qū)內(nèi)支渠眾多且未實現(xiàn)數(shù)字化監(jiān)測，本研究難以監(jiān)測全部渠道每日流量。在對灌區(qū)調(diào)研后認(rèn)為灌區(qū)上中下游對應(yīng)的支渠流量變化趨勢基本一致，研究區(qū)域內(nèi)其他支渠的流量使用已知流量的3 條支渠結(jié)合渠道設(shè)計流量推求。為保證一定的準(zhǔn)確性使用逐日觀測的干渠流量進行對比驗證，如支渠總流量與干渠流量不符則考慮閘門開閉、各渠道控制面積等因素加以調(diào)整。此外，研究區(qū)布設(shè)13 眼地下水觀測井，采用Hobo U20-001-01 自計水位計逐日6 次測定各觀測井地下水位（圖2）。模型運行所需氣象數(shù)據(jù)收集于毗鄰灌區(qū)的鐵力氣象站。同時基于灌區(qū)內(nèi)土壤質(zhì)地調(diào)查（64 處土壤剖面），獲取作物根區(qū)土壤條件及水分運動參數(shù)分布，用于灌區(qū)土壤水分運動模塊（圖3）。

圖2 研究區(qū)DEM 數(shù)據(jù)及水尺、觀測井布設(shè)位置示意

圖3 研究區(qū)土壤參數(shù)空間分布

3 灌區(qū)分布式水轉(zhuǎn)化模型

3.1 模型框架

3.1.1 模型結(jié)構(gòu) 該模型包括灌溉、排水、農(nóng)田作物生長與水轉(zhuǎn)化、地下水運動4 個模塊，基于Mat?lab R2018b 與ArcGIS 10.7 軟件實現(xiàn)灌區(qū)供耗排耦合過程定量模擬。其中灌溉與排水模塊是在對灌區(qū)灌溉渠道及排水溝道進行數(shù)字化的基礎(chǔ)上，進行灌溉渠道入滲損失、各水文響應(yīng)單元灌溉水量、井灌區(qū)抽水量以及排水單元排水量的計算；農(nóng)田作物生長及水循環(huán)模塊旨在耦合模擬農(nóng)田一維水轉(zhuǎn)化及作物生長過程；地下水運動模塊主要模擬計算地下水的垂向及水平向運動。模型將水文響應(yīng)單元劃分為均一網(wǎng)格，灌溉模塊將渠道來水分配至水文響應(yīng)單元，農(nóng)田作物生長與水轉(zhuǎn)化模塊對各水文響應(yīng)單元作物生長及水分轉(zhuǎn)化過程進行迭代計算，輸出作物生長指標(biāo)及水分分量運算結(jié)果，為排水模塊及地下水模塊提供地表產(chǎn)流量及地下水補給量。各模塊進行時空耦合，最終實現(xiàn)灌區(qū)供耗排過程的耦合模擬。模型具體結(jié)構(gòu)及關(guān)系見圖4。

圖4 模型整體結(jié)構(gòu)圖

3.1.2 邊界條件 依據(jù)研究區(qū)水文地質(zhì)條件及地表河流分布，西北以呼蘭河為定水頭邊界；東北以安邦河為定水頭邊界；南部以柯木克河為定水頭邊界；西部為依據(jù)呼蘭河與柯木克河插值確定的定水頭邊界；東南部為定流量邊界（圖5），邊界條件中水頭及流量均采用實測值。

圖5 研究區(qū)地下水隔水頂板高程及邊界條件分布示意

3.1.3 水文響應(yīng)單元劃分 基于ArcGIS 10.7 軟件對研究區(qū)進行水文響應(yīng)單元劃分，水文響應(yīng)單元大小為250 m×250 m，灌區(qū)內(nèi)有效水文響應(yīng)單元共2715 個。

3.2 各過程定量表征方法

3.2.1 農(nóng)田水轉(zhuǎn)化與作物生長過程 模型中作物生長過程采用EPIC（Erosion Productivity Impact Calcu?lator）模型進行描述［14-15］。為便于準(zhǔn)確量化稻田灌溉-土壤水-地下水轉(zhuǎn)化過程，模型將稻田在垂向分為5 層，由上而下分別為儲水區(qū)、泥漿層、傳輸層、飽和層與不透水層（圖6）。儲水區(qū)為稻田水層，當(dāng)儲水區(qū)高度超過田埂高度時會產(chǎn)生田間排水；泥漿層為犁底層之上的土壤；傳輸層為犁底層之下，地下水位之上的土壤，一般透水性強；飽和層為地下水所在，飽和層與傳輸層相鄰界面為地下水位；不透水層為隔水底板。

圖6 土壤垂向分層示意

水稻生育期內(nèi)，稻田土壤水與田面水層逐日均有補給與消耗。在稻田存在水層與無水層時，稻田水量平衡總方程分別可表示為：

稻田存在水層：

稻田不存在水層：

式中：S 為田面水層深度，mm；P 為降雨量，mm；I 為灌溉量，mm；ET 為蒸散發(fā)量，mm；D 為排水量，mm；J 為深層滲漏量，mm；m 為土壤含水量，mm；t 為日序數(shù)。

依據(jù)作物系數(shù)與水分脅迫系數(shù)計算農(nóng)田作物實際耗水量，其中參考作物耗水量ET0的計算采用FAO-56 推薦的Penman-Monteith 公式，作物潛在耗水量ETp依據(jù)參考作物耗水量與作物系數(shù)計算得出，同時依據(jù)植物葉冠生長程度劃分土壤潛在蒸發(fā)Ep與作物潛在蒸騰Tp，分別計算水分脅迫下的實際蒸發(fā)、蒸騰量Ea與Ta，最終獲得農(nóng)田實際蒸散發(fā)量［16-19］。在此過程中所需的作物生長指標(biāo)如葉面積指數(shù)等由EPIC 模型計算提供。稻田深層滲漏量依據(jù)田面水層深度及土壤質(zhì)地等值計算而得。

當(dāng)田面存在水層時，認(rèn)為稻田深層滲漏由田面水層深度與土壤質(zhì)地決定，其經(jīng)驗公式［20］可表示為：

式中a、b 為取決于稻田土質(zhì)的參數(shù)。

當(dāng)田面不存在水層且稻田土壤含水率大于田間持水率時，稻田深層滲漏依據(jù)式（4）計算［21］：

式中： ms為飽和含水率； md 為殘余含水率； ks為飽和導(dǎo)水率；C 為常數(shù)； RDmx為作物最大根深，mm； RDt為作物當(dāng)日根深，mm； mgt為泥漿層當(dāng)日含水率。

模擬時將泥漿層與儲水區(qū)視為整體進行水量分配。同時依據(jù)氣象、土壤質(zhì)地參數(shù)等數(shù)據(jù)計算逐日傳輸層向飽和層的滲漏量以及儲水區(qū)與泥漿層蒸散發(fā)量?；谥鹑盏乃址至颗c灌溉制度迭代計算地表徑流量與灌溉量，以更新下一日模型水分分量。計算流程見圖7。

圖7 農(nóng)田水轉(zhuǎn)化模擬流程圖

3.2.2 地下水運動過程 各水文響應(yīng)單元逐日地下水位變化量包括該水文響應(yīng)單元垂直方向與地上部分的水量交換和與相鄰水文響應(yīng)單元的水平徑流交換。地下水位的垂直補給與消耗一般取決于河流及渠道補給、深層滲漏等，對于井灌區(qū)地下水開采量為地下水重要的匯項。同時研究區(qū)地下水水力梯度較大，水平流動明顯。根據(jù)裘布依假設(shè)與達西定律，地下水水平徑流量可表示為［22］：

式中：JL為目標(biāo)水文響應(yīng)單元相鄰四個方向任意一個單元向該目標(biāo)單元輸入的逐日單位面積地下水徑流量，m/d；K 為目標(biāo)單元的滲透系數(shù)，m/d；h 為目標(biāo)單元潛水含水層厚度，m；B 為目標(biāo)單元邊長，m；Lx為目標(biāo)單元周圍四個方向上任意單元的地下水位，m；L 為中部目標(biāo)單元的地下水位，m；X 為中部目標(biāo)單元中心點與周圍相鄰單元中心點的距離，m。

3.2.3 灌溉輸配水過程 在ArcGIS 10.7 軟件中將灌區(qū)主要渠道進行數(shù)字化處理，劃分渠道控制區(qū)域（圖8）。根據(jù)水稻灌溉的實際情況，按照灌溉制度確定田間灌溉水量，并基于渠道來水量對各渠道控制的水文響應(yīng)單元進行灌溉水量分配。同時利用經(jīng)驗公式估算渠道輸水滲漏損失，渠道滲漏量公式如下［23］：

圖8 研究區(qū)進水渠道、排水溝道數(shù)字化示意

式中：σ為每千米渠道輸水損失（以渠道凈流量百分?jǐn)?shù)計）；A 和m 分別為渠床土壤透水系數(shù)、渠床土壤透水指數(shù)，二者依據(jù)渠道所在區(qū)域的土壤質(zhì)地確定；Qn為渠道凈流量，m3/s；Q1為渠道滲漏流量，m3/s；L 為渠道長度，km。

依據(jù)公式計算渠道滲漏的同時，根據(jù)渠道灌溉水量與滲漏損失量，實時更新渠道凈流量。

3.2.4 排水過程 在ArcGIS 10.7 軟件中將灌區(qū)主要排水溝道進行數(shù)字化處理，劃分溝道控制范圍（圖8）。研究區(qū)為水稻灌區(qū)，生育期水量充沛，渠道來水很大一部分在補給農(nóng)田灌溉后仍有剩余，故灌區(qū)排水水量由兩部分組成：一部分為田間水層高度超過田埂后產(chǎn)生的田間排水，另一部分為渠道剩余的水量產(chǎn)生的退水。田間排水由各個水文響應(yīng)單元稻田水平衡模型計算而得，退水產(chǎn)生的排水由水文響應(yīng)單元所屬支渠灌溉之后的剩余水量組成。

4 模型率定與驗證

4.1 作物生長過程分別采用2018年與2019年灌區(qū)定位監(jiān)測水稻葉面積指數(shù)（LAI）與株高數(shù)據(jù)進行作物生長模型部分參數(shù)（表1）率定與驗證（圖9）。選取R2與RMSE 作為模擬結(jié)果評價指標(biāo)，結(jié)果表明，EPIC 模型可以較好地模擬研究區(qū)水稻生長過程（表2）。 EPIC 模型基于積溫模擬作物生長，因此區(qū)域上的積溫不確定性會導(dǎo)致定位監(jiān)測數(shù)據(jù)與模擬結(jié)果存在一定偏差。如2019年葉面積指數(shù)模擬最大值出現(xiàn)時間先于實測值，且株高模擬值大多高于實測值應(yīng)是由積溫誤差導(dǎo)致。

圖9 2018、2019年作物生長指標(biāo)模擬結(jié)果與實測值對比

表1 農(nóng)田水轉(zhuǎn)化模型參數(shù)初始值及率定值

表2 農(nóng)田水轉(zhuǎn)化模型率定及驗證效果評價指標(biāo)

4.2 灌區(qū)蒸散發(fā)選取SEBAL 模型的灌區(qū)蒸散發(fā)遙感反演值進行區(qū)域蒸散發(fā)部分參數(shù)的率定驗證，將蒸散發(fā)遙感反演值與模型模擬結(jié)果二者的逐日區(qū)域平均值進行對比，率定蒸散發(fā)模塊參數(shù)（田間持水率、最大作物系數(shù)、蒸騰脅迫曲線影響因子等）。2018、2019年均為豐水年，降雨較多，由于云層覆蓋等原因造成遙感數(shù)據(jù)缺失嚴(yán)重，故選取2018 與2019年兩年生育期內(nèi)影像較為完整時段加以對比（圖10），分別為2018年83 d，2019年65 d。使用2018年數(shù)據(jù)進行模型率定R2=0.78，RMSE=0.98 mm/d；2019年數(shù)據(jù)進行模型驗證R2=0.75，RMSE=0.86 mm/d。認(rèn)為模擬結(jié)果符合精度要求。

圖10 2018、2019年區(qū)域蒸散發(fā)模型模擬值與遙感反演值對比

4.3 地下水埋深采用2018年和2019年地下水埋深數(shù)據(jù)用于地下水運動模塊參數(shù)（導(dǎo)水率、給水度）率定驗證（表3），取2018年地下水埋深數(shù)據(jù)用于模型率定，2019年地下水埋深數(shù)據(jù)用于模型驗證。進行地下水運動模塊模擬時，參考黑龍江水利勘測設(shè)計院編寫的慶安縣1∶5 萬水文地質(zhì)勘察報告（1987年）與黑龍江水文地質(zhì)圖（1∶350 萬）對研究區(qū)域斷裂情況的描述，將研究區(qū)劃分為7 個水文地質(zhì)分區(qū)進行模擬（圖11）。模擬結(jié)果如圖12，可見各點均勻分布在1∶1 相關(guān)線兩側(cè)附近，可認(rèn)為模型無系統(tǒng)性錯誤。由表4 可見各種模型評價指標(biāo)均在允許范圍內(nèi)，地下水埋深模擬值與實測值吻合較好，模型運行穩(wěn)定、可靠，可用于灌區(qū)供耗排耦合過程定量模擬。2018、2019年模擬結(jié)果評價指標(biāo)在灌區(qū)范圍內(nèi)呈現(xiàn)出下游優(yōu)于上游的情況，其原因在于研究區(qū)由上游向下游地下水埋深遞增，埋深較淺的地下水位觀測井?dāng)?shù)據(jù)更易受到偶然因素干擾，但埋深較深的地下水位觀測井實際情況與模型假設(shè)較為一致，因此模擬結(jié)果較為理想。2018、2019 兩年由于灌區(qū)生育期內(nèi)持續(xù)的田間深層滲漏以及渠道入滲補給，灌區(qū)整體地下水位呈上升趨勢，但井灌區(qū)由于存在抽取地下水用于灌溉的原因，其地下水位升高幅度低于渠灌區(qū)。

表4 地下水運動模塊率定及驗證效果評價指標(biāo)

圖12 率定期、驗證期地下水埋深計算值與觀測對比

表3 潛水給水度Sy 、滲透系數(shù)ks2 初始值與率定值

圖11 地質(zhì)分區(qū)示意

4.4 排水流量2018年對一支排及七支排溝排水流量數(shù)據(jù)進行監(jiān)測，對兩條排水溝排水流量模擬結(jié)果與實測結(jié)果進行了對比（圖13），其中一支排溝R2=0.67，RMSE=0.38 m3/s，七支排溝R2=0.61，RMSE=0.15 m3/s。2019年由于儀器故障，缺測該年排水溝道流量數(shù)據(jù)。結(jié)果表明，模型能較好地模擬灌區(qū)排水過程，但峰值時誤差較大，分析認(rèn)為誤差來源是將排水溝道控制范圍內(nèi)的田間排水與渠道退水等統(tǒng)一劃歸到該排水溝，未考慮溝渠聯(lián)通等情況。

圖13 2018年一支排、七支排流量模擬結(jié)果與實測值對比

5 討論

通過作物生長指標(biāo)、地下水位、蒸散發(fā)及排水流量等驗證表明，本文所建立的水稻灌區(qū)基于多過程耦合的分布式水轉(zhuǎn)化模型具有較高精度，可以客觀反映水稻灌區(qū)復(fù)雜的水轉(zhuǎn)化過程。相對于目前常用的灌區(qū)水轉(zhuǎn)化模型來說，該模型針對水稻灌區(qū)灌排復(fù)雜及農(nóng)田灌溉水與地下水轉(zhuǎn)化頻繁的特點，實現(xiàn)了灌溉過程、排水過程、農(nóng)田水分運動及作物生長過程、地下水運動過程的動態(tài)耦合。盡管眾多學(xué)者嘗試將流域水文模型（如SWAT 模型）與地下水運動模型（如MODFLOW 模型）進行耦合用以模擬灌區(qū)水轉(zhuǎn)化過程，但由于流域水文模型對于土壤垂向水分運動的概化，難以真實反應(yīng)農(nóng)田灌溉水入滲與蒸發(fā)過程［7］。特別對于稻田具有明顯水分分層特性條件下，本研究中所提出的稻田水轉(zhuǎn)化過程具有明顯的優(yōu)勢。進一步來說，稻田分層水轉(zhuǎn)化模型較傳統(tǒng)基于動力學(xué)過程對土壤水分運動模型參數(shù)依賴性減小，對灌區(qū)分布式水轉(zhuǎn)化模型具有較強的適用性。

此外，該模型將灌排渠溝進行數(shù)字化，突出了灌區(qū)灌排過程的時空差異性。尤其是對渠系輸配水過程中滲漏損失及對地下水的補給這一灌區(qū)水轉(zhuǎn)化過程特有的環(huán)節(jié)得以量化表征，一定程度提高了模型精度。相比較而言，目前基于農(nóng)田水文模型的分布式模擬及灌區(qū)水平衡分布式模擬模型均對復(fù)雜渠溝空間分布及其輸配水及排水過程缺乏必要的考慮［6，12-13］。從模型結(jié)構(gòu)來看，該模型由灌溉模塊、排水模塊、農(nóng)田作物生長與水轉(zhuǎn)化模塊以及地下水模塊耦合而成，模型中農(nóng)田水轉(zhuǎn)化過程、地下水運動過程參數(shù)較現(xiàn)有基于動力學(xué)過程的模擬方法大大減少，特別是模型所采取的土壤水與地下水過程統(tǒng)一的網(wǎng)格化水文響應(yīng)單元提高了灌區(qū)尺度分布式模型的模擬效率。

需要指出的是，本文所建立的水稻灌區(qū)分布式水轉(zhuǎn)化模型依賴于灌區(qū)渠溝分布、高分辨率土地利用及作物種植分布等空間數(shù)據(jù)。同時，相對于灌區(qū)水轉(zhuǎn)化客觀物理過程，模型對于排水過程模擬尚未考慮通過田埂側(cè)向滲流與溝渠聯(lián)通情況，一定程度引起了排水模擬的誤差；模型中灌溉模塊所涉及的灌溉制度需進一步基于渠道來水及作物需水進行優(yōu)化設(shè)計。

6 結(jié)論

本文針對水稻灌區(qū)水轉(zhuǎn)化特點，以黑龍江省和平灌區(qū)為例，建立了灌區(qū)基于多過程耦合的分布式水轉(zhuǎn)化模型?；诠鄥^(qū)實測水文過程對模型加以檢驗，結(jié)果表明模型可以較好模擬灌區(qū)主要水轉(zhuǎn)化過程。與傳統(tǒng)灌區(qū)水轉(zhuǎn)化模型相比，該模型實現(xiàn)了灌溉過程、排水過程、農(nóng)田水分運動及作物生長過程、地下水運動過程的動態(tài)耦合，可以有效表征灌區(qū)多尺度復(fù)雜水轉(zhuǎn)化過程。特別在農(nóng)田尺度，通過將稻田水轉(zhuǎn)化模型與作物生長模型耦合，克服了傳統(tǒng)灌區(qū)水轉(zhuǎn)化模型對作物生長過程模擬不足的劣勢。模型中將灌排渠溝進行數(shù)字化，突出了灌區(qū)灌排過程的時空差異性，實現(xiàn)了灌區(qū)水轉(zhuǎn)化的分布式模擬。該模型對于灌區(qū)用水效率時空動態(tài)評估、灌溉水資源管理等具有重要意義。后續(xù)研究中將對模型灌溉模塊和排水模塊優(yōu)化提升，進一步提升模型的合理性及普適性。