陳明洪，白 璐，劉宣冶，郭建松，汪超子

基于土壤微環(huán)境分層的平原水稻灌區(qū)磷污染模型

陳明洪，白 璐，劉宣冶，郭建松，汪超子

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

灌溉農(nóng)田產(chǎn)生的非點(diǎn)源磷是造成水體富營養(yǎng)化的主要原因之一，但是目前國內(nèi)外的磷污染模型對(duì)于平原灌區(qū)灌溉和排水管理下的水分運(yùn)動(dòng)過程和氧化還原條件下的磷轉(zhuǎn)化過程的定量表征還比較欠缺。該研究構(gòu)建了適用于平原水稻灌區(qū)水分運(yùn)動(dòng)和磷轉(zhuǎn)化遷移的機(jī)理性磷污染模型，模型根據(jù)稻田水量平衡和溝道運(yùn)動(dòng)波方程模擬灌區(qū)產(chǎn)匯流，采用考慮土壤微環(huán)境分層的磷轉(zhuǎn)化模型和對(duì)流擴(kuò)散方程模擬灌區(qū)產(chǎn)匯污。模型中將耕作層分為有氧層和無氧層，定量表征了由于水田干濕交替導(dǎo)致的土壤分層溶解氧變化和磷的轉(zhuǎn)化過程。為了驗(yàn)證模型的可靠性，利用黑龍江省和平灌區(qū)2018年試驗(yàn)田實(shí)測田間土壤水、積水及排水和2條支溝實(shí)測排水的水量和水質(zhì)數(shù)據(jù)對(duì)模型進(jìn)行了率定和驗(yàn)證。驗(yàn)證結(jié)果顯示試驗(yàn)田、一排和七排排水的徑流流量、磷濃度的模擬結(jié)果與實(shí)測結(jié)果都吻合較好。模擬排水流量的Nash-Sutcliffe效率系數(shù)（NSE）和決定系數(shù)（2）分別大于0.820和0.815；模擬總磷濃度的NSE和2分別大于0.811和0.821；考慮土壤微環(huán)境分層后得到的土壤可溶磷垂向分布結(jié)果比不考慮分層時(shí)與原位實(shí)測結(jié)果更接近。該磷污染模型被用于模擬和平水稻灌區(qū)的非點(diǎn)源磷污染。灌區(qū)磷污染濃度過程統(tǒng)計(jì)結(jié)果顯示，2018年水稻生育期內(nèi)通過排水和滲漏流失的磷為1.88 kg/hm2，約占施肥和灌溉磷輸入的5.7%。其中分蘗期和拔節(jié)孕穗期徑流磷輸出負(fù)荷最大，分別為0.85和0.60 kg/hm2；泡田期和分蘗期滲漏輸出負(fù)荷最大，分別為0.11和0.16 kg/hm2。灌區(qū)一排和四排由于控制面積大，輸出磷污染物總量大于其他排水溝。該研究所構(gòu)建的磷污染模型包含稻田灌溉和排水過程的水分運(yùn)動(dòng)及稻田干濕交替引起的氧化還原變化條件下的磷轉(zhuǎn)化和運(yùn)移，可為平原灌區(qū)水肥運(yùn)籌管理下的磷運(yùn)移模擬提供更準(zhǔn)確的方法。

水體污染；排水；氧化還原反應(yīng)；微環(huán)境分層；非點(diǎn)源磷污染；遷移轉(zhuǎn)化；灌區(qū)

0 引 言

水體污染已經(jīng)成為當(dāng)今社會(huì)人類面臨的主要環(huán)境問題之一，其中磷污染是造成水體富營養(yǎng)化的主要限制因子，為解決這一難題，各國都在控制營養(yǎng)物質(zhì)的來源上投入了大量的人力財(cái)力。大量的研究結(jié)果表明[1-2]，農(nóng)業(yè)非點(diǎn)源已經(jīng)成為當(dāng)今水環(huán)境的主要污染源。但農(nóng)業(yè)非點(diǎn)源污染具有形成過程隨機(jī)性大、影響因子復(fù)雜、分布范圍廣、危害大且產(chǎn)生機(jī)理模糊等特點(diǎn)，由此導(dǎo)致模擬和控制非點(diǎn)源污染難度較大。

定量化表征非點(diǎn)源污染負(fù)荷及過程是水環(huán)境污染治理的重要基礎(chǔ)工作。目前，模擬具有一定坡度的流域非點(diǎn)源污染模型比較完善，出現(xiàn)了SWAT（Soil and Water Assessment Tool）[3-4]、HSPF（Hydrological Simulation Program-Fortran）[5-6]、SWRRB（Simulator for Water Resources in Rural Basins）[7]、AGNPS（Agricultural Non-point Source Pollution）[8-9]、ANSWERS（A Real Nonpoint Source Watershed Environment Response Simulation）[10-11]等一系列模型，但針對(duì)平原灌區(qū)的非點(diǎn)源污染模型較少，大多數(shù)都是基于以上流域模型進(jìn)行改進(jìn)。仕玉治等[12]根據(jù)水量平衡方程和水稻田的灌溉方式修改了SWAT模型的自動(dòng)灌溉技術(shù)，模擬努敏河流域長閣灌區(qū)徑流量的效果較好；鄭捷等[13]在考慮灌區(qū)實(shí)際渠道的基礎(chǔ)上對(duì)子流域劃分以及作物耗水量計(jì)算模塊等方面進(jìn)行了改進(jìn)，成功地應(yīng)用于汾河灌區(qū)的水量平衡模擬；胡文慧等[14]根據(jù)汾河灌區(qū)實(shí)地的灌溉施肥調(diào)查資料以及實(shí)測的總氮、總磷、水文站監(jiān)測數(shù)據(jù)構(gòu)建SWAT模型，模擬了汾河灌區(qū)徑流和污染物負(fù)荷。SWAT模型雖然已有應(yīng)用于平原灌區(qū)的案例，但大多都采用模擬自然流域的SCS（Soil Conservation Service）徑流曲線法模擬產(chǎn)流，或者僅考慮農(nóng)田灌溉的人為影響調(diào)整灌溉模塊，并未考慮灌區(qū)作物實(shí)際排水模式，比如現(xiàn)行節(jié)水政策下水稻田的控制灌溉或干濕交替灌溉等。所以上述模型在平原水稻灌區(qū)人為管理下的水分運(yùn)動(dòng)過程模擬中仍顯能力不足，更影響了水文驅(qū)動(dòng)條件下的污染物遷移過程。

此外，對(duì)于磷污染物遷移模塊，一般非點(diǎn)源模型只考慮其質(zhì)量守恒和濃度平衡[15-16]，計(jì)算其平均輸出負(fù)荷，并沒有考慮磷的動(dòng)力遷移過程。然而水稻生育期內(nèi)干濕交替引起的土壤氧化還原條件變化較大，土壤磷素在有氧和無氧變化環(huán)境下的生化反應(yīng)過程差異較大，在上述模型中均未得到體現(xiàn)。因此，需要建立適合平原水稻灌區(qū)水分運(yùn)動(dòng)及磷污染過程機(jī)理的模型，以準(zhǔn)確定量地模擬磷在水稻生育期內(nèi)的遷移轉(zhuǎn)化過程，從而為評(píng)價(jià)平原水稻灌區(qū)磷污染對(duì)流域水環(huán)境的影響打下基礎(chǔ)。本研究建立基于土壤微環(huán)境分層的非點(diǎn)源磷污染模型，定量表征磷污染物在水土、土壤有氧和無氧分層界面的遷移轉(zhuǎn)化規(guī)律，并將模型應(yīng)用于黑龍江和平灌區(qū)非點(diǎn)源磷污染過程的模擬。

1 灌區(qū)磷污染模型原理及構(gòu)建

1.1 水文模塊

水分運(yùn)動(dòng)為污染物提供遷移的介質(zhì)和能量，水文路徑同時(shí)也是污染物遷移的路徑，合理并準(zhǔn)確描述水文過程直接影響整個(gè)模型模擬結(jié)果[17]。與SWAT、HSPF模型相似，根據(jù)灌區(qū)地形、土壤和水文特征將區(qū)域劃分為若干個(gè)水文響應(yīng)單元，在每個(gè)單元內(nèi)根據(jù)天然降雨和灌溉制度模擬地表水和土壤水的運(yùn)動(dòng)過程。

水稻田可能存在有和無淹沒水層2種情況。當(dāng)田面有淹沒水層且有排水時(shí)，田塊地表水流運(yùn)動(dòng)可用一維水動(dòng)力學(xué)水流運(yùn)動(dòng)控制方程[18]描述

1）田間灌溉

2）田間地表排水

當(dāng)田面有水層但無排水或無水層時(shí)，稻田土壤水運(yùn)動(dòng)過程均可簡化為二維垂向土壤水動(dòng)力學(xué)問題，模型利用Richards方程[20]計(jì)算

式中為土壤壓力水頭，cm；為垂向距離，向下為正，cm；()為土壤非飽和水力傳導(dǎo)度，cm/min；()為比水容量，cm。選用水流和溶質(zhì)運(yùn)移模擬軟件（HYDRUS）提供的van Genuchten-Mualem公式[21]計(jì)算土壤水力特性參數(shù)。

由于水稻田犁底層具有較強(qiáng)的阻水作用，本文暫不考慮犁底層向底土層和地下水的垂直下滲作用，只考慮犁底層的側(cè)滲作用。在耕作層僅考慮向犁底層的垂直入滲，側(cè)面為零通量邊界。田面存在淹沒水層時(shí)，()=0，()=s，s為土壤飽和水力傳導(dǎo)度，cm/min。無淹沒水層時(shí)，采用蒸發(fā)強(qiáng)度作為上邊界。

稻田地表排水根據(jù)排水溝道匯流順序匯入各級(jí)排水溝道，排水在溝道水流演進(jìn)采用略去慣性項(xiàng)和壓力項(xiàng)的一維圣維南方程[22]模擬。

1.2 磷的遷移轉(zhuǎn)化過程

水稻在不同生育期需水狀況不同，導(dǎo)致水田呈現(xiàn)干濕交替的水土環(huán)境，從而影響稻田的磷轉(zhuǎn)化運(yùn)移過程。干濕交替會(huì)造成稻田土壤水土界面有氧-無氧環(huán)境頻繁變化，而磷轉(zhuǎn)化的生物化學(xué)過程又受到有氧-無氧環(huán)境的強(qiáng)烈影響，因此需要考慮土壤和水層、有氧層與無氧層之間的界面效應(yīng)。根據(jù)土壤中含氧量差異將稻田土壤分為有氧層和無氧層。一般河道沉積物有氧層為地表以下3 mm厚[23]，但水稻土壤有氧層可以取1 cm[24]，無氧層取有氧層以下30 cm厚。在田面積水和無人為干擾的情況下，有氧層和無氧層的厚度保持不變，但田面無淹沒水層時(shí)，有氧層會(huì)隨著土壤含水層下降而下降。本研究提出的稻田土壤磷轉(zhuǎn)化過程見圖1。

注：R1為施肥磷通量；R2、R3分別為有氧層和上覆水、有氧層和無氧層的磷擴(kuò)散通量；R4為顆粒混摻磷通量；R5為稻田磷衰減通量；R6、R7分別為水體層向有氧層、有氧層向無氧層的磷滲透通量；R8為灌溉水?dāng)y帶磷通量。

1.2.1 磷在稻田的循環(huán)過程

由圖1可以看出，磷在稻田的循環(huán)由外界輸入、各層之間的交換、水層的生化衰減反應(yīng)3部分構(gòu)成。磷在各層之間的交換主要包括3個(gè)過程：水體之間的擴(kuò)散作用、生物擾動(dòng)造成的顆?；鞊?、滲透作用。

1）磷的來源主要是施肥、降雨灌溉外界帶入，施入土壤中的有機(jī)磷肥或動(dòng)植物殘?bào)w的有機(jī)磷在微生物的作用下會(huì)產(chǎn)生礦化分解變成無機(jī)磷，主要是正磷酸鹽。施肥和灌溉帶入的磷通量分別用1、8表示。

2）水體之間的擴(kuò)散作用。用2、3分別表示有氧層和上覆水、有氧層和無氧層之間的磷擴(kuò)散通量。擴(kuò)散是在濃度梯度作用下發(fā)生的，擴(kuò)散通量與濃度梯度成正比，表達(dá)式如下[23]

式中C為水中污染物濃度，mg/L；C為有氧層中污染物濃度，mg/L；C為無氧層中污染物濃度，mg/L；L01、L12分別是水體層與有氧層、有氧層與無氧層間的水體擴(kuò)散傳質(zhì)系數(shù)，m/s；d1、d2分別表示溶解磷在有氧層和無氧層中所占的比例。

3）有氧層和無氧層之間顆粒混摻作用。在生物擾動(dòng)作用下，有氧層和無氧層間的土壤顆粒發(fā)生混摻，引起磷交換。顆粒混摻交換通量用4表示[23]

式中12表示土壤顆粒的混摻速率，m/s；p1、p2分別表示吸附磷在有氧層和無氧層中所占的比例。

4）滲透作用。由于重力作用，水體層和有氧層、有氧層和無氧層之間會(huì)發(fā)生滲透，滲透水?dāng)y帶著磷向下遷移。有氧層和無氧層的滲透通量分別用6、7表示[23]

式中L01和L02分別是水體層與有氧層、有氧層與無氧層之間的滲透作用傳質(zhì)系數(shù)，m/s。

5）衰減作用。稻田水層中的磷由于自凈作用會(huì)隨著時(shí)間衰減，滿足一階反應(yīng)動(dòng)力學(xué)。衰減通量用5表示

式中f為衰減系數(shù)，d-1。

1.2.2 考慮土壤微環(huán)境分層的磷分配過程

磷在固液相的分配比例是模型的重要參數(shù)，磷的吸附解吸過程是分配比例確定的依據(jù)。土壤顆粒中的鐵（氫）氧化物表面存在著大量表面羥基，對(duì)無機(jī)磷的吸附作用很強(qiáng)，鐵及其氧化物的存在形式?jīng)Q定了磷在土壤中的形態(tài)[25]。在有氧條件下，溶解態(tài)的Fe2+被氧化為FeOOH(s)，對(duì)磷的吸附作用較強(qiáng)；在無氧或相對(duì)缺氧的條件下，F(xiàn)eOOH(s)被還原為Fe2+，對(duì)磷的吸附作用減弱?？梢娙芙饬缀臀搅自谕寥乐兴嫉谋壤齞j和pj是溶氧濃度的函數(shù)。

有氧條件

無氧條件

分配系數(shù)的最早提出是一個(gè)經(jīng)驗(yàn)常數(shù)，缺少基本的理論推導(dǎo)。自20世紀(jì)70年代初期Stumm等[26]提出表面絡(luò)合模型后，其被廣泛應(yīng)用于吸附過程的研究。表面絡(luò)合模型是基于顆粒表面電荷，用溶液中絡(luò)合平衡類似的方法處理顆粒物界面上的結(jié)合問題。它根據(jù)表面質(zhì)子反應(yīng)、與磷酸根的配體交換確定2個(gè)表面酸度常數(shù)和3個(gè)表面絡(luò)合常數(shù)，再根據(jù)表面絡(luò)合模型確定吸附態(tài)和溶解態(tài)的磷含量[25,27]，進(jìn)而可確定溶解相和顆粒相的比例dj

稻田表層土壤及水稻根際土壤有相對(duì)較高的溶解氧含量，而非根際土壤和深層土壤溶解氧含量極低[28]。劉錦濤等[29]研究稻田在常規(guī)灌溉和控制灌溉條件下土壤溶解氧的分布，發(fā)現(xiàn)溶解氧濃度在不同生育期和不同土壤深度均有所差異。因此，本文根據(jù)不同時(shí)期農(nóng)田水層深度及氣象資料計(jì)算表層水和土壤有氧層的溶氧濃度變化過程，依據(jù)土壤初始有效鐵含量確定各層土壤中FeOOH(s)的含量，使用絡(luò)合模型再根據(jù)式（12）確定磷在固液相的分配比例。表層水/土壤層氧氣濃度計(jì)算如下[30]

1.2.3 磷遷移轉(zhuǎn)化整體控制方程

上節(jié)內(nèi)容即為水稻田磷遷移轉(zhuǎn)化的過程及其定量表征，包括田面積水、有氧層和無氧層，由此得到磷的總控制方程，表達(dá)式如下

式中1和2分別為土壤有氧層和無氧層厚度，m；D為無氧層污染物擴(kuò)散系數(shù)，m2/s；為土壤層垂向坐標(biāo)，m。

污染物進(jìn)入溝道后應(yīng)用一維非保守性污染物遷移擴(kuò)散方程[31]模擬，基本方程如下

式中d為溝道水深，m；為溝道污染物濃度，mg/L；為溝道流速，m/s；為溝道污染物擴(kuò)散系數(shù)，m2/s；d為溝道綜合衰減系數(shù)，d-1；為溝道水平方向坐標(biāo)，m。

水稻田排水及溝道水體中如存在懸浮顆粒物，應(yīng)考慮顆粒物對(duì)磷的攜帶作用。一般采用分配系數(shù)K表征磷在顆粒物和水相的分配。水體總磷濃度計(jì)算如下[32]

式中C為水體總磷濃度，mg/L；為固體顆粒物濃度，kg/L；K為磷素的固液分配系數(shù)，L/kg。本研究在野外試驗(yàn)時(shí)采集了試驗(yàn)田及主要排水溝的水樣，實(shí)測懸浮固體顆粒物濃度極低，故后文模擬中忽略水體中固體含磷量。

本文所列偏微分方程均采用隱式有限差分方法進(jìn)行離散，并給出初始條件和邊界條件求解計(jì)算域內(nèi)所有變量。

2 研究區(qū)概況

和平灌區(qū)位于呼蘭河流域的呼蘭河沖積平原上，灌區(qū)大部分面積屬于慶安縣，少部分屬于鐵力市，灌區(qū)總面積約為6 700 hm2。灌區(qū)位置和主要水系如圖2所示，和平灌區(qū)和勞模灌區(qū)分別位于呼蘭河兩岸。和平灌區(qū)地理坐標(biāo)為127°17'～127°49'E，46°34'～47°07'N。慶安縣屬于中溫帶半干旱半濕潤的大陸性季風(fēng)氣候，灌區(qū)多年平均降水量545.3 mm，變化范圍在450～700 mm之間，降水年際年內(nèi)分配不均，空間分布趨勢(shì)是由東北向西南遞減。該區(qū)蒸發(fā)量及其時(shí)空分布變化較大，多年平均水面蒸發(fā)量為664.5 mm。

和平灌區(qū)的種植作物95%以上為水稻，灌區(qū)水稻生育期為5—9月，期間施基肥、分蘗肥和穗肥，具體情況見表1。該區(qū)為自流引水灌區(qū)，灌區(qū)從呼蘭河及安邦河、拉林清河等支流引水，灌區(qū)內(nèi)有1條干渠，18條支渠，10條排水支溝，但9排和10排位于井灌區(qū)，生育期基本無排水，本文模擬未考慮。

注：為保持匯水流域完整性，除和平灌區(qū)外，圖中流域還包括了呼蘭河西北側(cè)的勞模灌區(qū)。

表1 和平灌區(qū)水稻生育期及施肥制度

為了解灌區(qū)土壤背景屬性及其空間異質(zhì)性分布，在整個(gè)灌區(qū)中選取了44個(gè)典型采樣點(diǎn)，分別取表層0～20 cm土壤樣品回實(shí)驗(yàn)室進(jìn)行測定。土壤質(zhì)地和背景資料見表2。實(shí)測數(shù)據(jù)顯示，研究區(qū)土壤機(jī)械組成差異較小，根據(jù)美國制分級(jí)標(biāo)準(zhǔn)，灌區(qū)土壤基本為粉壤土。土壤總磷平均質(zhì)量分?jǐn)?shù)約為0.80～1.08 g/kg，東北部總磷背景含量較高，整體趨勢(shì)向西南遞減。土壤有效鐵的空間差異性較大。

表2 和平灌區(qū)農(nóng)田土壤的基本性質(zhì)

3 模型應(yīng)用

3.1 模型構(gòu)建與數(shù)據(jù)來源

灌區(qū)磷污染物模型采用的基礎(chǔ)數(shù)字高程模型（DEM）是SRTM 90 m精度的數(shù)據(jù)（圖3a）。借助GIS平臺(tái)根據(jù)DEM高程圖和人工burn-in[33]功能生成符合灌區(qū)實(shí)際渠系情況的河網(wǎng)，將灌區(qū)調(diào)研的實(shí)際河網(wǎng)疊加在已有的DEM上，保持河道所在柵格的高程值不變，而將其他非河道所在位置的柵格整體增加一個(gè)微小值，加大河道所在柵格和周邊柵格之間的坡度，使得河道所在柵格的匯水能力增強(qiáng)，提高河道提取的準(zhǔn)確程度。土地利用采用2010年黑龍江省遙感土地利用數(shù)據(jù)（圖3b）。土壤類型分布采用南京土壤所1∶100萬土壤空間數(shù)據(jù)，并利用灌區(qū)采樣數(shù)據(jù)進(jìn)行修正（圖3c）。將模型基礎(chǔ)數(shù)據(jù)導(dǎo)入ARCSWAT軟件劃分水文響應(yīng)單元和子流域，本流域共劃分出64個(gè)子流域（圖3d）。為保持匯水流域完整性，除和平灌區(qū)外，DEM及流域劃分還包括了呼蘭河西北側(cè)的勞模灌區(qū)（見圖2和圖3），但本研究數(shù)據(jù)分析僅含和平灌區(qū)模擬數(shù)據(jù)。氣象數(shù)據(jù)來自黑龍江省水科院提供的慶安縣氣象站資料，數(shù)據(jù)系列時(shí)間長度2000—2018 年，數(shù)據(jù)主要包括逐日的觀測數(shù)據(jù)，包括最高氣溫、最低氣溫、降水、相對(duì)濕度、風(fēng)速和日照時(shí)數(shù)等。此外，和平灌區(qū)試驗(yàn)田氣象數(shù)據(jù)采用試驗(yàn)站小型氣象站的2015—2018年數(shù)據(jù)。灌區(qū)各子流域灌溉數(shù)據(jù)利用和平灌區(qū)灌溉定額和輪灌制度擬定。

圖3 和平灌區(qū)磷污染模型的基礎(chǔ)資料及子流域劃分

3.2 模型驗(yàn)證

3.2.1 野外監(jiān)測試驗(yàn)

本研究分別在2017和2018年5—9月在和平灌區(qū)開展了野外監(jiān)測試驗(yàn)。

試驗(yàn)田：根據(jù)典型性、代表性、數(shù)據(jù)可得性和便利性原則選擇與灌區(qū)灌溉模式相似的常規(guī)灌溉典型試驗(yàn)田，田塊大小為100 m×50 m，位于慶安縣水稻灌溉試驗(yàn)站。在排水口處安裝三角堰，三角堰前安裝路格自動(dòng)水位計(jì)記錄水位變化，計(jì)數(shù)間隔10 min，以測量排水量；田間安裝2處自動(dòng)水位計(jì)，記錄田面水深；田塊灌溉水量采用進(jìn)水管的水表計(jì)數(shù)。試驗(yàn)田每5 d取水樣1次，排水期適當(dāng)加密，人工取樣200 mL，原狀水樣送實(shí)驗(yàn)室分析，測定磷濃度。利用DGT（Diffusive Gradients in Thin-films）和HR-Peeper裝置原位測量垂向土壤和孔隙水磷和鐵含量。監(jiān)測數(shù)據(jù)主要用于試驗(yàn)田尺度模型水文和水質(zhì)驗(yàn)證。

排水支溝：灌區(qū)區(qū)域尺度上選擇2條控制范圍明確的排水支溝，分別為一排水和七排水，具體位置見圖2。在排水支溝末端布設(shè)超聲波多普勒流量計(jì)監(jiān)測其排水量，監(jiān)測時(shí)間間隔為1 h，并采用TC-8000D自動(dòng)取樣器（青島溯源環(huán)保設(shè)備有限公司）定期采取水樣。采樣頻率為7 d 1次，有暴雨時(shí)加密采樣頻率。原狀水樣送實(shí)驗(yàn)室分析，測定磷濃度。監(jiān)測數(shù)據(jù)主要用于區(qū)域尺度模型水文和水質(zhì)驗(yàn)證。

3.2.2 試驗(yàn)田尺度模型驗(yàn)證

模型模擬了試驗(yàn)田生育期各層土壤含水率、田面積水深度和地表排水量。初期，灌溉用水應(yīng)用于干旱的試驗(yàn)田，形成田面積水。除強(qiáng)制排水期（穗形成期）外，積水深度隨降雨和灌溉而變化。為了促使模型穩(wěn)定，將2017年獲得的實(shí)測數(shù)據(jù)（有缺失）和2018年泡田期前的土壤含水率數(shù)據(jù)作為模型的初始數(shù)據(jù)進(jìn)行預(yù)熱，但未考慮冬春季凍融對(duì)土壤水分運(yùn)動(dòng)的影響。田塊飽和含水率和田間持水量分別取為55.9%和39.6%。計(jì)算步長為1 d，通過調(diào)整模型參數(shù)使模擬值與實(shí)測值吻合，得到模擬田間水深和田塊排水的變化過程，如圖4所示。根據(jù)圖4可以看出，6月5日—6月20日為水稻分蘗前期，需水量較大，一般不需要排水，初步確定三角堰安裝高度使田面水位高于10 cm時(shí)排水；6月20日水稻開始進(jìn)入分蘗中期和后期，這時(shí)水稻田需要控制水層深度保持在較低水位，進(jìn)而抑制無效分蘗，促進(jìn)抽穗開花，故6月20日之后調(diào)整三角堰高度為水位高于5 cm時(shí)排水，因此之后降雨較多時(shí)，水稻田排水量遠(yuǎn)大于分蘗前期的排水量?？傮w上，模型輸出符合實(shí)測資料，較好地預(yù)測了生育期稻田水分平衡狀況。

試驗(yàn)田積水層及排水的磷濃度模擬結(jié)果與實(shí)測結(jié)果對(duì)比如圖5所示。磷濃度變化主要受施肥影響，圖5a中5月和7月2個(gè)峰值是由于基肥和穗肥釋放到上覆水中所致。灌溉水的水質(zhì)對(duì)積水層的磷濃度也有一定影響，在曲線中形成一些小凸起。因田間排水是間斷性的，排水濃度的模擬值和實(shí)測值散點(diǎn)圖（圖5b）也顯示出了很好的相關(guān)性。

模型還模擬了田塊不同土壤層顆粒物和孔隙水中磷的濃度分布，圖6所示為分蘗期和拔節(jié)孕穗期土壤水可溶磷濃度分布結(jié)果。同時(shí)，將不考慮土壤微環(huán)境變化的模型模擬結(jié)果與本模型模擬、實(shí)測結(jié)果進(jìn)行了對(duì)比。兩模型的主要區(qū)別是不考慮微環(huán)境變化的模型中不考慮各層土壤溶解氧和pH值的差異，所以模型中磷在不同層土壤—水中的分配系數(shù)是相同的，且不隨時(shí)間發(fā)生變化。由圖6可知，不考慮土壤微環(huán)境變化模型對(duì)耕作層下部（無氧層）水中磷濃度是低估的，而對(duì)有氧層土壤水和上覆水中磷濃度是略微高估的。圖6中均方根誤差（RMSE）和決定系數(shù)（2）表明本模型給出了與原位實(shí)測結(jié)果更接近的結(jié)果，可見考慮土壤微環(huán)境變化及分層對(duì)磷污染模擬是必要的。

圖4 田面積水深度和排水流量驗(yàn)證

圖5 試驗(yàn)田積水層及排水的總磷濃度驗(yàn)證

3.2.3 灌區(qū)區(qū)域尺度模型驗(yàn)證

灌區(qū)典型排溝一排和七排的排水過程與實(shí)測排水量對(duì)比如圖7所示。由圖可以看出，一排、七排排水流量峰值較一致，分別是6月29日、7月13日、7月20日、以及7月25日，以上幾次峰值都是由于降雨較多產(chǎn)生的集中排水流量較大。一排和七排的排水總磷（TP）濃度與實(shí)測濃度對(duì)比如圖8所示。區(qū)域排水中TP濃度波動(dòng)較大，整體呈波動(dòng)減小趨勢(shì)，第1次峰值是由于基肥施用后首次由于降雨造成田塊排水形成的徑流攜帶了較多的磷，第2次峰值是6月下旬早施穗肥后降雨產(chǎn)生的排水?dāng)y帶了一定量的磷，經(jīng)過幾次排水的沖刷，磷流失較多，后期沒有施用磷肥，因此排水中的磷含量平穩(wěn)的減小。

注：SMS表示土壤微環(huán)境分層。

圖7 典型支溝排水流量驗(yàn)證

圖8 典型支溝排水總磷濃度驗(yàn)證

3.2.4 模型參數(shù)率定和模擬效果評(píng)價(jià)

經(jīng)過參數(shù)敏感性分析，磷污染模型主要對(duì)表3中所列參數(shù)較為敏感。為取得最好的模型驗(yàn)證結(jié)果，率定的參數(shù)取值如表3所示。采用Nash-Suttcliffe效率系數(shù)（NSE）、RMSE和2作為模型模擬效果的評(píng)價(jià)指標(biāo)，可得模型對(duì)灌區(qū)的水文水質(zhì)模擬效果見表4?？梢姡哪K的NSE、2均大于0.82，水質(zhì)模塊的NSE、2均大于0.81。模型性能總體上在預(yù)期之內(nèi)，考慮到流域面積和土地使用的復(fù)雜性，模型效率大于0.8被認(rèn)為是可以接受的。灌區(qū)磷污染物模型模擬結(jié)果可以反映水稻生育期間灌區(qū)水分運(yùn)動(dòng)和污染物遷移過程，模擬結(jié)果可靠。

表3 模型參數(shù)率定結(jié)果

表4 模型水文水質(zhì)模擬效果評(píng)價(jià)

3.3 灌區(qū)磷污染負(fù)荷模擬統(tǒng)計(jì)

灌區(qū)2018年5—9月水稻田的磷產(chǎn)出負(fù)荷根據(jù)區(qū)域內(nèi)水稻田的排水量和磷濃度過程計(jì)算，磷滲漏輸出負(fù)荷根據(jù)各子流域滲漏量和滲漏磷濃度計(jì)算得到，并據(jù)此分析各個(gè)生育期的磷輸入輸出量，如表5所示。分蘗期和拔節(jié)孕穗期徑流輸出負(fù)荷分別為0.85 kg/hm2和0.60 kg/hm2。分蘗期由于降雨導(dǎo)致磷輸出負(fù)荷較大，拔節(jié)抽穗期雖排水較多，但是經(jīng)過分蘗期排水的沖刷，排水磷濃度較小，輸出負(fù)荷較分蘗期小。其他生育期無排水。但泡田期和分蘗期滲漏輸出負(fù)荷最大，分別為0.11 kg/hm2和0.16 kg/hm2。主要是由于基肥和早施穗肥導(dǎo)致土壤水濃度較高所致。水稻生育期內(nèi)通過排水和滲漏流失的磷為1.88 kg/hm2，約占施肥和灌溉磷輸入的5.7%。各主要排水支溝控制農(nóng)田地表排水負(fù)荷、滲漏輸出負(fù)荷及排溝末端輸出負(fù)荷如圖9所示。由圖可以看出灌區(qū)一排和四排輸出的磷負(fù)荷總量最大，分別為1.40 t和1.39 t。這是因?yàn)橐慌藕退呐趴刂泼娣e較大，水田和居民點(diǎn)較多，污染物輸出大于其余排水支溝。

表5 水稻生育期內(nèi)磷的輸入輸出平衡

圖9 和平灌區(qū)2018年5—9月各排水支溝及控制區(qū)域磷負(fù)荷

4 結(jié) 論

1）本文構(gòu)建了適合平原灌區(qū)的非點(diǎn)源磷污染機(jī)理模型，該模型考慮水稻田實(shí)際水文循環(huán)過程和基于土壤微環(huán)境分層的磷運(yùn)移過程，建立的水分運(yùn)動(dòng)方程和磷遷移轉(zhuǎn)化方程比較符合平原灌區(qū)實(shí)際情況，能準(zhǔn)確模擬平原水稻灌區(qū)的水量和水質(zhì)變化。

2）利用黑龍江省和平灌區(qū)2018年試驗(yàn)田土壤水、積水層排水和典型支溝的排水流量、總磷濃度過程對(duì)模型進(jìn)行率定與驗(yàn)證，試驗(yàn)田積水水深、排水、一排和七排徑流量模擬效果指標(biāo)Nash-Sutcliffe效率系數(shù)分別為0.967、0.980、0.820和0.861，決定系數(shù)分別為0.976、0.981、0.815和0.877；試驗(yàn)田、一排和七排排水磷濃度模擬效果指標(biāo)Nash-Sutcliffe效率系數(shù)分別為0.811、0.890、0.826，決定系數(shù)分別為0.821、0.930、0.910，所有指標(biāo)均大于0.81，模擬效果較好，表明該模型適用于以水稻為主要作物的平原灌區(qū)。

3）根據(jù)模擬的灌區(qū)水文水質(zhì)過程計(jì)算了整個(gè)生育期稻田輸入輸出量及各排水支溝輸出負(fù)荷，分蘗期和拔節(jié)孕穗期徑流輸出負(fù)荷分別為0.85、0.60 kg/hm2，泡田期和分蘗期滲漏輸出磷負(fù)荷最大，分別為0.11、0.16 kg/hm2。水稻生育期內(nèi)通過排水和滲漏流失的磷約占施肥和灌溉磷輸入的5.7%。其中，灌區(qū)一排和四排輸出的磷負(fù)荷總量最大。

[1] 李懷恩，沈晉，劉玉生. 流域非點(diǎn)源污染模型的建立與應(yīng)用實(shí)例[J]. 環(huán)境科學(xué)學(xué)報(bào)，1997，17(2)：12-18. Li Huai'en, Shen Jin, Liu Yusheng. The establishment and case study of the model for nonpoint source pollution for watershed[J]. Acta Scientiae Circumstantiae, 1997, 17(2): 12-18. (in Chinese with English abstract)

[2] 陶?qǐng)@，王少麗，管孝艷，等. 青海省農(nóng)業(yè)面源污染源特征分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019，35(10)：164-172. Tao Yuan, Wang Shaoli, Guan Xiaoyan, et al. Characteristic analysis of non-point source pollution in Qinghai Province[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(10): 164-172. (in Chinese with English abstract)

[3] Neitsch S L, Arnold J G, Kiniry J R, et al. Soil and water assessment tool theoretical documentation, Version 2000: TWRI Report TR-191[R]. Texas: Texas Water Resources Institute, College Station, 2002.

[4] Wu D, Cui Y, Xie X, et al. Improvement and testing of SWAT for multi-source irrigation systems with paddy rice[J]. Journal of Hydrology, 2019, 568: 1031-1041.

[5] Ouyang Y, Higman J, Hatten J. Estimation of dynamic load of mercury in a river with BASINS-HSPF model[J]. Journal of Soils and Sediments, 2012, 12(2): 207-216.

[6] 耿潤哲，張鵬飛，龐樹江，等. 不同氣候模式對(duì)密云水庫流域非點(diǎn)源污染負(fù)荷的影響[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2015，31(22)：240-249. Geng Runzhe, Zhang Pengfei, Pang Shujiang, et al. Impact of different climate change scenarios on non-point source pollution losses in Miyun Reservoir[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(22): 240-249. (in Chinese with English abstract)

[7] Arnold J G, Williams J R, Nicks A D, et al. SWRRB: A Basin Scale Simulation Model for Soil and Water Resources Management[M]. College Station: Texas A&M University Press, 1990.

[8] Perrone J. Hydrologic Modeling of an Agricultural Watershed in Quebec Using AGNPS[D]. Montreal: McGill University, 1997.

[9] 鐘科元，陳興偉，陳瑩，等. 地形和土壤/土地利用空間聚合對(duì)徑流和輸沙量影響的AnnAGNPS模擬[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2016，32(8)：127-135. Zhong Keyuan, Chen Xingwei, Chen Ying, et al. Simulation of effects of topography and soil/land use spatial aggregation on sediment yield and runoff using AnnAGNPS[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(8): 127-135. (in Chinese with English abstract)

[10] Bouraoui F, Dillaha T A. ANSWERS-2000: Non- point-source nutrient planning model[J]. Journal of Environmental Engineering, 2000, 126(11): 1045- 1055.

[11] Ahmadi S H, Amin S, Keshavarzi A, et al. Simulating watershed outlet sediment concentration using the ANSWERS model by applying two sediment transport capacity equations[J]. Biosystems Engineering, 2006, 94(4): 615-626.

[12] 仕玉治，張弛，周惠成，等. SWAT模型在水稻灌區(qū)的改進(jìn)及應(yīng)用研究[J]. 水電能源科學(xué)，2010，28(7)：18-22. Shi Yuzhi, Zhang Chi, Zhou Huicheng, et al. Development and application of SWAT Model to paddy district in watershed scale[J]. Water Resources and Power, 2010, 28(7): 18-22. (in Chinese with English abstract)

[13] 鄭捷，李光永，韓振中，等. 改進(jìn)的SWAT模型在平原灌區(qū)的應(yīng)用[J]. 水利學(xué)報(bào)，2011，42(1)：88-97. Zheng Jie, Li Guangyong, Han Zhenzhong, et al. Application of modified SWAT model in plain irrigation district[J]. Journal of Hydraulic Engineering, 2011, 42(1): 88-97. (in Chinese with English abstract)

[14] 胡文慧，李光永，郭亞潔，等. 汾河灌區(qū)農(nóng)業(yè)面源污染經(jīng)驗(yàn)統(tǒng)計(jì)模型的構(gòu)建與驗(yàn)證[J]. 中國農(nóng)業(yè)大學(xué)學(xué)報(bào)，2015，20(2)：207-215. Hu Wenhui, Li Guangyong, Guo Yajie, et al. Construction and verification of an empirical model for agricultural nonpoint source pollution in the Fenhe irrigation district[J]. Journal of China Agricultural University, 2015, 20(2): 207-215. (in Chinese with English abstract)

[15] 徐保利，代俊峰，俞陳文炅，等. 漓江流域氮磷排放對(duì)水肥管理和下墊面屬性變化的響應(yīng)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(2)：245-254. Xu Baoli, Dai Junfeng, Yuchen Wenjiong, et al. Responses of nitrogen and phosphorus emissions to water and fertilizer management and underlying surface property changes in Lijiang River Basin[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(2): 245-254. (in Chinese with English abstract)

[16] 郝守寧，普布次仁，董飛. 林芝畜禽養(yǎng)殖糞便排放時(shí)空演變及耕地污染負(fù)荷分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019，35(16)：225-232. Hao Shouning, Pubu Ciren, Dong Fei. Analysis on spatial-temporal evolution of livestock manure emission and pollution load of cultivated land in Nyingchi[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(16): 225-232. (in Chinese with English abstract)

[17] 李燕. 基于HSPF模型的水文水質(zhì)過程模擬研究--以中田河流域?yàn)槔齕D]. 南京：南京農(nóng)業(yè)大學(xué)，2013. Li Yan. Study on Hydrology and Water Quality Process Simulation Based on HSPF Model: A Case Study in Zhongtianhe River Watershed[D]. Nanjing: Nanjing Agricultural University, 2013. (in Chinese with English abstract)

[18] Yan R, Gao J, Huang J. WALRUS-paddy model for simulating the hydrological processes of lowland polders with paddy fields and pumping stations[J]. Agricultural Water Management, 2016, 169: 148-161.

[19] Djaman K, Rudnick D, Mel V C, et al. Evaluation of valiantzas’ simplified forms of the FAO-56 Penman- Monteith reference evapotranspiration model in a humid climate[J]. Journal of Irrigation and Drainage Engineering, 2017, 143(8): 06017005.

[20] 楊霞，邵東國，徐保利. 東北寒區(qū)黑土稻田土壤水分剖面二維運(yùn)動(dòng)規(guī)律研究[J]. 水利學(xué)報(bào)，2018，49(8)：1017-1026. Yang Xia, Shao Dongguo, Xu Baoli. Indoor experiment and simulation of soil water two-dimensional movement of the paddy fields in the northeast frigid of China[J]. Journal of Hydraulic Engineering, 2018, 49(8): 1017-1026. (in Chinese with English abstract)

[21] van Genuchtem M T H. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(5): 892-898.

[22] Borah D K. Sediment Discharge Model for Small Watersheds[J]. Transactions of the ASAE, 1989, 32(3): 874-880.

[23] Di Toro D M. Sediment Flux Modeling[M]. New York: Wiley-Interscience, 2001.

[24] Roelcke M, Han Y, Cai Z, et al. Nitrogen mineralization in paddy soils of the Chinese Taihu Region under aerobic conditions[J]. Nutrient Cycling in Agroecosystems, 2002, 63(2/3): 255-266.

[25] Huang L, Fang H, Reible D. Mathematical model for interactions and transport of phosphorus and sediment in the Three Gorges Reservoir[J]. Water Research, 2015, 85: 393-403.

[26] Stumm W, Huang C P, Jenkins S R. Specific chemical interaction affecting the stability of dispersed system[J]. Croatica Chemica Acta, 1970, 42(2): 223-245.

[27] 湯鴻霄，錢易，文湘華，等. 水體顆粒物和難降解有機(jī)物的特性與控制技術(shù)原理：上卷水體顆粒物[M]. 北京：中國環(huán)境科學(xué)出版社，2000.

[28] Zhang J, Zhou X, Chen L, et al. Comparison of the abundance and community structure of ammonia oxidizing prokaryotes in rice rhizosphere under three different irrigation cultivation modes[J]. World Journal of Microbiology & Biotechnology, 2016, 32: 85.

[29] 劉錦濤，黃萬勇，楊士紅，等. 加氣灌溉模式下稻田土壤水溶解氧的變化規(guī)律[J]. 江蘇農(nóng)業(yè)科學(xué)，2015，43(2)：389-392. Liu Jintao, Huang Wanyong, Yang Shihong, et al. The Change law of dissolved oxygen in paddy soil under aerated irrigation mode[J]. Jiangsu Agricultural Sciences, 2015, 43(2): 389-392. (in Chinese with English abstract)

[30] Berlin M, Suresh Kumar G, Nambi I M. Numerical modelling on transport of nitrogen from wastewater and fertilizer applied on paddy fields[J]. Ecological Modelling. 2014, 278: 85-99.

[31] 董志勇. 環(huán)境流體力學(xué)[M]. 北京：科學(xué)出版社，2015.

[32] Pu J, Wang S, Ni Z, et al. Implications of phosphorus partitioning at the suspended particle-water interface for lake eutrophication in China’s largest freshwater lake, Poyang Lake[J]. Chemosphere, 2021, 263: 128334.

[33] 李碩，賴正清，王橋，等. 基于SWAT模型的平原河網(wǎng)區(qū)水文過程分布式模擬[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2013，29(6)：106-112. Li Shuo, Lai Zhengqing, Wang Qiao, et al. Distributed simulation for hydrological process in Plain River network region using SWAT model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(6): 106-112. (in Chinese with English abstract)

Phosphorus pollution model for plain paddy irrigation district based on soil microenvironment sub-stratification

Chen Minghong, Bai Lu, Liu Xuanye, Guo Jiansong, Wang Chaozi

(,,100083,)

Non-point source phosphorus pollution generated from irrigated farmlands is one of the main causes of local and regional eutrophication. However, current phosphorus pollution models either do not include the water movement in irrigation and drainage process or do not consider the phosphorus transformation under the exchanging aerobic and anaerobic conditions. Therefore, we developed a physically based phosphorus pollution model to quantitatively describe the water movement and phosphorus fate and transport processes in irrigated paddy fields in plain areas. The simulation of the runoff yield in an irrigation area was based on the water balance equations describing the water input and output of the paddy fields and the motion wave equations describing the water movement in the drainage channel networks. The simulation of the excess phosphorus yield was based on the convection diffusion equations and a phosphorus transformation model considering the soil sub-stratification-the cultivated horizon was sub-divided into aerobic and anaerobic layers. In this way, the changes in dissolved oxygen and the processes of phosphorus transformation in different soil layers caused by the alternating wet and dry conditions could be quantitatively described in details. The phosphorus flux of diffusion, particle mixing and infiltration between the water layer, the aerobic soil layer and the anaerobic soil layer were also quantified. The model was calibrated and verified with the observed ponding water depth, drainage discharge, and phosphorus concentrations in the runoff and soil water in one experimental paddy field and two typical drainage ditches in Heping Irrigation District, Heilongjiang, China in 2018. The simulated drainage discharge and phosphorus concentrations of the experimental paddy field and the drain ditches agreed well with the observations. The Nash-Sutcliffe efficiency coefficient (NSE) and coefficient of determination (2) of the simulated drainage discharge were greater than 0.820 and 0.815, respectively. And the NSE and2of simulated total phosphorus concentration were greater than 0.811 and 0.821, respectively. The simulated vertical distribution of the soil soluble phosphorus obtained by considering the aerobic and anaerobic layers of the cultivated horizon were closer to the in situ observation than the results obtained with the same model but do not consider the soil sub-stratification. Then, the verified phosphorus pollution model was used to estimate the non-point source phosphorus pollution in the whole Heping Irrigation District. The phosphorus loss through drainage and leakage during the growth stages of rice was 1.88 kg/hm2, which was about 5.7% of the phosphorus input from fertilization and irrigation. Among the 1.88 kg/hm2phosphorus loss, the phosphorus output load of runoff at the tillering stage (0.85 kg/hm2), and the jointing and booting stage (0.60 kg/hm2) was the first and second largest loss, due to rainfall washout of soil phosphorus. The loss by leakage output load was the second and first largest at the soaking stage (0.11 kg/hm2) and the tillering stage (0.16 kg/hm2), due to the basic fertilizer and the early booting stage fertilizer. For the whole Heping Irrigation District, the total excess phosphorus exported from the first ditch (1.40 t) and the fourth ditch (1.39 t) were the first and second largest, due to their larger control area of the irrigation district. Overall, the physically based phosphorus pollution model developed in this study included the water movement in irrigation and drainage process, considered the phosphorus transformation under the exchanging aerobic and anaerobic conditions caused by the alternating wet and dry conditions, and provided more accurate estimation of phosphorus fate and transport in irrigated paddy fields in plain areas.

water pollution; drainage; redox reactions; microenvironment sub-stratification; non-point source phosphorus pollution; migration and transformation; irrigation district

陳明洪，白璐，劉宣冶，等. 基于土壤微環(huán)境分層的平原水稻灌區(qū)磷污染模型[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(20)：17-26.doi：10.11975/j.issn.1002-6819.2020.20.003 http://www.tcsae.org

Chen Minghong, Bai Lu, Liu Xuanye, et al. Phosphorus pollution model for plain paddy irrigation district based on soil microenvironment sub-stratification[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(20): 17-26. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.20.003 http://www.tcsae.org

2019-05-05

2020-10-10

國家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（2016YFC0400107）；國家自然科學(xué)基金面上項(xiàng)目（52079137）

陳明洪，副教授，主要從事流域水環(huán)境模擬與調(diào)度研究。Email：chenminghong@cau.edu.cn

10.11975/j.issn.1002-6819.2020.20.003

X52

A

1002-6819(2020)-20-0017-10