基于土壤微环境分层的平原水稻灌区磷污染模型

陈明洪，白 璐，刘宣冶，郭建松，汪超子

基于土壤微环境分层的平原水稻灌区磷污染模型

陈明洪，白 璐，刘宣冶，郭建松，汪超子

（中国农业大学水利与土木工程学院，北京 100083）

灌溉农田产生的非点源磷是造成水体富营养化的主要原因之一，但是目前国内外的磷污染模型对于平原灌区灌溉和排水管理下的水分运动过程和氧化还原条件下的磷转化过程的定量表征还比较欠缺。该研究构建了适用于平原水稻灌区水分运动和磷转化迁移的机理性磷污染模型，模型根据稻田水量平衡和沟道运动波方程模拟灌区产汇流，采用考虑土壤微环境分层的磷转化模型和对流扩散方程模拟灌区产汇污。模型中将耕作层分为有氧层和无氧层，定量表征了由于水田干湿交替导致的土壤分层溶解氧变化和磷的转化过程。为了验证模型的可靠性，利用黑龙江省和平灌区2018年试验田实测田间土壤水、积水及排水和2条支沟实测排水的水量和水质数据对模型进行了率定和验证。验证结果显示试验田、一排和七排排水的径流流量、磷浓度的模拟结果与实测结果都吻合较好。模拟排水流量的Nash-Sutcliffe效率系数（NSE）和决定系数（2）分别大于0.820和0.815；模拟总磷浓度的NSE和2分别大于0.811和0.821；考虑土壤微环境分层后得到的土壤可溶磷垂向分布结果比不考虑分层时与原位实测结果更接近。该磷污染模型被用于模拟和平水稻灌区的非点源磷污染。灌区磷污染浓度过程统计结果显示，2018年水稻生育期内通过排水和渗漏流失的磷为1.88 kg/hm2，约占施肥和灌溉磷输入的5.7%。其中分蘖期和拔节孕穗期径流磷输出负荷最大，分别为0.85和0.60 kg/hm2；泡田期和分蘖期渗漏输出负荷最大，分别为0.11和0.16 kg/hm2。灌区一排和四排由于控制面积大，输出磷污染物总量大于其他排水沟。该研究所构建的磷污染模型包含稻田灌溉和排水过程的水分运动及稻田干湿交替引起的氧化还原变化条件下的磷转化和运移，可为平原灌区水肥运筹管理下的磷运移模拟提供更准确的方法。

水体污染；排水；氧化还原反应；微环境分层；非点源磷污染；迁移转化；灌区

0 引 言

水体污染已经成为当今社会人类面临的主要环境问题之一，其中磷污染是造成水体富营养化的主要限制因子，为解决这一难题，各国都在控制营养物质的来源上投入了大量的人力财力。大量的研究结果表明[1-2]，农业非点源已经成为当今水环境的主要污染源。但农业非点源污染具有形成过程随机性大、影响因子复杂、分布范围广、危害大且产生机理模糊等特点，由此导致模拟和控制非点源污染难度较大。

定量化表征非点源污染负荷及过程是水环境污染治理的重要基础工作。目前，模拟具有一定坡度的流域非点源污染模型比较完善，出现了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]等一系列模型，但针对平原灌区的非点源污染模型较少，大多数都是基于以上流域模型进行改进。仕玉治等[12]根据水量平衡方程和水稻田的灌溉方式修改了SWAT模型的自动灌溉技术，模拟努敏河流域长阁灌区径流量的效果较好；郑捷等[13]在考虑灌区实际渠道的基础上对子流域划分以及作物耗水量计算模块等方面进行了改进，成功地应用于汾河灌区的水量平衡模拟；胡文慧等[14]根据汾河灌区实地的灌溉施肥调查资料以及实测的总氮、总磷、水文站监测数据构建SWAT模型，模拟了汾河灌区径流和污染物负荷。SWAT模型虽然已有应用于平原灌区的案例，但大多都采用模拟自然流域的SCS（Soil Conservation Service）径流曲线法模拟产流，或者仅考虑农田灌溉的人为影响调整灌溉模块，并未考虑灌区作物实际排水模式，比如现行节水政策下水稻田的控制灌溉或干湿交替灌溉等。所以上述模型在平原水稻灌区人为管理下的水分运动过程模拟中仍显能力不足，更影响了水文驱动条件下的污染物迁移过程。

此外，对于磷污染物迁移模块，一般非点源模型只考虑其质量守恒和浓度平衡[15-16]，计算其平均输出负荷，并没有考虑磷的动力迁移过程。然而水稻生育期内干湿交替引起的土壤氧化还原条件变化较大，土壤磷素在有氧和无氧变化环境下的生化反应过程差异较大，在上述模型中均未得到体现。因此，需要建立适合平原水稻灌区水分运动及磷污染过程机理的模型，以准确定量地模拟磷在水稻生育期内的迁移转化过程，从而为评价平原水稻灌区磷污染对流域水环境的影响打下基础。本研究建立基于土壤微环境分层的非点源磷污染模型，定量表征磷污染物在水土、土壤有氧和无氧分层界面的迁移转化规律，并将模型应用于黑龙江和平灌区非点源磷污染过程的模拟。

1 灌区磷污染模型原理及构建

1.1 水文模块

水分运动为污染物提供迁移的介质和能量，水文路径同时也是污染物迁移的路径，合理并准确描述水文过程直接影响整个模型模拟结果[17]。与SWAT、HSPF模型相似，根据灌区地形、土壤和水文特征将区域划分为若干个水文响应单元，在每个单元内根据天然降雨和灌溉制度模拟地表水和土壤水的运动过程。

水稻田可能存在有和无淹没水层2种情况。当田面有淹没水层且有排水时，田块地表水流运动可用一维水动力学水流运动控制方程[18]描述

1）田间灌溉

2）田间地表排水

当田面有水层但无排水或无水层时，稻田土壤水运动过程均可简化为二维垂向土壤水动力学问题，模型利用Richards方程[20]计算

式中为土壤压力水头，cm；为垂向距离，向下为正，cm；()为土壤非饱和水力传导度，cm/min；()为比水容量，cm。选用水流和溶质运移模拟软件（HYDRUS）提供的van Genuchten-Mualem公式[21]计算土壤水力特性参数。

由于水稻田犁底层具有较强的阻水作用，本文暂不考虑犁底层向底土层和地下水的垂直下渗作用，只考虑犁底层的侧渗作用。在耕作层仅考虑向犁底层的垂直入渗，侧面为零通量边界。田面存在淹没水层时，()=0，()=s，s为土壤饱和水力传导度，cm/min。无淹没水层时，采用蒸发强度作为上边界。

稻田地表排水根据排水沟道汇流顺序汇入各级排水沟道，排水在沟道水流演进采用略去惯性项和压力项的一维圣维南方程[22]模拟。

1.2 磷的迁移转化过程

水稻在不同生育期需水状况不同，导致水田呈现干湿交替的水土环境，从而影响稻田的磷转化运移过程。干湿交替会造成稻田土壤水土界面有氧-无氧环境频繁变化，而磷转化的生物化学过程又受到有氧-无氧环境的强烈影响，因此需要考虑土壤和水层、有氧层与无氧层之间的界面效应。根据土壤中含氧量差异将稻田土壤分为有氧层和无氧层。一般河道沉积物有氧层为地表以下3 mm厚[23]，但水稻土壤有氧层可以取1 cm[24]，无氧层取有氧层以下30 cm厚。在田面积水和无人为干扰的情况下，有氧层和无氧层的厚度保持不变，但田面无淹没水层时，有氧层会随着土壤含水层下降而下降。本研究提出的稻田土壤磷转化过程见图1。

注：R1为施肥磷通量；R2、R3分别为有氧层和上覆水、有氧层和无氧层的磷扩散通量；R4为颗粒混掺磷通量；R5为稻田磷衰减通量；R6、R7分别为水体层向有氧层、有氧层向无氧层的磷渗透通量；R8为灌溉水携带磷通量。

1.2.1 磷在稻田的循环过程

由图1可以看出，磷在稻田的循环由外界输入、各层之间的交换、水层的生化衰减反应3部分构成。磷在各层之间的交换主要包括3个过程：水体之间的扩散作用、生物扰动造成的颗粒混掺、渗透作用。

1）磷的来源主要是施肥、降雨灌溉外界带入，施入土壤中的有机磷肥或动植物残体的有机磷在微生物的作用下会产生矿化分解变成无机磷，主要是正磷酸盐。施肥和灌溉带入的磷通量分别用1、8表示。

2）水体之间的扩散作用。用2、3分别表示有氧层和上覆水、有氧层和无氧层之间的磷扩散通量。扩散是在浓度梯度作用下发生的，扩散通量与浓度梯度成正比，表达式如下[23]

式中C为水中污染物浓度，mg/L；C为有氧层中污染物浓度，mg/L；C为无氧层中污染物浓度，mg/L；L01、L12分别是水体层与有氧层、有氧层与无氧层间的水体扩散传质系数，m/s；d1、d2分别表示溶解磷在有氧层和无氧层中所占的比例。

3）有氧层和无氧层之间颗粒混掺作用。在生物扰动作用下，有氧层和无氧层间的土壤颗粒发生混掺，引起磷交换。颗粒混掺交换通量用4表示[23]

式中12表示土壤颗粒的混掺速率，m/s；p1、p2分别表示吸附磷在有氧层和无氧层中所占的比例。

4）渗透作用。由于重力作用，水体层和有氧层、有氧层和无氧层之间会发生渗透，渗透水携带着磷向下迁移。有氧层和无氧层的渗透通量分别用6、7表示[23]

式中L01和L02分别是水体层与有氧层、有氧层与无氧层之间的渗透作用传质系数，m/s。

5）衰减作用。稻田水层中的磷由于自净作用会随着时间衰减，满足一阶反应动力学。衰减通量用5表示

式中f为衰减系数，d-1。

1.2.2 考虑土壤微环境分层的磷分配过程

磷在固液相的分配比例是模型的重要参数，磷的吸附解吸过程是分配比例确定的依据。土壤颗粒中的铁（氢）氧化物表面存在着大量表面羟基，对无机磷的吸附作用很强，铁及其氧化物的存在形式决定了磷在土壤中的形态[25]。在有氧条件下，溶解态的Fe2+被氧化为FeOOH(s)，对磷的吸附作用较强；在无氧或相对缺氧的条件下，FeOOH(s)被还原为Fe2+，对磷的吸附作用减弱。可见溶解磷和吸附磷在土壤中所占的比例dj和pj是溶氧浓度的函数。

有氧条件

无氧条件

分配系数的最早提出是一个经验常数，缺少基本的理论推导。自20世纪70年代初期Stumm等[26]提出表面络合模型后，其被广泛应用于吸附过程的研究。表面络合模型是基于颗粒表面电荷，用溶液中络合平衡类似的方法处理颗粒物界面上的结合问题。它根据表面质子反应、与磷酸根的配体交换确定2个表面酸度常数和3个表面络合常数，再根据表面络合模型确定吸附态和溶解态的磷含量[25,27]，进而可确定溶解相和颗粒相的比例dj

稻田表层土壤及水稻根际土壤有相对较高的溶解氧含量，而非根际土壤和深层土壤溶解氧含量极低[28]。刘锦涛等[29]研究稻田在常规灌溉和控制灌溉条件下土壤溶解氧的分布，发现溶解氧浓度在不同生育期和不同土壤深度均有所差异。因此，本文根据不同时期农田水层深度及气象资料计算表层水和土壤有氧层的溶氧浓度变化过程，依据土壤初始有效铁含量确定各层土壤中FeOOH(s)的含量，使用络合模型再根据式（12）确定磷在固液相的分配比例。表层水/土壤层氧气浓度计算如下[30]

1.2.3 磷迁移转化整体控制方程

上节内容即为水稻田磷迁移转化的过程及其定量表征，包括田面积水、有氧层和无氧层，由此得到磷的总控制方程，表达式如下

式中1和2分别为土壤有氧层和无氧层厚度，m；D为无氧层污染物扩散系数，m2/s；为土壤层垂向坐标，m。

污染物进入沟道后应用一维非保守性污染物迁移扩散方程[31]模拟，基本方程如下

式中d为沟道水深，m；为沟道污染物浓度，mg/L；为沟道流速，m/s；为沟道污染物扩散系数，m2/s；d为沟道综合衰减系数，d-1；为沟道水平方向坐标，m。

水稻田排水及沟道水体中如存在悬浮颗粒物，应考虑颗粒物对磷的携带作用。一般采用分配系数K表征磷在颗粒物和水相的分配。水体总磷浓度计算如下[32]

式中C为水体总磷浓度，mg/L；为固体颗粒物浓度，kg/L；K为磷素的固液分配系数，L/kg。本研究在野外试验时采集了试验田及主要排水沟的水样，实测悬浮固体颗粒物浓度极低，故后文模拟中忽略水体中固体含磷量。

本文所列偏微分方程均采用隐式有限差分方法进行离散，并给出初始条件和边界条件求解计算域内所有变量。

2 研究区概况

和平灌区位于呼兰河流域的呼兰河冲积平原上，灌区大部分面积属于庆安县，少部分属于铁力市，灌区总面积约为6 700 hm2。灌区位置和主要水系如图2所示，和平灌区和劳模灌区分别位于呼兰河两岸。和平灌区地理坐标为127°17'～127°49'E，46°34'～47°07'N。庆安县属于中温带半干旱半湿润的大陆性季风气候，灌区多年平均降水量545.3 mm，变化范围在450～700 mm之间，降水年际年内分配不均，空间分布趋势是由东北向西南递减。该区蒸发量及其时空分布变化较大，多年平均水面蒸发量为664.5 mm。

和平灌区的种植作物95%以上为水稻，灌区水稻生育期为5—9月，期间施基肥、分蘖肥和穗肥，具体情况见表1。该区为自流引水灌区，灌区从呼兰河及安邦河、拉林清河等支流引水，灌区内有1条干渠，18条支渠，10条排水支沟，但9排和10排位于井灌区，生育期基本无排水，本文模拟未考虑。

注：为保持汇水流域完整性，除和平灌区外，图中流域还包括了呼兰河西北侧的劳模灌区。

表1 和平灌区水稻生育期及施肥制度

为了解灌区土壤背景属性及其空间异质性分布，在整个灌区中选取了44个典型采样点，分别取表层0～20 cm土壤样品回实验室进行测定。土壤质地和背景资料见表2。实测数据显示，研究区土壤机械组成差异较小，根据美国制分级标准，灌区土壤基本为粉壤土。土壤总磷平均质量分数约为0.80～1.08 g/kg，东北部总磷背景含量较高，整体趋势向西南递减。土壤有效铁的空间差异性较大。

表2 和平灌区农田土壤的基本性质

3 模型应用

3.1 模型构建与数据来源

灌区磷污染物模型采用的基础数字高程模型（DEM）是SRTM 90 m精度的数据（图3a）。借助GIS平台根据DEM高程图和人工burn-in[33]功能生成符合灌区实际渠系情况的河网，将灌区调研的实际河网叠加在已有的DEM上，保持河道所在栅格的高程值不变，而将其他非河道所在位置的栅格整体增加一个微小值，加大河道所在栅格和周边栅格之间的坡度，使得河道所在栅格的汇水能力增强，提高河道提取的准确程度。土地利用采用2010年黑龙江省遥感土地利用数据（图3b）。土壤类型分布采用南京土壤所1∶100万土壤空间数据，并利用灌区采样数据进行修正（图3c）。将模型基础数据导入ARCSWAT软件划分水文响应单元和子流域，本流域共划分出64个子流域（图3d）。为保持汇水流域完整性，除和平灌区外，DEM及流域划分还包括了呼兰河西北侧的劳模灌区（见图2和图3），但本研究数据分析仅含和平灌区模拟数据。气象数据来自黑龙江省水科院提供的庆安县气象站资料，数据系列时间长度2000—2018 年，数据主要包括逐日的观测数据，包括最高气温、最低气温、降水、相对湿度、风速和日照时数等。此外，和平灌区试验田气象数据采用试验站小型气象站的2015—2018年数据。灌区各子流域灌溉数据利用和平灌区灌溉定额和轮灌制度拟定。

图3 和平灌区磷污染模型的基础资料及子流域划分

3.2 模型验证

3.2.1 野外监测试验

本研究分别在2017和2018年5—9月在和平灌区开展了野外监测试验。

试验田：根据典型性、代表性、数据可得性和便利性原则选择与灌区灌溉模式相似的常规灌溉典型试验田，田块大小为100 m×50 m，位于庆安县水稻灌溉试验站。在排水口处安装三角堰，三角堰前安装路格自动水位计记录水位变化，计数间隔10 min，以测量排水量；田间安装2处自动水位计，记录田面水深；田块灌溉水量采用进水管的水表计数。试验田每5 d取水样1次，排水期适当加密，人工取样200 mL，原状水样送实验室分析，测定磷浓度。利用DGT（Diffusive Gradients in Thin-films）和HR-Peeper装置原位测量垂向土壤和孔隙水磷和铁含量。监测数据主要用于试验田尺度模型水文和水质验证。

排水支沟：灌区区域尺度上选择2条控制范围明确的排水支沟，分别为一排水和七排水，具体位置见图2。在排水支沟末端布设超声波多普勒流量计监测其排水量，监测时间间隔为1 h，并采用TC-8000D自动取样器（青岛溯源环保设备有限公司）定期采取水样。采样频率为7 d 1次，有暴雨时加密采样频率。原状水样送实验室分析，测定磷浓度。监测数据主要用于区域尺度模型水文和水质验证。

3.2.2 试验田尺度模型验证

模型模拟了试验田生育期各层土壤含水率、田面积水深度和地表排水量。初期，灌溉用水应用于干旱的试验田，形成田面积水。除强制排水期（穗形成期）外，积水深度随降雨和灌溉而变化。为了促使模型稳定，将2017年获得的实测数据（有缺失）和2018年泡田期前的土壤含水率数据作为模型的初始数据进行预热，但未考虑冬春季冻融对土壤水分运动的影响。田块饱和含水率和田间持水量分别取为55.9%和39.6%。计算步长为1 d，通过调整模型参数使模拟值与实测值吻合，得到模拟田间水深和田块排水的变化过程，如图4所示。根据图4可以看出，6月5日—6月20日为水稻分蘖前期，需水量较大，一般不需要排水，初步确定三角堰安装高度使田面水位高于10 cm时排水；6月20日水稻开始进入分蘖中期和后期，这时水稻田需要控制水层深度保持在较低水位，进而抑制无效分蘖，促进抽穗开花，故6月20日之后调整三角堰高度为水位高于5 cm时排水，因此之后降雨较多时，水稻田排水量远大于分蘖前期的排水量。总体上，模型输出符合实测资料，较好地预测了生育期稻田水分平衡状况。

试验田积水层及排水的磷浓度模拟结果与实测结果对比如图5所示。磷浓度变化主要受施肥影响，图5a中5月和7月2个峰值是由于基肥和穗肥释放到上覆水中所致。灌溉水的水质对积水层的磷浓度也有一定影响，在曲线中形成一些小凸起。因田间排水是间断性的，排水浓度的模拟值和实测值散点图（图5b）也显示出了很好的相关性。

模型还模拟了田块不同土壤层颗粒物和孔隙水中磷的浓度分布，图6所示为分蘖期和拔节孕穗期土壤水可溶磷浓度分布结果。同时，将不考虑土壤微环境变化的模型模拟结果与本模型模拟、实测结果进行了对比。两模型的主要区别是不考虑微环境变化的模型中不考虑各层土壤溶解氧和pH值的差异，所以模型中磷在不同层土壤—水中的分配系数是相同的，且不随时间发生变化。由图6可知，不考虑土壤微环境变化模型对耕作层下部（无氧层）水中磷浓度是低估的，而对有氧层土壤水和上覆水中磷浓度是略微高估的。图6中均方根误差（RMSE）和决定系数（2）表明本模型给出了与原位实测结果更接近的结果，可见考虑土壤微环境变化及分层对磷污染模拟是必要的。

图4 田面积水深度和排水流量验证

图5 试验田积水层及排水的总磷浓度验证

3.2.3 灌区区域尺度模型验证

灌区典型排沟一排和七排的排水过程与实测排水量对比如图7所示。由图可以看出，一排、七排排水流量峰值较一致，分别是6月29日、7月13日、7月20日、以及7月25日，以上几次峰值都是由于降雨较多产生的集中排水流量较大。一排和七排的排水总磷（TP）浓度与实测浓度对比如图8所示。区域排水中TP浓度波动较大，整体呈波动减小趋势，第1次峰值是由于基肥施用后首次由于降雨造成田块排水形成的径流携带了较多的磷，第2次峰值是6月下旬早施穗肥后降雨产生的排水携带了一定量的磷，经过几次排水的冲刷，磷流失较多，后期没有施用磷肥，因此排水中的磷含量平稳的减小。

注：SMS表示土壤微环境分层。

图7 典型支沟排水流量验证

图8 典型支沟排水总磷浓度验证

3.2.4 模型参数率定和模拟效果评价

经过参数敏感性分析，磷污染模型主要对表3中所列参数较为敏感。为取得最好的模型验证结果，率定的参数取值如表3所示。采用Nash-Suttcliffe效率系数（NSE）、RMSE和2作为模型模拟效果的评价指标，可得模型对灌区的水文水质模拟效果见表4。可见，水文模块的NSE、2均大于0.82，水质模块的NSE、2均大于0.81。模型性能总体上在预期之内，考虑到流域面积和土地使用的复杂性，模型效率大于0.8被认为是可以接受的。灌区磷污染物模型模拟结果可以反映水稻生育期间灌区水分运动和污染物迁移过程，模拟结果可靠。

表3 模型参数率定结果

表4 模型水文水质模拟效果评价

3.3 灌区磷污染负荷模拟统计

灌区2018年5—9月水稻田的磷产出负荷根据区域内水稻田的排水量和磷浓度过程计算，磷渗漏输出负荷根据各子流域渗漏量和渗漏磷浓度计算得到，并据此分析各个生育期的磷输入输出量，如表5所示。分蘖期和拔节孕穗期径流输出负荷分别为0.85 kg/hm2和0.60 kg/hm2。分蘖期由于降雨导致磷输出负荷较大，拔节抽穗期虽排水较多，但是经过分蘖期排水的冲刷，排水磷浓度较小，输出负荷较分蘖期小。其他生育期无排水。但泡田期和分蘖期渗漏输出负荷最大，分别为0.11 kg/hm2和0.16 kg/hm2。主要是由于基肥和早施穗肥导致土壤水浓度较高所致。水稻生育期内通过排水和渗漏流失的磷为1.88 kg/hm2，约占施肥和灌溉磷输入的5.7%。各主要排水支沟控制农田地表排水负荷、渗漏输出负荷及排沟末端输出负荷如图9所示。由图可以看出灌区一排和四排输出的磷负荷总量最大，分别为1.40 t和1.39 t。这是因为一排和四排控制面积较大，水田和居民点较多，污染物输出大于其余排水支沟。

表5 水稻生育期内磷的输入输出平衡

图9 和平灌区2018年5—9月各排水支沟及控制区域磷负荷

4 结 论

1）本文构建了适合平原灌区的非点源磷污染机理模型，该模型考虑水稻田实际水文循环过程和基于土壤微环境分层的磷运移过程，建立的水分运动方程和磷迁移转化方程比较符合平原灌区实际情况，能准确模拟平原水稻灌区的水量和水质变化。

2）利用黑龙江省和平灌区2018年试验田土壤水、积水层排水和典型支沟的排水流量、总磷浓度过程对模型进行率定与验证，试验田积水水深、排水、一排和七排径流量模拟效果指标Nash-Sutcliffe效率系数分别为0.967、0.980、0.820和0.861，决定系数分别为0.976、0.981、0.815和0.877；试验田、一排和七排排水磷浓度模拟效果指标Nash-Sutcliffe效率系数分别为0.811、0.890、0.826，决定系数分别为0.821、0.930、0.910，所有指标均大于0.81，模拟效果较好，表明该模型适用于以水稻为主要作物的平原灌区。

3）根据模拟的灌区水文水质过程计算了整个生育期稻田输入输出量及各排水支沟输出负荷，分蘖期和拔节孕穗期径流输出负荷分别为0.85、0.60 kg/hm2，泡田期和分蘖期渗漏输出磷负荷最大，分别为0.11、0.16 kg/hm2。水稻生育期内通过排水和渗漏流失的磷约占施肥和灌溉磷输入的5.7%。其中，灌区一排和四排输出的磷负荷总量最大。

[1] 李怀恩，沈晋，刘玉生. 流域非点源污染模型的建立与应用实例[J]. 环境科学学报，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] 陶园，王少丽，管孝艳，等. 青海省农业面源污染源特征分析[J]. 农业工程学报，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] 耿润哲，张鹏飞，庞树江，等. 不同气候模式对密云水库流域非点源污染负荷的影响[J]. 农业工程学报，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] 钟科元，陈兴伟，陈莹，等. 地形和土壤/土地利用空间聚合对径流和输沙量影响的AnnAGNPS模拟[J]. 农业工程学报，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模型在水稻灌区的改进及应用研究[J]. 水电能源科学，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] 郑捷，李光永，韩振中，等. 改进的SWAT模型在平原灌区的应用[J]. 水利学报，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] 胡文慧，李光永，郭亚洁，等. 汾河灌区农业面源污染经验统计模型的构建与验证[J]. 中国农业大学学报，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] 徐保利，代俊峰，俞陈文炅，等. 漓江流域氮磷排放对水肥管理和下垫面属性变化的响应[J]. 农业工程学报，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] 郝守宁，普布次仁，董飞. 林芝畜禽养殖粪便排放时空演变及耕地污染负荷分析[J]. 农业工程学报，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模型的水文水质过程模拟研究--以中田河流域为例[D]. 南京：南京农业大学，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] 杨霞，邵东国，徐保利. 东北寒区黑土稻田土壤水分剖面二维运动规律研究[J]. 水利学报，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] 汤鸿霄，钱易，文湘华，等. 水体颗粒物和难降解有机物的特性与控制技术原理：上卷水体颗粒物[M]. 北京：中国环境科学出版社，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] 刘锦涛，黄万勇，杨士红，等. 加气灌溉模式下稻田土壤水溶解氧的变化规律[J]. 江苏农业科学，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] 董志勇. 环境流体力学[M]. 北京：科学出版社，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模型的平原河网区水文过程分布式模拟[J]. 农业工程学报，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

陈明洪，白璐，刘宣冶，等. 基于土壤微环境分层的平原水稻灌区磷污染模型[J]. 农业工程学报，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

国家重点研发计划项目（2016YFC0400107）；国家自然科学基金面上项目（52079137）

陈明洪，副教授，主要从事流域水环境模拟与调度研究。Email：chenminghong@cau.edu.cn

10.11975/j.issn.1002-6819.2020.20.003

X52

A

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