佘冬立,胡磊,夏永秋,劉文娟,李虹,馬建軍,馬琨
(1.河海大學農(nóng)業(yè)科學與工程學院,南京 211100;2.中國科學院南京土壤研究所,南京 210008;3.寧夏大學農(nóng)學院,銀川750021;4.寧夏農(nóng)業(yè)環(huán)境保護監(jiān)測站,銀川 750002;5.寧夏大學生態(tài)環(huán)境學院,銀川 750021;6.寧夏大學西北土地退化與生態(tài)恢復國家重點實驗室培育基地/西北退化生態(tài)系統(tǒng)恢復與重建教育部重點實驗室,銀川 750021)
寧夏引黃灌區(qū)位于我國西北內(nèi)陸干旱半干旱地區(qū),是我國四大古灌區(qū)之一,糧食產(chǎn)量和農(nóng)業(yè)產(chǎn)值占自治區(qū)的三分之二以上,在區(qū)域農(nóng)業(yè)經(jīng)濟發(fā)展中發(fā)揮著重要作用。灌區(qū)農(nóng)業(yè)生產(chǎn)集約化程度高、管理方式粗放,大量氮磷污染物質(zhì)隨退水經(jīng)各級排水溝排入黃河[1],對黃河上游水質(zhì)具有直接影響[2]。近年來,隨著國家對農(nóng)田化肥和規(guī)?;B(yǎng)殖投入的增加以及農(nóng)村生活污水盲目排放規(guī)模的擴大,寧夏引黃灌區(qū)面源污染呈現(xiàn)加劇的趨勢[3],以農(nóng)田為核心的灌區(qū)種植業(yè)退水污染是主要污染負荷來源[4],春灌和冬灌是灌區(qū)種植業(yè)面源污染流失高發(fā)期。大量施肥后,農(nóng)田土壤通過退水向周圍水體遷移排放氮磷等污染物質(zhì),經(jīng)過沿程農(nóng)溝、斗溝、支溝消納后,最終進入干溝和黃河,在造成土壤養(yǎng)分流失的同時,也加劇了地表水和地下水水質(zhì)惡化的風險。
農(nóng)業(yè)面源污染過程難以直接監(jiān)測,模型成為解析農(nóng)業(yè)面源污染負荷的主要工具。比較著名的過程模型是由美國農(nóng)業(yè)部農(nóng)業(yè)研究所開發(fā)的AGNPS[5]和美國國家環(huán)保局開發(fā)的SWAT和BASINS等模型[6],經(jīng)驗模型有Johnes[7]輸出系數(shù)模型。目前,引黃灌區(qū)農(nóng)業(yè)面源污染模型研究較少。胡亞偉等[8]以青銅峽灌區(qū)第一排水溝為典型小流域,結合3S 技術,根據(jù)灌區(qū)水系及灌溉特點對SWAT 模型進行改進,建立了寧夏引黃灌區(qū)面源污染負荷模型,但上述過程模型本質(zhì)上是建立在西方以旱作為主的農(nóng)業(yè)流域基礎之上,其過程模擬基于超滲產(chǎn)流理論,不能反映帶有田埂水田退水的實際情況,并且模擬所需參數(shù)眾多,獲取難度大。李強坤等[9]基于單元分析的觀點,提出了負荷貢獻率的概念,建立了寧夏引黃灌區(qū)農(nóng)業(yè)非點源改進輸出系數(shù)估算模型,但此類經(jīng)驗模型氮磷輸出系數(shù)的確定依賴于點尺度實驗或者樣地調(diào)查,很難體現(xiàn)不同區(qū)域、不同時段輸出系數(shù)的差異,并且不能反映農(nóng)業(yè)管理措施與面源負荷的響應關系,應用性與擴展性較差。因此,面對引黃灌區(qū)水力聯(lián)系的復雜性和不確定性,為了降低模型應用難度、提高農(nóng)業(yè)管理效率,亟需構建適合寧夏引黃灌區(qū)面源污染流失特點的輕簡化模型。
本研究基于寧夏引黃灌區(qū)典型種植業(yè)面源污染流失過程的理論依據(jù)和觀測數(shù)據(jù),構建了水田和旱地氮磷流失量模型,揭示了灌區(qū)種植業(yè)氮磷流失時空特征與空間格局,為改善灌區(qū)溝道及黃河水質(zhì)環(huán)境提供科學有效的措施,對實現(xiàn)寧夏回族自治區(qū)社會經(jīng)濟可持續(xù)發(fā)展具有科學指導意義。
寧夏引黃灌區(qū)(37°49′~39°23′ N,105°37′~106°39′E)位于寧夏回族自治區(qū)中北部,東至鄂爾多斯臺地西緣,西至賀蘭山東,南至寧夏青銅峽水利樞紐,北至石嘴山,屬于溫帶大陸性干旱氣候,全年干旱少雨。灌區(qū)含括了銀川市、石嘴山市、吳忠市在內(nèi)的3 個地級市及永寧、惠農(nóng)、青銅峽等11 個區(qū)縣(圖1)。灌溉面積為 33 萬 hm2[11],占總土地面積(6 239 km2)[10]的53%,2021 年引黃水量為 37.58 億 m3,排水量為 16.17億m3,多年平均引黃水量為40.95 億m3,排水量為22.09 億m3。灌區(qū)北高南低,地勢較為平坦,坡降約1/4 000[12]。
圖1 寧夏引黃灌區(qū)位置及土地利用[13]Figure 1 The location and land use of the Yellow River irrigation area of Ningxia[13]
面源污染過程觀測主要在引黃灌區(qū)境內(nèi)寧夏大學實驗農(nóng)場開展。實驗農(nóng)場位于寧夏中部黃河沖積平原永寧縣境內(nèi),該地年降水量180~200 mm,年蒸發(fā)量高達2 800~3 200 mm,年均氣溫8.5 ℃,年均≥10 ℃積溫3 300 ℃,無霜期140~160 d,平均日照時數(shù)3 000 h,日溫差13 ℃。觀測區(qū)包含兩塊地塊及兩條主要溝道,西邊玉米地塊面積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]。稻田農(nóng)溝長254 m,分布有3個水質(zhì)監(jiān)測點,兩塊玉米地之間的支溝長500 m,分布有4個水質(zhì)監(jiān)測點,首尾設置2個自動流量計(圖2)。
圖2 氮磷濃度與流量觀測點所在位置Figure 2 The location of nitrogen and phosphorus concentration and flow observation points
從2021 年7 月至11 月退水期間,在各水質(zhì)監(jiān)測點每日連續(xù)觀測溝道水體總氮(TN)和總磷(TP)濃度,采用清洗過的采水器在溝道中央水面下20 cm 處采集溝道水樣,采集完的水樣立即轉(zhuǎn)移到預先清洗過的PE 取樣瓶中,過0.45μm 濾膜后送至實驗室分析。水樣TN 濃度采用總有機碳/總氮分析儀(Multi N/C 3100,德國)分析,TP 濃度采用電感耦合等離子體發(fā)射光譜儀(Optima 8000,美國)分析。其中農(nóng)溝監(jiān)測點數(shù)據(jù)用于模型參數(shù)標定,支溝監(jiān)測點數(shù)據(jù)用于模型結果驗證。自動流量計每5 min 記錄一次流量,然后通過2 個自動流量計差值確定區(qū)域退水量,用于退水量模擬結果驗證。
面源氮磷流失量取決于退水量和退水中氮磷濃度。本研究基于水量平衡與土壤物理化學吸附理論,分別構建寧夏引黃灌區(qū)旱地與水田退水量、退水中氮磷濃度模型,估算農(nóng)田氮磷流失量及流失系數(shù)。退水量模擬按旱地(玉米、小麥)和水田(水稻)的不同灌溉模式,利用水量平衡方程分別計算。農(nóng)田退水中氮磷濃度是通過構建肥料施用量、退水發(fā)生時間與土壤氮磷平衡濃度模型進行模擬。
1.2.1 模型原理與構建
一般而言,徑流和淋溶是農(nóng)田養(yǎng)分流失的基本途徑,對于東部平原地區(qū),降水是造成農(nóng)業(yè)面源污染的主要驅(qū)動力,通過地表徑流和地下滲漏進入水體。而在寧夏引黃灌區(qū),由于常年干旱少雨,地表蒸發(fā)強烈,一般降雨難以有效發(fā)生產(chǎn)流。農(nóng)田氮磷主要在灌溉水的驅(qū)動下,隨退水、淋溶和側(cè)滲等途徑進入受納水體。故構建模型(1)計算區(qū)域農(nóng)田氮磷流失量:
式中:L為區(qū)域農(nóng)田氮磷流失總量,kg;Ai為耕作類型i的灌溉面積,hm2;R為退水量,m3·hm-2;C為退水中氮磷濃度,mg·L-1;0.001為量綱轉(zhuǎn)換系數(shù)。
(1)退水量計算
由于旱地(玉米、小麥)是分不同灌水期進行灌溉,而水田(水稻)在全生育期均有灌水,故針對旱地和水田不同的灌溉模式,通過水量平衡方程分別計算退水量。
取單位面積(hm2)的包氣帶柱體作為水量平衡計算單元(圖3),上界面為地表,下界面選為潛水面,土壤水量的輸入項包括降水量、灌溉水量、側(cè)向補給量和深層補給量。土壤水量的輸出項包括蒸發(fā)量、側(cè)向排泄量、深層滲漏量以及植物的蒸騰量。由于灌區(qū)土壤水水平流動十分微弱,故不考慮側(cè)向補給和排泄[15]。由此退水量可以用如下模型表示:
圖3 土壤水量平衡概念圖Figure 3 A conceptual figure of soil water balance
式中:R為退水量,mm;P為降水量,mm;I為灌溉水量,mm;E為土壤蒸發(fā)量,mm;T為作物蒸騰量,mm;Du為地下水補給量,mm;Dd為深層滲漏量,mm;ΔSW為土壤儲水量變化,mm。
由于稻田大部分土壤水分飽和,地下水位高,土壤水分變化小,因此水分平衡只考慮田面蒸發(fā)、作物蒸騰(合為蒸散ET)和灌水、降雨[16],退水量通過下式計算:
只有在稻田成熟期灌水結束后,由于蒸散和重力作用且缺乏灌水補充導致地下水位迅速下降,土壤含水量急劇下降,此時退水量為0。
(2)氮磷濃度模擬
退水中的氮磷濃度與肥料施用量、退水發(fā)生時間及土壤平衡濃度相關。在施肥初期,土壤中氮磷濃度最高,但是此時土壤對氮磷的吸附也最快,氮磷轉(zhuǎn)化過程最快。隨著土壤吸附量的增加,吸附速率逐步下降,氮磷轉(zhuǎn)化過程和吸附過程逐步變慢,最后達到一種平衡狀態(tài)(圖4),整個過程可以用冪指數(shù)下降函數(shù)模擬[17]:
圖4 土壤中氮磷濃度變化過程概念化模型Figure 4 A conceptual model for the change process of nitrogen and phosphorus concentrations in soil
式中:C為退水中氮磷濃度,mg·L-1,F(xiàn)為施肥量,kg·hm-2;a1、d1為參數(shù);t為產(chǎn)流距離灌溉的時間,d;c1為平衡濃度,mg·L-1。
上述模型中a1、c1和d1受管理措施、土壤性質(zhì)和土地利用的影響,因此該模型反映了農(nóng)田管理,如施肥和灌溉、降雨、土壤條件等綜合過程,并且具有參數(shù)少、物理機制明確的特點,通過參數(shù)標定,可以用于估算灌區(qū)農(nóng)田氮磷流失量。
1.2.2 模型參數(shù)獲取
玉米于6月8日進行春灌,7月5日和7月27日進行兩次夏澆,11月7日進行冬灌。水稻在生育期的前中期需水量較大,全生育期需要進行大量灌水以維持水稻正常生長發(fā)育的需求。葉廷平等[18]關于引黃灌區(qū)農(nóng)作物地面灌溉制度的研究表明,小麥于4月25日灌春灌頭水,5 月 10 日、5 月 25 日灌 2、3 水,6 月 15 日灌4 水,11 月7 日進行冬灌。各次灌水情況見表1。根據(jù)王靜等[19]的研究(表2),研究區(qū)3 種作物全生育期灌水量均大于需水量,灌溉后都將產(chǎn)生退水。
表1 3種主要作物灌水量(mm)Table 1 Irrigation amount of the three main crops(mm)
表2 寧夏引黃灌區(qū)主要作物全生育期灌溉需水量Table 2 Irrigation water requirement of main crops in entire growth stages in the Yellow River irrigation area of Ningxia
對于玉米地,根據(jù)楊麗虎等[15]的研究和本研究關于降水量的監(jiān)測,水量平衡參數(shù)如表3 所示。對于小麥地,參考王靜等[19]關于小麥生育期需水量的研究以及葉廷平等[18]關于引黃灌區(qū)農(nóng)作物地面灌溉制度的研究,制定水量平衡參數(shù)(表3)。對于稻田系統(tǒng),按照生育期計算土壤水分平衡和退水量,由于不同年度水稻灌溉量和蒸散量比較固定,因此引用易軍[16]在寧夏引黃灌區(qū)稻田的試驗結果,具體見表4。
表3 寧夏引黃灌區(qū)玉米和小麥灌水水量平衡參數(shù)(mm)Table 3 Parameters of maize and wheat irrigation water balance in the Yellow River irrigation area of Ningxia(mm)
表4 寧夏引黃灌區(qū)水稻不同生育期水量平衡參數(shù)(mm)Table 4 Parameters of water balance in different growth stages of rice in the yellow river irrigation area of Ningxia(mm)
采用蒙特卡洛方法分析模型(1)的不確定性[20]。作為一種常見的模型不確定性分析方法,蒙特卡洛法可將參數(shù)的不確定性轉(zhuǎn)化為模型結果的不確定性[21]。本研究根據(jù)已有數(shù)據(jù)的統(tǒng)計分析,假設退水量服從均勻分布,退水氮磷濃度服從標準正態(tài)分布,應用Excel生成服從其分布規(guī)律的隨機數(shù)各400 個作為參數(shù)樣本,分析3 種作物氮磷流失模型計算結果的不確定性范圍。蒙特卡洛模擬所需參數(shù)變化范圍見表5。
表5 蒙特卡洛模擬所需參數(shù)變化范圍Table 5 Variability of model input parameters for Monte Carlo simulation
2.1.1 氮磷濃度模擬參數(shù)標定與驗證
應用7 月28 日至8 月14 日農(nóng)溝水質(zhì)連續(xù)監(jiān)測值標定模型(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 氮磷濃度模型參數(shù)標定Figure 5 Model parameter calibration
2.1.2 退水量與氮磷流失量模擬結果驗證
應用觀測區(qū)水自動流量計及支溝水質(zhì)監(jiān)測點數(shù)據(jù),與冬灌退水量、氮磷流失量的模擬結果進行對比。模型驗證結果如圖6 所示,退水量、TN 流失量和總磷流失量的模擬誤差分別為11%、3%和12%,模型模擬精度較高。
圖6 模擬結果驗證Figure 6 Simulation result verification
將水量平衡與氮磷濃度模型和農(nóng)田氮磷流失量模型進行耦合,分別估算玉米、小麥和水稻氮磷流失量(圖7)及流失系數(shù)。玉米各次灌水退水TN 總流失量(以N計,下同)為3.98 kg·hm-2,流失系數(shù)為0.74%,TP 流失量合計為0.17 kg·hm-2,流失系數(shù)為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,流 失 系 數(shù) 為1.74%,TP 流失量合計為 0.65 kg·hm-2,流失系數(shù)為0.40%,其中春灌頭水的氮磷流失量最高,分別貢獻了73.6%的TN流失和84.6%的TP流失。
水稻全生育期及冬季灌水TN 流失量合計為20.46 kg·hm-2,流失系數(shù)為 5.76%,TP 流失量合計為1.03 kg·hm-2,流失系數(shù)為0.38%,其中分蘗期氮磷流失量最大,返青期次之,二者共占79.6%的TN 流失和61.4%的TP 流失。其他生育期灌水氮磷流失量較低,自孕穗期開始氮磷流失逐漸減少。
結合各區(qū)縣3 種主要作物種植面積模擬各行政單元農(nóng)田TN、TP 流失量空間分布及成分組成如圖8所示。全灌區(qū)種植業(yè)TN 總流失量(以N 計,下同)為887.51 t,玉米、小麥和水稻分別貢獻了25%、8%和67%,平均氮流失系數(shù)為1.99%;TP 總流失量為48.05 t,玉米、小麥和水稻分別貢獻了19%、18%和63%,平均磷流失系數(shù)為0.15%。從氮磷流失空間分布可以看出,TN 和TP 的流失分布情況基本一致,均呈現(xiàn)以金鳳區(qū)為中心向四周遞增的空間分布特征。從氮磷流失量角度分析,銀川市種植業(yè)氮磷流失量及流失系數(shù)為三市中最高,TN 流失量為431.6 t,氮流失系數(shù)為2.15%,TP流失量為23.4 t,磷流失系數(shù)為0.16%,其流失量占比接近灌區(qū)氮磷流失總量的50%。種植業(yè)氮磷流失最大的區(qū)縣為平羅縣,其TN 流失量為316 t,TP 流失量為17.4 t,其次為賀蘭縣和永寧縣,流失量最小的區(qū)縣為大武口區(qū),其氮磷流失分別僅占全灌區(qū)總流失量的0.3%和0.4%;氮磷流失系數(shù)最大的區(qū)縣為興慶區(qū),系數(shù)分別為3.69%和0.26%,氮流失系數(shù)最小的區(qū)縣為惠農(nóng)區(qū),為1.06%,磷流失系數(shù)最小的區(qū)縣為利通區(qū),僅為0.10%。從氮磷流失組成成分來看,惠農(nóng)區(qū)的玉米和大武口區(qū)的水稻在當?shù)剞r(nóng)田TN流失中占主要貢獻,占比分別為57%和45%?;蒉r(nóng)區(qū)和大武口區(qū)水稻在當?shù)剞r(nóng)田TP 流失中占主要貢獻,占比分別為49%和47%,其他9 個區(qū)縣的水稻在當?shù)剞r(nóng)田氮磷流失中的貢獻均超過了50%。
進一步將上述灌區(qū)農(nóng)田氮磷流失量分配到1 km分辨率柵格圖像上,得到引黃灌區(qū)農(nóng)田氮磷流失強度空間格局圖。如圖9 所示,灌區(qū)農(nóng)田TN 和TP 流失強度空間格局基本一致,永寧縣中部以及河東地區(qū)東部是氮磷流失的熱點區(qū)域。
圖9 寧夏引黃灌區(qū)農(nóng)田氮磷流失強度空間格局(1 km)Figure 9 Spatial pattern of nitrogen and phosphorus loss intensity in farmland of the Yellow River irrigation area of Ningxia
由于退水量和退水氮磷濃度的模擬誤差或波動會導致氮磷流失量發(fā)生變化,由水量平衡方程(2)和方程(3)及氮磷濃度模型(4),可確定氮磷流失量模型(1)的參數(shù)及流失量不確定性范圍,如表6 所示。玉米、小麥和水稻置信度為95%的氮素流失量(以N計)置信區(qū)間分別為[1.19,1.33]kg·hm-2、[1.18,1.33]kg·hm-2和[2.68,3.12]kg·hm-2,磷素流失量置信區(qū)間分別為[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
春季既是農(nóng)作物的重要生長期,也是面源污染流失的高風險期[22]。玉米春灌氮磷流失分別占其總流失量的84.2%和76.5%,小麥春灌頭水氮磷流失分別占其總流失量的73.6%和84.6%。水稻返青分蘗期氮磷流失分別占總流失的79.6%和61.4%,這與陳玲等[23]關于稻季播種插栽期至分蘗盛期是氮磷流失高風險期的研究結論一致。這是因為春季田間施肥量較大,灌水初期較低的氣溫對水體氮磷的降解有一定抑制作用[24],再加上農(nóng)作物處于生育期前期,作物肥料利用少,因此水體氮磷濃度高。同時,灌溉量大導致地下水位偏高,土壤深層滲漏量小,退水和氮磷流失也隨之增加[8]。
引黃灌區(qū)氮磷平均流失系數(shù)為1.99%和0.15%,遠低于南方濕潤平原區(qū)稻-麥系統(tǒng)(氮素的徑流損失系數(shù)為9.1%~22.5%,磷素的徑流損失系數(shù)為0.948%~3.394%)[25]。不確定性分析表明,退水量是氮磷流失的主要影響因素。而寧夏引黃灌區(qū)氣候干旱、降水匱乏,農(nóng)田灌溉產(chǎn)生的退水流失遠低于南方濕潤區(qū)主要由降水驅(qū)動的徑流損失。磷素流失系數(shù)普遍低于氮素說明磷素較氮素更難流失,這是因為寧夏土壤具有很強的固磷能力[26]。土壤中磷素隨退水流失的主要載體是泥沙[27],磷肥進入土壤后難以溶解,溶解后的磷素容易被土壤團聚體吸附,形成穩(wěn)定結構,不易在土壤中遷移。
氮磷流失空間分布顯示,引黃灌區(qū)三市中,銀川市種植業(yè)氮磷污染流失量及流失系數(shù)最大,流失量占比接近灌區(qū)氮磷流失總量的50%,張愛平[4]在引黃灌區(qū)運用平均濃度法估算灌區(qū)農(nóng)田退水氮磷負荷的研究中也有相似的結論。氮磷流失最大的區(qū)縣為平羅縣,其次為賀蘭縣和永寧縣,主要原因是三縣占有的耕地面積大,分別占全灌區(qū)耕地的34%、14%和19%。各區(qū)縣農(nóng)田氮磷流失來源多以稻田流失為主,其中有9 個區(qū)縣水稻氮磷流失貢獻超過了50%,這與張學軍等[28]的結論相近。永寧縣中部以及河東地區(qū)東部是氮磷流失的熱點區(qū)域,這說明兩地氮磷流失風險更大,面源污染治理形勢嚴峻。
(1)本研究根據(jù)寧夏灌區(qū)種植業(yè)水田和旱地的流失過程,建立了農(nóng)田氮磷流失量的輕簡化模型,反映了農(nóng)田管理、降水和土壤條件等過程對退水量和退水中氮磷濃度的影響。該模型參數(shù)少,物理機制明確。
(2)應用該模型分別模擬了灌區(qū)種植業(yè)主要作物(玉米、小麥和水稻)氮磷流失量的時間變異。各作物退水氮磷流失量的時間變異大,春季是面源污染流失的高風險期。
(3)應用該模型進一步模擬了寧夏灌區(qū)種植業(yè)氮磷流失空間分布,全灌區(qū)總氮總流失量(以N 計)為887.51 t,平均氮流失系數(shù)為1.99%;總磷總流失量為48.05 t,平均磷流失系數(shù)為0.15%。種植業(yè)氮磷流失量最大的區(qū)縣為平羅縣,其次為賀蘭縣和永寧縣,流失量較小的區(qū)縣為金鳳區(qū)和大武口區(qū)。