基于DRAINMOD模型的不同灌排模式稻田水氮運移模擬

謝陽村，徐 敏*，高世凱，沈連起

（1.生態(tài)環(huán)境部 環(huán)境規(guī)劃院，北京 100000；2.華北水利水電大學(xué) 水利學(xué)院， 鄭州 450046；3.壽光市水利局，山東 壽光 262700）

0 引 言

【研究意義】氮素是農(nóng)作物生長發(fā)育的關(guān)鍵因素，在傳統(tǒng)灌排模式下, 由于灌溉和降雨造成的排水較多，大量氮素直接進入水體，導(dǎo)致農(nóng)業(yè)面源污染[1-4]。已有研究表明，單一水稻控制灌溉與控制排水技術(shù)可高效利用養(yǎng)分和水分，減輕環(huán)境負(fù)荷，實現(xiàn)節(jié)水高產(chǎn)、減排、控污的目標(biāo)[5-6]。但將水稻控制灌溉與控制排水技術(shù)進行有機結(jié)合，形成水稻全生育期的控制灌排技術(shù)[4,7]，進行節(jié)水減排控污效果的研究卻較少。此外，農(nóng)業(yè)面源污染區(qū)域發(fā)生具有隨意性、污染排放不確定性以及空間異質(zhì)性等特點，對其影響過程、影響因素以及產(chǎn)生的機理研究造成難度[8-9]?！狙芯窟M展】輔助于完善農(nóng)田排水管理以及氮素流失的農(nóng)田排水模型已經(jīng)引起國內(nèi)外的廣泛關(guān)注[10-12]，但對于稻田的排水和氮素運移模型較為少見。DRAINMOD 模型是美國農(nóng)業(yè)部推薦的農(nóng)田排水模型，模型不斷發(fā)展更新，并逐漸得到推廣應(yīng)用[13-15]。部分學(xué)者利用DRAINMOD模型嘗試應(yīng)用于稻田，認(rèn)為該模型能夠較好地模擬灌區(qū)稻田徑流排水過程和地下水位動態(tài)變化[16-18]，也有學(xué)者利用DRAINMOD 模型對暴雨排水過程氮素變化進行探究[19-20]。高學(xué)睿等[21]將DRAINMOD 模型應(yīng)用到氣候相對濕潤的湖北漳河灌區(qū)，對灌區(qū)稻田的排水和氮素流失規(guī)律進行模擬，發(fā)現(xiàn)模型可有效地模擬預(yù)測南方典型灌區(qū)稻田排水特性及氮素運移特性，相對誤差范圍為2.8%～16.0%。Hashemi 等[22]發(fā)現(xiàn)DRAINMOD 模型可以較好地模擬不同排水深度和間距組合下稻田硝態(tài)氮損失和排水流量，相對誤差范圍為1.0%～35.0%?！厩腥朦c】綜上，對于DRAINMOD模型在稻田的應(yīng)用研究多基于傳統(tǒng)灌排模式，但對于不同灌排模式下稻田的應(yīng)用研究較少。【擬解決的關(guān)鍵問題】采用Morris 方法對DRAINMOD 模型的水氮運移參數(shù)進行靈敏度分析，利用DRAINMOD 模擬不同灌排模式稻田水氮動態(tài)變化規(guī)律，探討DRAINMOD模型在水稻不同灌排模式下的適用性，以期為我國南方稻作區(qū)灌排實踐提供科學(xué)依據(jù)。

1 材料與方法

1.1 試驗概況及農(nóng)藝措施

試驗研究于2017年和2018年6―10月在漣水縣水利科學(xué)研究站的試驗大田完成，試驗區(qū)位于江蘇省淮安市漣水縣朱碼鎮(zhèn)王二莊村境內(nèi)，地理坐標(biāo)東經(jīng)119°16′，北緯33°50′，屬于亞熱帶濕潤性氣候區(qū)，多年平均氣溫14.4 ℃，年降水979.1 mm，降水量時程分布不均，汛期降水量占全年降水量比重較大，年均蒸發(fā)1 385.4 mm，年均日照時間2 280 h，年均無霜期240 d。試驗區(qū)0～100 cm 土層土壤為壤土，土壤體積質(zhì)量為1.42 g/cm3，pH 值為6.82，田間持水率為27.9%（質(zhì)量），有機質(zhì)量為2.19%，全氮量、全磷量分別為0.98、1.12 g/kg，速效氮量為21.5 mg/kg。采用格田小區(qū)試驗，格田規(guī)格為90 m×27 m，格田田埂高30 cm以上，底部寬40 cm，上部寬30 cm，除農(nóng)溝一側(cè)的田埂外，格田周邊嵌入50 cm 的薄膜并覆蓋至田梗，避免了各小區(qū)間水位影響。試驗區(qū)布置見圖1。

圖1 試驗區(qū)布置 Fig.1 The layout of the test area

供試水稻品種為兩優(yōu)9918，該品種株型松緊適中，群體整齊，莖稈堅韌，葉色淡綠，葉片內(nèi)卷挺直，熟期落色好，抗病蟲害性能較好，耐高溫能力強。該品種有效穗數(shù)241.5 萬穗/hm2，穗總粒數(shù)175.7粒，結(jié)實率85.5%，千粒質(zhì)量27.5 g，株高可達122.7 cm 左右，全生育期平均為130.8 d，多年平均產(chǎn)量為597.1 kg。水稻農(nóng)藝措施見表1，復(fù)合肥N、P、K 質(zhì)量比為15∶15∶15，尿素含氮量為46.4%。

表1 2017―2018年水稻農(nóng)藝措施 Table 1 Agronomic measures of rice in 2017―2018

1.2 試驗設(shè)計

考慮到南方地區(qū)水稻生長期與汛期同季，控制灌排模式可以在保持較低灌水下限的同時，提高雨后蓄水上限，擴大稻田的儲水庫容，達到充分利用天然降水，減少灌排水頻率、灌排定額和氮磷排放量，并依據(jù)前人連續(xù)幾年的試驗結(jié)果[7,23,24]及江蘇省水稻節(jié)水灌溉技術(shù)規(guī)范，各處理水位調(diào)控方案見表2。

表2 各處理農(nóng)田水位調(diào)控方案 Table 2 Water control program of each treatment

選取水稻分蘗期、拔節(jié)孕穗期、抽穗開花期和乳熟期4 個生育階段進行控水試驗，共3 個處理，每個處理設(shè)3 個重復(fù)。傳統(tǒng)灌排（CK）農(nóng)田水位達到灌水下限值-20 cm 立即灌水至6 cm，雨后水層深度超過6 cm 即排水?？刂乒嗯欧譃檩p旱控排（LCID）和重旱控排（HCID），LCID 達到灌水下限值-20 cm 立即灌水至6 cm，分蘗期雨后水層深度超過10 cm 即排水，其他生育期雨后水層深度超過20 cm 即排水；HCID 達到灌水下限值-50 cm 立即灌水至6 cm，雨后排水方式與LCID 相同。

1.3 試驗觀測指標(biāo)及方法

1）氣象資料。由田間試驗小區(qū)附近的自動氣象站收集，主要包括降雨時間、降雨歷時及降雨量和大氣溫度、大氣壓強、大氣相對濕度、風(fēng)速、水面蒸發(fā)量、有效輻射量等。

2）農(nóng)田水位。每天09:00 對農(nóng)田水位進行觀測。當(dāng)農(nóng)田地面有水層時，通過豎尺在固定觀測點測量農(nóng)田地面水層深度。無水層時，通過在試驗田塊中間等距離安裝的3 個地下水位觀測井記錄各小區(qū)農(nóng)田水位的埋深。

3）灌水量。通過在農(nóng)渠上安裝的三角量水堰測定各格田灌溉水量。

4）排水量。采用水位差法，記錄排水前后格田水層深度，根據(jù)二者差值計算排泄水量，其中當(dāng)降雨引起的徑流時應(yīng)加上排水時段內(nèi)的降雨量，日內(nèi)發(fā)生多次排水按1 次計算。

5）水樣提取及分析。在格田小區(qū)選擇4 個取樣點，地表水用塑料注射器不擾動土層分別在水面、1/2水深、底部取水，混合裝入采樣瓶，一般5～7 d 取1 次水樣，但施肥后以及暴雨后1 周內(nèi)隔日取樣；滲漏水通過地下水觀測井，深度為1.5 m，井管底部密封，管壁打孔并用無紡布過濾，井口加蓋防雨。取樣時先用真空泵清空井內(nèi)積水，1～2 h 后吸取井管中滲漏水，取樣頻率同地表水。排水時在排水出口取水，每30 min 取1 次，當(dāng)產(chǎn)流均勻時可適當(dāng)延長步長。所有水樣采樣后貼好標(biāo)簽，并進行冷藏（3 ℃）處理，在24 h 內(nèi)進行水質(zhì)分析。水樣中銨態(tài)氮、硝態(tài)氮分別采用絮凝沉淀納氏試劑光度法、紫外分光光度法，所用儀器主要為UV2800 島津紫外分光光度儀。

2 模型原理及主要參數(shù)

2.1 模型原理

DRAINMOD 模型以簡單的土壤剖面范圍內(nèi)水平衡原理為基礎(chǔ)來描述農(nóng)田水文變化過程，可以準(zhǔn)確地預(yù)測出農(nóng)田水位、地表徑流和地下排水量以及作物產(chǎn)量等。在時間增量Δt內(nèi)，土壤水量平衡方程可表示為：

式中：ΔV為土體水量的變化量（cm）；D為側(cè)向排水量（cm）；ET為蒸發(fā)蒸騰量（cm）；DS為深層滲漏量（cm）。

在單元Δt時段內(nèi)，地表水量平衡方程為：

式中：ΔS為地表蓄水量的變化量（cm）；P為降雨量（含灌水量）（cm）；F為入滲量（cm）；RO為地表徑流量（cm）。

DRAINMOD 模型可以用于模擬濕潤地區(qū)淺水位條件下農(nóng)田非飽和區(qū)（一維，垂向）及飽和區(qū)（二維，垂向和側(cè)向）的土壤水和氮素轉(zhuǎn)化運移，模型采用多相一維對流彌散方程表達式來模擬計算農(nóng)田中的氮素運移，并用一階有限差分法求解。多相一維對流彌散方程式為：

式中：θa為土壤液相體積分?jǐn)?shù)；θg為土壤氣相體積分?jǐn)?shù)；ρb為土壤固相干體積質(zhì)量（M/L3）；Ca為某種形式氮素液相濃度（M/L3）；Cg為某種氮素形式氣相濃度（M/L3）；Cs為某種形式氮素固相濃度（M/M）；Da為水動力擴散系數(shù)（L2/T）；dg為分子擴散系數(shù)；Va液相體積通量；S為源匯項；t為時間（T）；z為空間坐標(biāo)（L）。

表3 土壤參數(shù)和排水參數(shù) Table 3 Main input parameters of soil and drainage

表4 氮素參數(shù) Table 4 Nitrogen parameters

2.2 模型參數(shù)

土壤參數(shù)通過入滲試驗、野外鉆孔法、環(huán)刀法等方法測定。由于原系統(tǒng)針對的是暗管排水，在應(yīng)用于明溝排水時需要進行相應(yīng)的參數(shù)調(diào)整，排水深度和排水間距試驗測定；排水管有效半徑和地面到不透水層深度通過模型率定；由于土壤質(zhì)地和試驗條件相似，排水系數(shù)和Kirkham積水深度參考文獻[25]確定。氮素相關(guān)參數(shù)通過模型率定和參考文獻[26-27]確定。作物參數(shù)有效根系深度根據(jù)2017年和2018年實測數(shù)據(jù)以及相關(guān)經(jīng)驗進行設(shè)定。模型初始參數(shù)見表3—表5。采用2017年試驗觀測數(shù)據(jù)進行參數(shù)率定、2018年數(shù)據(jù)對模型進行驗證。

表5 作物參數(shù) Table 5 Growth parameters

3 模型參數(shù)敏感性分析

3.1 模型參數(shù)選擇及分布

敏感性分析分為局部敏感性及全局敏感性，其中全局敏感性考慮了參數(shù)之間的相互影響，更加適用于DRAINMOD 模型。將參數(shù)設(shè)定為均勻分布，其參數(shù)分布范圍參考文獻[25]，具體見表6。

表6 模型參數(shù)服從的分布 Table 6 Obedient probability distribution of model parameters

水分運移模擬中水位變化選擇通過衡量模擬值與觀測值的擬合程度進而評估參數(shù)敏感性，即Nash Sutcliffe 效率（Nash）、相對誤差總量RE以及相關(guān)系數(shù)CC，計算式為：

式中：Qi為農(nóng)田水位實測值（cm）；Mi為農(nóng)田水位模擬量（cm）；n為模擬次數(shù)；為農(nóng)田水位實測平均值（cm）；為農(nóng)田水位模擬平均值（cm）。

其他參數(shù)采用定性分析的Morris 法[28]進行敏感性分析，該方法通過參數(shù)變化效應(yīng)的平均值μ和標(biāo)準(zhǔn)差σ來判斷參數(shù)的敏感性，其中μ值反映參數(shù)敏感性大小，μ值越大，代表參數(shù)對模型輸出結(jié)果的影響越大；σ值反映參數(shù)間交互作用的大小或者表示參數(shù)為非線性影響，σ值越大，代表參數(shù)交互作用越大。

3.2 水分運移模擬敏感性分析

對于選取的水分運移參數(shù)，依據(jù)Morris 方法共進行20 次抽樣，得到100 組模型參數(shù)的抽樣樣本，運行模型得出各個參數(shù)的靈敏度，結(jié)果如圖2 所示。

由目標(biāo)函數(shù)均值μ可以得出，I.L對模型農(nóng)田水位輸出結(jié)果影響最小，De影響最大，三層側(cè)向飽和導(dǎo)水率在不同目標(biāo)函數(shù)上有所區(qū)別，但Ksat2影響最大，尤其在Nash和CC上響應(yīng)較明顯；對于地表排水量和地下滲漏量，Ksat2與De的影響最大，I.L影響最小。

圖2 水文參數(shù)均值與標(biāo)準(zhǔn)差散點 Fig.2 Mean value and standard value of hydrological parameters

由目標(biāo)函數(shù)標(biāo)準(zhǔn)差σ可以得出，整體呈現(xiàn)出“山谷”型，De受其他參數(shù)的影響（或非線性作用）最為明顯，而Ksat2基本不受其他參數(shù)的影響；對于地表排水量和地下滲漏量，I.L處于較高水平，De和Ksat2受其他因素影響較小。

3.3 氮素運移模擬敏感性分析

對于選取的氮素運移參數(shù)，依據(jù)Morris方法共進行20次抽樣，采樣220組數(shù)據(jù)得到220組模型參數(shù)的抽樣樣本，運行模型得出各個參數(shù)的靈敏度，如圖3所示。

由目標(biāo)函數(shù)均值μ可以得出，對于地表銨態(tài)氮徑流損失，λ對輸出結(jié)果的影響最明顯，因為擴散是銨態(tài)氮運移的主要機理，當(dāng)銨態(tài)氮垂向運移距離較短時會導(dǎo)致較多的銨態(tài)氮通過徑流方式排出，其次為Toptdec，因為當(dāng)土壤溫度與分解溫度相近時，有利于有機質(zhì)分解為銨態(tài)氮；對于地表硝態(tài)氮徑流損失，反硝化反應(yīng)起主導(dǎo)作用，反硝化參數(shù)Vmaxden、Kmden、Toptden對輸出結(jié)果影響最大，有機質(zhì)分解參數(shù)Toptdec、Kdeca、Kdecs影響最??；對于地下銨態(tài)氮滲漏損失，輸出結(jié)果基本只受到有機質(zhì)分解和硝化反應(yīng)參數(shù)影響，依次為Kmnit、Toptdec、Vmaxnit、Kdeca、Topt-nit、Kdecs；對于地下硝態(tài)氮滲漏損失，輸出結(jié)果基本只受到硝化反應(yīng)及反硝化反應(yīng)參數(shù)影響，依次為Vmaxden＞Kmnit＞Vmaxnit＞Kmden＞Toptden＞Toptnit。

由目標(biāo)函數(shù)標(biāo)準(zhǔn)差σ可以得出，整體呈“山峰”型趨勢，地表銨態(tài)氮徑流損失參數(shù)Toptdec、Toptnit、Kdeca、λ相互影響較大；地表硝態(tài)氮徑流損失參數(shù)中僅Kmden略受其他參數(shù)影響；地下銨態(tài)氮滲漏損失Toptnit和Kdeca相互影響較大；地下硝態(tài)氮滲漏損失Kmnit受其他參數(shù)影響大。

圖3 氮素參數(shù)均值與標(biāo)準(zhǔn)差散點 Fig.3 Mean and standard values of nitrogen

圖4 農(nóng)田水位變化模擬值和實測值對比 Fig.4 Simulated and measured values of farmland water level

表7 農(nóng)田水位擬合性能 Table 7 Fitting performance of farmland water level

4 模擬結(jié)果與分析

4.1 不同灌排模式下稻田水分動態(tài)

各處理農(nóng)田水位變化模擬值和實測值對比如圖4所示。從圖4 可以看出，實測的農(nóng)田水位與模型模擬結(jié)果之間吻合較好，模擬結(jié)果能較好地反映田間農(nóng)田水位的變化情況。其中，正水位模擬程度較好，但負(fù)水位（即水位低于土壤表面）則模擬結(jié)果存在較大差別，但屬于可接受范圍。由表7 可以得出，農(nóng)田水位擬合性能參數(shù)Nash在0.81～0.97 之間，均大于0.8，模擬效果較好；Re絕對值在0.01～0.53 之間，且主要為地下水位模擬誤差；CC在0.86～0.98 之間，模擬值與實測值變化趨勢較為一致。

不同灌排模式灌排水量模擬值和實測值見表7。各處理灌水量和排水量模擬值與實際值差別小于10%，模擬效果較好。與CK 相比，2017年LCID 和HCID 灌水量分別減少16.1%、37.1%，排水量分別減少68.6%、72.6%；2018年LCID 和HCID 灌水量分別減少9.7%、26.2%，排水量分別減少33.0%、41.9%。說明控制灌排模式通過合理調(diào)蓄雨水資源，減少稻田排水次數(shù)和排水量，對于節(jié)約灌溉用水有重要意義。而HCID 與LCID 相比，灌水量更少，說明節(jié)水省工效果更加顯著。

4.2 不同灌排模式下稻田氮素動態(tài)

各處理氮素徑流損失累計量如圖5 所示。從整體趨勢看，各處理模擬值與實測值在前期徑流損失量較大，在移栽后幾天發(fā)生較大變化是由于生育前期連續(xù)遇雨，而返青期和分蘗初期需要保持淺水層，因此地表排水量較大，且此時基肥的施入引起地表水中氮素濃度較高。其中，CK 徑流損失量最大，主要是由于前期排水量較多，而控制灌排攔蓄雨水，減少了氮素的排放。與實測值相比，銨態(tài)氮徑流損失模擬較好，雖然在2018年存在些許差異，主要是由于排水影響，而非氮素濃度導(dǎo)致，硝態(tài)氮徑流損失模擬存在較大差別。

圖5 氮素徑流損失累計量 Fig.5 Cumulative nitrogen runoff loss

氮素徑流損失累計總量模擬值與實測值對比見表8。各處理氮素徑流損失累計總量模擬值與實際值差別小于20%，模擬效果較好。從銨態(tài)氮徑流損失累計總量看，銨態(tài)氮模擬值與實際值相差不大，且與地表排水量模擬值變化情況保持一致；硝態(tài)氮徑流損失累計總量在2017年各處理均有一定程度的升高，而2018年CK 模式下有所降低，這主要是模型對排水中硝態(tài)氮質(zhì)量濃度模擬存在偏差所致。與CK 相比，2017年LCID 和HCID 銨態(tài)氮徑流損失累計總量分別減少41.2%、45.0%，硝態(tài)氮徑流損失累計總量分別減少49.6%、53.1%；2018年LCID 和HCID 銨態(tài)氮徑流損失累計總量分別減少43.6%、47.4%，硝態(tài)氮徑流損失累計總量分別減少29.8%、32.8%。說明控制灌排模式約束了地表水高質(zhì)量濃度氮素時期的稻田排水，實現(xiàn)了減排控污的目標(biāo)。

表8 不同處理氮素徑流損失累計總量 Table 8 Fitting performance of nitrogen runoff loss

5 討 論

模型參數(shù)的合理選取對于模擬結(jié)果的精確性起到關(guān)鍵作用，而大多數(shù)參數(shù)是依靠試驗數(shù)據(jù)率定間接取得的，采用靈敏度分析來評價各個參數(shù)的不確定性對模擬結(jié)果的影響，確定出比較敏感的參數(shù)，有針對性地對這些參數(shù)進行測量和調(diào)節(jié)，從而選取合理的模型參數(shù)[29]。通過DRAINMOD 模型進行模擬時，20～40 cm 土層側(cè)向飽和導(dǎo)水率Ksat2的變化對稻田水分運移結(jié)果影響較大，但與其他參數(shù)相互作用不明顯，進行模擬輸入時，可配合其他參數(shù)在適當(dāng)范圍內(nèi)進行微調(diào)，來保證模擬精度。不透水層深度I.L和排水管有效管徑De極易受其他參數(shù)影響，故其初始輸入值的精度直接影響模型的模擬效果。銨態(tài)氮模擬的準(zhǔn)確性取決于彌散系數(shù)、硝化反應(yīng)參數(shù)、有機質(zhì)適宜分解溫度，硝態(tài)氮模擬的準(zhǔn)確性取決于反硝化參數(shù)和硝化參數(shù)[19,21]，應(yīng)該加強試驗手段保證現(xiàn)場測量的精度。

DRAINMOD 模型進行模擬計算時，當(dāng)農(nóng)田水位為負(fù)值時，模擬值與實測值之間存在較大偏差，這可能是因為地下水位較淺時，土壤水力學(xué)特性對地下水位的預(yù)測比較敏感[14]，因此模擬效果存在偏差，但并未顯著影響到灌水及灌水后地表水層深度。與采用 HYDRUS-1D 模型模擬稻田中水分運移的研究相比較[24]，DRAINMOD 模型表現(xiàn)良好。銨態(tài)氮模擬值與實際值相差不大，且與地表排水量模擬值變化情況保持一致，說明銨態(tài)氮徑流負(fù)荷損失模擬偏差主要是由排水量所引起的。而DRAINMOD 模型低估了生育前期排水中硝態(tài)氮濃度，高估了中后期硝態(tài)氮濃度，使得曲線呈現(xiàn)差別，還需進一步研究。

本研究對不同灌排模式下灌排水量和氮素徑流負(fù)荷進行模擬，結(jié)果發(fā)現(xiàn)控制灌排模式將水稻控制灌溉技術(shù)與控制排水相結(jié)合，較傳統(tǒng)灌排模式，可高效利用養(yǎng)分和水分，充分發(fā)揮稻田的濕地效應(yīng)，減少灌排水量和稻田氮磷污染物負(fù)荷，實現(xiàn)節(jié)水減排的目標(biāo)，與前人研究結(jié)果[28]一致。同時，已有研究發(fā)現(xiàn)控制灌排通過調(diào)控稻田水分狀況，可以實現(xiàn)水稻高產(chǎn)，與傳統(tǒng)灌排相比，減產(chǎn)少于10%，且差異不顯著[7,30-31]。但是，對水稻田來說，由于田間環(huán)境有很大差異，稻作區(qū)的溝、田水體具有相互影響、協(xié)同控制的特點，DRAINMOD 模型對于不同灌排模式的適用性需要進一步在不同田間尺度、不同水文年型和不同地理分區(qū)中進行檢驗與推廣。

6 結(jié) 論

1）20～40 cm 土層側(cè)向飽和導(dǎo)水率Ksat2對稻田水分運移模擬結(jié)果影響最為明顯，同時各個參數(shù)間存在非線性作用，以不透水層深度I.L和排水管有效管徑De受其他參數(shù)影響最為明顯；彌散系數(shù)λ、硝化反應(yīng)參數(shù)（Kmnit、Vmaxnit、Toptnit）、有機質(zhì)適宜分解溫度Toptdec對稻田銨態(tài)氮運移模擬結(jié)果影響較大；反硝化參數(shù)（Vmaxden、Kmden、Toptden）和硝化參數(shù)（Kmnit、Vmaxnit、Toptnit）對稻田硝態(tài)氮運移模擬結(jié)果影響較大。

2）DRAINMOD 模型能夠較好模擬不同灌排模式的稻田水位，擬合度在正水位時較高，但在負(fù)水位時較低；DRAINMOD 模型能夠較好地模擬不同灌排模式氮素運移，其中銨態(tài)氮的擬合度較好。控制灌排模式通過合理調(diào)蓄雨水資源，約束地表水高濃度氮素時期的稻田排水，可以實現(xiàn)節(jié)水減排的目的。