水稻灌區基于多過程耦合的分布式水轉化模型

王 璞，盧 炤，霍再林

（中國農業大學 水利與土木工程學院，北京 100083）

1 研究背景

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

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

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

2 研究區域及水文監測

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

圖1 研究區位置

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

圖2 研究區DEM 數據及水尺、觀測井布設位置示意

圖3 研究區土壤參數空間分布

3 灌區分布式水轉化模型

3.1 模型框架

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

圖4 模型整體結構圖

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

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

3.1.3 水文響應單元劃分 基于ArcGIS 10.7 軟件對研究區進行水文響應單元劃分，水文響應單元大小為250 m×250 m，灌區內有效水文響應單元共2715 個。

3.2 各過程定量表征方法

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

圖6 土壤垂向分層示意

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

稻田存在水層：

稻田不存在水層：

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

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

當田面存在水層時，認為稻田深層滲漏由田面水層深度與土壤質地決定，其經驗公式［20］可表示為：

式中a、b 為取決于稻田土質的參數。

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

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

模擬時將泥漿層與儲水區視為整體進行水量分配。同時依據氣象、土壤質地參數等數據計算逐日傳輸層向飽和層的滲漏量以及儲水區與泥漿層蒸散發量。基于逐日的水分分量與灌溉制度迭代計算地表徑流量與灌溉量，以更新下一日模型水分分量。計算流程見圖7。

圖7 農田水轉化模擬流程圖

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

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

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

圖8 研究區進水渠道、排水溝道數字化示意

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

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

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

4 模型率定與驗證

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

圖9 2018、2019年作物生長指標模擬結果與實測值對比

表1 農田水轉化模型參數初始值及率定值

表2 農田水轉化模型率定及驗證效果評價指標

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

圖10 2018、2019年區域蒸散發模型模擬值與遙感反演值對比

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

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

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

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

圖11 地質分區示意

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

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

5 討論

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

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

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

6 結論

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