寧夏引黃灌區種植業面源污染流失量模擬

佘冬立，胡磊，夏永秋，劉文娟，李虹，馬建軍，馬琨

（1.河海大學農業科學與工程學院，南京 211100；2.中國科學院南京土壤研究所，南京 210008；3.寧夏大學農學院，銀川750021；4.寧夏農業環境保護監測站，銀川 750002；5.寧夏大學生態環境學院，銀川 750021；6.寧夏大學西北土地退化與生態恢復國家重點實驗室培育基地/西北退化生態系統恢復與重建教育部重點實驗室，銀川 750021）

寧夏引黃灌區位于我國西北內陸干旱半干旱地區，是我國四大古灌區之一，糧食產量和農業產值占自治區的三分之二以上，在區域農業經濟發展中發揮著重要作用。灌區農業生產集約化程度高、管理方式粗放，大量氮磷污染物質隨退水經各級排水溝排入黃河[1]，對黃河上游水質具有直接影響[2]。近年來，隨著國家對農田化肥和規?；B殖投入的增加以及農村生活污水盲目排放規模的擴大，寧夏引黃灌區面源污染呈現加劇的趨勢[3]，以農田為核心的灌區種植業退水污染是主要污染負荷來源[4]，春灌和冬灌是灌區種植業面源污染流失高發期。大量施肥后，農田土壤通過退水向周圍水體遷移排放氮磷等污染物質，經過沿程農溝、斗溝、支溝消納后，最終進入干溝和黃河，在造成土壤養分流失的同時，也加劇了地表水和地下水水質惡化的風險。

農業面源污染過程難以直接監測，模型成為解析農業面源污染負荷的主要工具。比較著名的過程模型是由美國農業部農業研究所開發的AGNPS[5]和美國國家環保局開發的SWAT和BASINS等模型[6]，經驗模型有Johnes[7]輸出系數模型。目前，引黃灌區農業面源污染模型研究較少。胡亞偉等[8]以青銅峽灌區第一排水溝為典型小流域，結合3S 技術，根據灌區水系及灌溉特點對SWAT 模型進行改進，建立了寧夏引黃灌區面源污染負荷模型，但上述過程模型本質上是建立在西方以旱作為主的農業流域基礎之上，其過程模擬基于超滲產流理論，不能反映帶有田埂水田退水的實際情況，并且模擬所需參數眾多，獲取難度大。李強坤等[9]基于單元分析的觀點，提出了負荷貢獻率的概念，建立了寧夏引黃灌區農業非點源改進輸出系數估算模型，但此類經驗模型氮磷輸出系數的確定依賴于點尺度實驗或者樣地調查，很難體現不同區域、不同時段輸出系數的差異，并且不能反映農業管理措施與面源負荷的響應關系，應用性與擴展性較差。因此，面對引黃灌區水力聯系的復雜性和不確定性，為了降低模型應用難度、提高農業管理效率，亟需構建適合寧夏引黃灌區面源污染流失特點的輕簡化模型。

本研究基于寧夏引黃灌區典型種植業面源污染流失過程的理論依據和觀測數據，構建了水田和旱地氮磷流失量模型，揭示了灌區種植業氮磷流失時空特征與空間格局，為改善灌區溝道及黃河水質環境提供科學有效的措施，對實現寧夏回族自治區社會經濟可持續發展具有科學指導意義。

1 材料與方法

1.1 研究區概況

寧夏引黃灌區（37°49′~39°23′ N，105°37′~106°39′E）位于寧夏回族自治區中北部，東至鄂爾多斯臺地西緣，西至賀蘭山東，南至寧夏青銅峽水利樞紐，北至石嘴山，屬于溫帶大陸性干旱氣候，全年干旱少雨。灌區含括了銀川市、石嘴山市、吳忠市在內的3 個地級市及永寧、惠農、青銅峽等11 個區縣（圖1）。灌溉面積為 33 萬 hm2[11]，占總土地面積（6 239 km2）[10]的53%，2021 年引黃水量為 37.58 億 m3，排水量為 16.17億m3，多年平均引黃水量為40.95 億m3，排水量為22.09 億m3。灌區北高南低，地勢較為平坦，坡降約1/4 000[12]。

圖1 寧夏引黃灌區位置及土地利用[13]Figure 1 The location and land use of the Yellow River irrigation area of Ningxia[13]

面源污染過程觀測主要在引黃灌區境內寧夏大學實驗農場開展。實驗農場位于寧夏中部黃河沖積平原永寧縣境內，該地年降水量180~200 mm，年蒸發量高達2 800~3 200 mm，年均氣溫8.5 ℃，年均≥10 ℃積溫3 300 ℃，無霜期140~160 d，平均日照時數3 000 h，日溫差13 ℃。觀測區包含兩塊地塊及兩條主要溝道，西邊玉米地塊面積23.35 hm2，東邊水稻田地塊4.53 hm2。玉米地合計施用氮肥（以N計）540 kg·hm-2，磷肥 400 kg·hm-2；水稻田合計施用氮肥 355 kg·hm-2，磷肥 270 kg·hm-2；小麥地合計施用氮肥 294 kg·hm-2，磷肥162 kg·hm-2[14]。稻田農溝長254 m，分布有3個水質監測點，兩塊玉米地之間的支溝長500 m，分布有4個水質監測點，首尾設置2個自動流量計（圖2）。

圖2 氮磷濃度與流量觀測點所在位置Figure 2 The location of nitrogen and phosphorus concentration and flow observation points

從2021 年7 月至11 月退水期間，在各水質監測點每日連續觀測溝道水體總氮（TN）和總磷（TP）濃度，采用清洗過的采水器在溝道中央水面下20 cm 處采集溝道水樣，采集完的水樣立即轉移到預先清洗過的PE 取樣瓶中，過0.45μm 濾膜后送至實驗室分析。水樣TN 濃度采用總有機碳/總氮分析儀（Multi N/C 3100，德國）分析，TP 濃度采用電感耦合等離子體發射光譜儀（Optima 8000，美國）分析。其中農溝監測點數據用于模型參數標定，支溝監測點數據用于模型結果驗證。自動流量計每5 min 記錄一次流量，然后通過2 個自動流量計差值確定區域退水量，用于退水量模擬結果驗證。

1.2 模型構建方法

面源氮磷流失量取決于退水量和退水中氮磷濃度。本研究基于水量平衡與土壤物理化學吸附理論，分別構建寧夏引黃灌區旱地與水田退水量、退水中氮磷濃度模型，估算農田氮磷流失量及流失系數。退水量模擬按旱地（玉米、小麥）和水田（水稻）的不同灌溉模式，利用水量平衡方程分別計算。農田退水中氮磷濃度是通過構建肥料施用量、退水發生時間與土壤氮磷平衡濃度模型進行模擬。

1.2.1 模型原理與構建

一般而言，徑流和淋溶是農田養分流失的基本途徑，對于東部平原地區，降水是造成農業面源污染的主要驅動力，通過地表徑流和地下滲漏進入水體。而在寧夏引黃灌區，由于常年干旱少雨，地表蒸發強烈，一般降雨難以有效發生產流。農田氮磷主要在灌溉水的驅動下，隨退水、淋溶和側滲等途徑進入受納水體。故構建模型（1）計算區域農田氮磷流失量：

式中：L為區域農田氮磷流失總量，kg；Ai為耕作類型i的灌溉面積，hm2；R為退水量，m3·hm-2；C為退水中氮磷濃度，mg·L-1；0.001為量綱轉換系數。

（1）退水量計算

由于旱地（玉米、小麥）是分不同灌水期進行灌溉，而水田（水稻）在全生育期均有灌水，故針對旱地和水田不同的灌溉模式，通過水量平衡方程分別計算退水量。

取單位面積（hm2）的包氣帶柱體作為水量平衡計算單元（圖3），上界面為地表，下界面選為潛水面，土壤水量的輸入項包括降水量、灌溉水量、側向補給量和深層補給量。土壤水量的輸出項包括蒸發量、側向排泄量、深層滲漏量以及植物的蒸騰量。由于灌區土壤水水平流動十分微弱，故不考慮側向補給和排泄[15]。由此退水量可以用如下模型表示：

圖3 土壤水量平衡概念圖Figure 3 A conceptual figure of soil water balance

式中：R為退水量，mm；P為降水量，mm；I為灌溉水量，mm；E為土壤蒸發量，mm；T為作物蒸騰量，mm；Du為地下水補給量，mm；Dd為深層滲漏量，mm；ΔSW為土壤儲水量變化，mm。

由于稻田大部分土壤水分飽和，地下水位高，土壤水分變化小，因此水分平衡只考慮田面蒸發、作物蒸騰（合為蒸散ET）和灌水、降雨[16]，退水量通過下式計算：

只有在稻田成熟期灌水結束后，由于蒸散和重力作用且缺乏灌水補充導致地下水位迅速下降，土壤含水量急劇下降，此時退水量為0。

（2）氮磷濃度模擬

退水中的氮磷濃度與肥料施用量、退水發生時間及土壤平衡濃度相關。在施肥初期，土壤中氮磷濃度最高，但是此時土壤對氮磷的吸附也最快，氮磷轉化過程最快。隨著土壤吸附量的增加，吸附速率逐步下降，氮磷轉化過程和吸附過程逐步變慢，最后達到一種平衡狀態（圖4），整個過程可以用冪指數下降函數模擬[17]：

圖4 土壤中氮磷濃度變化過程概念化模型Figure 4 A conceptual model for the change process of nitrogen and phosphorus concentrations in soil

式中：C為退水中氮磷濃度，mg·L-1，F為施肥量，kg·hm-2；a1、d1為參數；t為產流距離灌溉的時間，d；c1為平衡濃度，mg·L-1。

上述模型中a1、c1和d1受管理措施、土壤性質和土地利用的影響，因此該模型反映了農田管理，如施肥和灌溉、降雨、土壤條件等綜合過程，并且具有參數少、物理機制明確的特點，通過參數標定，可以用于估算灌區農田氮磷流失量。

1.2.2 模型參數獲取

玉米于6月8日進行春灌，7月5日和7月27日進行兩次夏澆，11月7日進行冬灌。水稻在生育期的前中期需水量較大，全生育期需要進行大量灌水以維持水稻正常生長發育的需求。葉廷平等[18]關于引黃灌區農作物地面灌溉制度的研究表明，小麥于4月25日灌春灌頭水，5 月 10 日、5 月 25 日灌 2、3 水，6 月 15 日灌4 水，11 月7 日進行冬灌。各次灌水情況見表1。根據王靜等[19]的研究（表2），研究區3 種作物全生育期灌水量均大于需水量，灌溉后都將產生退水。

表1 3種主要作物灌水量（mm）Table 1 Irrigation amount of the three main crops（mm）

表2 寧夏引黃灌區主要作物全生育期灌溉需水量Table 2 Irrigation water requirement of main crops in entire growth stages in the Yellow River irrigation area of Ningxia

對于玉米地，根據楊麗虎等[15]的研究和本研究關于降水量的監測，水量平衡參數如表3 所示。對于小麥地，參考王靜等[19]關于小麥生育期需水量的研究以及葉廷平等[18]關于引黃灌區農作物地面灌溉制度的研究，制定水量平衡參數（表3）。對于稻田系統，按照生育期計算土壤水分平衡和退水量，由于不同年度水稻灌溉量和蒸散量比較固定，因此引用易軍[16]在寧夏引黃灌區稻田的試驗結果，具體見表4。

表3 寧夏引黃灌區玉米和小麥灌水水量平衡參數（mm）Table 3 Parameters of maize and wheat irrigation water balance in the Yellow River irrigation area of Ningxia（mm）

表4 寧夏引黃灌區水稻不同生育期水量平衡參數（mm）Table 4 Parameters of water balance in different growth stages of rice in the yellow river irrigation area of Ningxia（mm）

1.3 不確定性分析

采用蒙特卡洛方法分析模型（1）的不確定性[20]。作為一種常見的模型不確定性分析方法，蒙特卡洛法可將參數的不確定性轉化為模型結果的不確定性[21]。本研究根據已有數據的統計分析，假設退水量服從均勻分布，退水氮磷濃度服從標準正態分布，應用Excel生成服從其分布規律的隨機數各400 個作為參數樣本，分析3 種作物氮磷流失模型計算結果的不確定性范圍。蒙特卡洛模擬所需參數變化范圍見表5。

表5 蒙特卡洛模擬所需參數變化范圍Table 5 Variability of model input parameters for Monte Carlo simulation

2 結果與分析

2.1 模型參數標定與驗證

2.1.1 氮磷濃度模擬參數標定與驗證

應用7 月28 日至8 月14 日農溝水質連續監測值標定模型（4）中a1、d1和c1（圖5）。對于水體中TN 濃度，a1為 5.0、d1為 0.15、c1為 0.5，對于水體中 TP 濃度，a1為0.9、d1為0.23、c1為0.05，模型解釋度分別為0.74和0.49，可用于水體中氮磷濃度模擬。

圖5 氮磷濃度模型參數標定Figure 5 Model parameter calibration

2.1.2 退水量與氮磷流失量模擬結果驗證

應用觀測區水自動流量計及支溝水質監測點數據，與冬灌退水量、氮磷流失量的模擬結果進行對比。模型驗證結果如圖6 所示，退水量、TN 流失量和總磷流失量的模擬誤差分別為11%、3%和12%，模型模擬精度較高。

圖6 模擬結果驗證Figure 6 Simulation result verification

2.2 典型種植業氮磷流失量時間變異模擬

將水量平衡與氮磷濃度模型和農田氮磷流失量模型進行耦合，分別估算玉米、小麥和水稻氮磷流失量（圖7）及流失系數。玉米各次灌水退水TN 總流失量（以N計，下同）為3.98 kg·hm-2，流失系數為0.74%，TP 流失量合計為0.17 kg·hm-2，流失系數為0.04%，其中春灌分別貢獻了84.2%的TN流失和76.5%的TP流失，遠高于其他各次灌水，是氮磷流失的高風險期。

圖7 3種作物單位面積氮磷流失量時間變化Figure 7 Time variation of nitrogen and phosphorus loss per unit area of three crops

小麥各次灌水退水氮磷流失趨勢與玉米基本一致 ，TN 流 失 量 合 計 為 5.12 kg·hm-2，流 失 系 數 為1.74%，TP 流失量合計為 0.65 kg·hm-2，流失系數為0.40%，其中春灌頭水的氮磷流失量最高，分別貢獻了73.6%的TN流失和84.6%的TP流失。

水稻全生育期及冬季灌水TN 流失量合計為20.46 kg·hm-2，流失系數為 5.76%，TP 流失量合計為1.03 kg·hm-2，流失系數為0.38%，其中分蘗期氮磷流失量最大，返青期次之，二者共占79.6%的TN 流失和61.4%的TP 流失。其他生育期灌水氮磷流失量較低，自孕穗期開始氮磷流失逐漸減少。

2.3 寧夏引黃灌區種植業氮磷流失量空間分布估算

結合各區縣3 種主要作物種植面積模擬各行政單元農田TN、TP 流失量空間分布及成分組成如圖8所示。全灌區種植業TN 總流失量（以N 計，下同）為887.51 t，玉米、小麥和水稻分別貢獻了25%、8%和67%，平均氮流失系數為1.99%；TP 總流失量為48.05 t，玉米、小麥和水稻分別貢獻了19%、18%和63%，平均磷流失系數為0.15%。從氮磷流失空間分布可以看出，TN 和TP 的流失分布情況基本一致，均呈現以金鳳區為中心向四周遞增的空間分布特征。從氮磷流失量角度分析，銀川市種植業氮磷流失量及流失系數為三市中最高，TN 流失量為431.6 t，氮流失系數為2.15%，TP流失量為23.4 t，磷流失系數為0.16%，其流失量占比接近灌區氮磷流失總量的50%。種植業氮磷流失最大的區縣為平羅縣，其TN 流失量為316 t，TP 流失量為17.4 t，其次為賀蘭縣和永寧縣，流失量最小的區縣為大武口區，其氮磷流失分別僅占全灌區總流失量的0.3%和0.4%；氮磷流失系數最大的區縣為興慶區，系數分別為3.69%和0.26%，氮流失系數最小的區縣為惠農區，為1.06%，磷流失系數最小的區縣為利通區，僅為0.10%。從氮磷流失組成成分來看，惠農區的玉米和大武口區的水稻在當地農田TN流失中占主要貢獻，占比分別為57%和45%?；蒉r區和大武口區水稻在當地農田TP 流失中占主要貢獻，占比分別為49%和47%，其他9 個區縣的水稻在當地農田氮磷流失中的貢獻均超過了50%。

進一步將上述灌區農田氮磷流失量分配到1 km分辨率柵格圖像上，得到引黃灌區農田氮磷流失強度空間格局圖。如圖9 所示，灌區農田TN 和TP 流失強度空間格局基本一致，永寧縣中部以及河東地區東部是氮磷流失的熱點區域。

圖9 寧夏引黃灌區農田氮磷流失強度空間格局（1 km）Figure 9 Spatial pattern of nitrogen and phosphorus loss intensity in farmland of the Yellow River irrigation area of Ningxia

2.4 模型不確定性分析

由于退水量和退水氮磷濃度的模擬誤差或波動會導致氮磷流失量發生變化，由水量平衡方程（2）和方程（3）及氮磷濃度模型（4），可確定氮磷流失量模型（1）的參數及流失量不確定性范圍，如表6 所示。玉米、小麥和水稻置信度為95%的氮素流失量（以N計）置信區間分別為[1.19，1.33]kg·hm-2、[1.18，1.33]kg·hm-2和[2.68，3.12]kg·hm-2，磷素流失量置信區間分別為[0.07，0.08]kg·hm-2、[0.13，0.15]kg·hm-2和[0.16，0.19]kg·hm-2。

表6 3種作物氮磷流失量、退水量和退水氮磷濃度的不確定性Table 6 Uncertainties of nitrogen and phosphorus loss，drainage water volume and nitrogen and phosphorus concentration of drainage water in the three crops

3 討論

春季既是農作物的重要生長期，也是面源污染流失的高風險期[22]。玉米春灌氮磷流失分別占其總流失量的84.2%和76.5%，小麥春灌頭水氮磷流失分別占其總流失量的73.6%和84.6%。水稻返青分蘗期氮磷流失分別占總流失的79.6%和61.4%，這與陳玲等[23]關于稻季播種插栽期至分蘗盛期是氮磷流失高風險期的研究結論一致。這是因為春季田間施肥量較大，灌水初期較低的氣溫對水體氮磷的降解有一定抑制作用[24]，再加上農作物處于生育期前期，作物肥料利用少，因此水體氮磷濃度高。同時，灌溉量大導致地下水位偏高，土壤深層滲漏量小，退水和氮磷流失也隨之增加[8]。

引黃灌區氮磷平均流失系數為1.99%和0.15%，遠低于南方濕潤平原區稻-麥系統（氮素的徑流損失系數為9.1%~22.5%，磷素的徑流損失系數為0.948%~3.394%）[25]。不確定性分析表明，退水量是氮磷流失的主要影響因素。而寧夏引黃灌區氣候干旱、降水匱乏，農田灌溉產生的退水流失遠低于南方濕潤區主要由降水驅動的徑流損失。磷素流失系數普遍低于氮素說明磷素較氮素更難流失，這是因為寧夏土壤具有很強的固磷能力[26]。土壤中磷素隨退水流失的主要載體是泥沙[27]，磷肥進入土壤后難以溶解，溶解后的磷素容易被土壤團聚體吸附，形成穩定結構，不易在土壤中遷移。

氮磷流失空間分布顯示，引黃灌區三市中，銀川市種植業氮磷污染流失量及流失系數最大，流失量占比接近灌區氮磷流失總量的50%，張愛平[4]在引黃灌區運用平均濃度法估算灌區農田退水氮磷負荷的研究中也有相似的結論。氮磷流失最大的區縣為平羅縣，其次為賀蘭縣和永寧縣，主要原因是三縣占有的耕地面積大，分別占全灌區耕地的34%、14%和19%。各區縣農田氮磷流失來源多以稻田流失為主，其中有9 個區縣水稻氮磷流失貢獻超過了50%，這與張學軍等[28]的結論相近。永寧縣中部以及河東地區東部是氮磷流失的熱點區域，這說明兩地氮磷流失風險更大，面源污染治理形勢嚴峻。

4 結論

（1）本研究根據寧夏灌區種植業水田和旱地的流失過程，建立了農田氮磷流失量的輕簡化模型，反映了農田管理、降水和土壤條件等過程對退水量和退水中氮磷濃度的影響。該模型參數少，物理機制明確。

（2）應用該模型分別模擬了灌區種植業主要作物（玉米、小麥和水稻）氮磷流失量的時間變異。各作物退水氮磷流失量的時間變異大，春季是面源污染流失的高風險期。

（3）應用該模型進一步模擬了寧夏灌區種植業氮磷流失空間分布，全灌區總氮總流失量（以N 計）為887.51 t，平均氮流失系數為1.99%；總磷總流失量為48.05 t，平均磷流失系數為0.15%。種植業氮磷流失量最大的區縣為平羅縣，其次為賀蘭縣和永寧縣，流失量較小的區縣為金鳳區和大武口區。