由于普遍增施化肥以及不合理的农业管理措施，造成氮肥以挥发和淋溶等形式损失增大，硝态氮的淋失是农田氮素损失的重要途径之一，由此导致的对地下水资源污染问题日益突出^{[1, 2]}。土壤剖面残留较高的氮素易形成垂直向下运移的趋势^[3]，并伴随灌溉与降雨下渗造成地下水污染^{[4, 5]}。如何在保护地下水水质和防控非点源污染的同时合理进行农田施肥已成为国内外研究热点^{[1, 6, 7, 14, 15]}。目前，计算和评估农田淋溶的方法包括直接测定法和其他模型法：前者包括淋溶盘或土水势法，由于其仅能测定点上数据，区域代表性差；后者在国内正在推广改进以及需要对国内区域环境进行耦合调整，实用性才能得到提高^{[8, 9, 12]}。模型模拟法因具有操作简单、模拟准确和代表性广等优点而成为评估和预测淋溶的核心方法，目前应用较多的模型有GLEAMS（Groundwater Loading Effects of Agricultural Management Systems）^[8]、SWAT（Soil Water Assessment Tool）^{[9, 10]}等，其在综合分析作物生产和土壤氮素运移方面还有待加强^{[11, 12]}。本文选取Root Zone Water Quality Model（RZWQM）模型偶合了农业生产管理以及环境影响的模块，成为预测和评估农田氮素淋溶的新工具^{[12, 13]}。

RZWQM模型是美国农业部大平原系统研究所（USDA-ARS，Great Plain System Research Unit）于1992年推出的农业系统作物和环境管理模型。该模型整合了作物根区所有影响作物生长的物理、生物和化学过程^{[15, 16]}，在我国华北地区已经进行了一些应用尝试，如Wang等^[17]用其优化污水灌溉条件下冬小麦-夏玉米施肥措施，Ma和Fang等^{[18, 19, 20]}模拟了不同灌溉制度下作物产量及水分利用效率。已有的研究多侧重于灌溉和施肥量的利用效率研究方面，应用模型对农田施肥过程硝态氮淋溶的评估虽然已有研究，但还处于初级阶段^{[21, 22]}。该模型对田间实测数据指标要求高，根据实测数据对模型参数进行率定和验证后，可以在华北平原其他数据缺乏的田间进行更广泛的验证或推广^[23]，为华北地区农田硝态氮淋溶预测开辟新的思路和估算方法。

冬小麦-夏玉米轮作试验在河北省农林科学院大河试验站进行。该地区属于中温带大陆性季风气候，年均温度为12.9 ℃。年均降水量为496.6 mm左右，其中春季降雨量较少，7—9月降水量占全年总降水量的70%。

田间施肥试验始于2003年，但从2009年开始气象数据比较齐全，故本文以2009年田间实测数据进行模型参数率定，以2010和2011年试验数据进行模型参数验证。土壤类型为粘壤质洪冲积石灰性褐土，土壤剖面基本理化性状见表 1。其中初始硝态氮和初始铵态氮数据在小麦收获后，玉米播种前测量，作为模型的本底值。设置4个冬小麦-夏玉米轮作周期施氮量，分别为0（N0）、215（N1）、400（N2）、575 kg·hm^-2（N3）。每个处理设4次重复，随机区组排列。小区面积为44 m2。氮、磷、钾肥分别为尿素、过磷酸钙和氯化钾。小麦季氮肥施用量分别为0、150、250、350 kg·hm^-2，磷肥（P₂O₅）225 kg·hm^-2，钾肥（K₂O）105 kg·hm^-2，氮肥于小麦播前和拔节期分两次施用（底追各半），磷钾肥做底肥一次施用。玉米季对应施氮量分别为0、65、150、225 kg·hm^-2，钾肥（K2O）施用量52.5 kg·hm^-2，氮肥在玉米大喇叭口期和抽雄期分两次施用，钾肥在大喇叭口期一次施用。玉米收获时以小区计算产量，且采集1 m深土壤剖面样品，土样采集时间分别为2009年10月12日、2010年10月9日和2011年10月12日。挖取小区1 m深土层剖面，每20 cm为一层，分别用环刀取单位体积土样进行数据测定，分别测定土壤含水率以及土壤硝态氮含量，采用 1 mol·L^-1KCl浸提比色法测定硝态氮。

表 1 土壤基本参数数据 Table 1 Soil basic parameters used

表选项

RZWQM模型是美国农业部大平原系统研究所研发的能够对土壤根区水质、土壤养分运移以及作物生长进行综合模拟的模型，由6个模块组成：物理模块、化学模块、养分模块、作物生长模块、杀虫剂模块和管理模块。RZWQM模型包括“天”和“时”两个时间尺度：以“天”为尺度计算离子、肥料、灌溉水、耕地措施，以及潜在蒸散发和蒸腾量；水分迁移及营养化学过程以“时”尺度计算，包括土壤水分的再分配、营养盐的迁移转化、渗滤、径流、杀虫剂淋洗、热量损耗、实际蒸散发、植物养分吸收过程等，然后计算杀虫剂迁移转化过程（如果涉及的话）、碳氮迁移转化过程和土壤物质平衡过程，最后运行植物生长模块。RZWQM模型能模拟土壤氮转化的主要过程，包括作物氮残留、有机态的矿化、氮素固定、氮素淋溶、氨挥发、硝化和反硝化等。应用该模型可以较好地模拟氮素的迁移转化，获得氮素淋溶数值和评估淋溶情况更加方便快捷，从而有针对性地实行优化施肥措施以削减损失。

模型初始输入数据包括气象数据和土壤基本数据。其中气象数据包括模拟期间日降水、日最高温、日最低温、风速、相对空气湿度等；土壤基本数据包括土层容重、田间持水量、土壤pH值、土壤含水率、剖面硝态氮和铵态氮初始含量（表 1）。

模型参数的率定遵循先水分模块，然后养分模块，最后作物模块的顺序^[24]。模型率定效果的评价是判定参数优化的关键，不同统计指标各有优缺点^[25]，本文在模型参数率定过程中采用归一化误差（均方根差/平均值）进行参数优化^[11]，即不断改变参数值以减小模拟结果和实测结果的差异，以归一化误差数值最小时模型选定的参数作为最终率定结果。选用两个指标评价模型运行的结果：均方误差（RMSE），属于绝对误差指标，反映模拟效果的绝对无偏性和极值效应；平均相对误差（MRE），属于相对误差指标，反映模拟效果的相对无偏性。模型参数率定和验证效果通过RMSE和MRE体现，当RMSE达到最小值而MRE趋近于0，表明模拟效果较优。计算公式如下：

应用2009—2011年N3处理剖面实测水分数据率定模型参数，以N2、N1和N0处理实测数据验证模型参数（图 1）。2009—2011年玉米收获时土壤剖面含水率模拟值与实测值变化趋势相同。2009年在0～20、20～40、40～60、60～80、80～100 cm土层的率定结果RMSE分别为0.001 9、0.006 0、0.002 5、0.007 6、0.000 2 cm3·cm-3，MRE分别为-1.13%、-3.43%、-1.35%、4.22%、-0.11%。在2010—2011年N2、N1和N0处理剖面含水率验证结果表明RMSE分别为0.012 2～0.012 4、0.011 5～0.011 7、0.011 0～0.011 3、0.004 6～0.006 1、0.010 2～0.014 8 cm3·cm-3，MRE分别为-7.81%～ -7.53%、-8.54%～8.23%、-8.80%～-7.42%、-4.00%～3.08%和7.07%～10.11%。

图 1 N3处理土壤剖面含水率模拟值与实测值对比 Figure 1 Comparison of simulated and measured soil water contents in N3 treatment

图选项

方差分析结果表明，上部土层含水率实测值和模拟值能够较好地吻合，效果优于下部土层，可能原因为模型没有良好地模拟地表水分蒸散量与水分下渗作用造成的。但实测值与模拟值之间差异不显著，模型水分模块可以较好地反映土壤实际含水率。

2009—2011年N3处理土壤剖面硝态氮模拟值与实测值变化趋势相同（图 2），0～20、20～40、40～60、60～80、80～100 cm土层RMSE变化范围分别为4.32～7.30、1.45～7.47、2.46～6.16、2.91～5.84、3.94～6.01 mg·kg-1，MRE变化范围分别为 -11.7%～27.7%、-4.9%～44.7%、-10.5%～33.8%、-16.2%～0.1%、-17.9%～26.7%。N3处理土壤剖面硝态氮模拟值可以较好地反映土壤剖面硝态氮实测值。

图 2 N3处理土壤剖面硝态氮模拟值和实测值对比 Figure 2 Comparison of simulated and measured nitrate concentrations in N3 treatment

图选项

应用2009—2011年N2、N1和N0处理剖面土壤硝态氮实测值对模型参数进行了验证。N2处理土壤硝态氮模拟值与实测值变化趋势与N3处理相似（图 3），且方差分析结果表明实测值与模拟值之间差异不显著，对土壤剖面硝态氮含量模拟效果可以接受。此外，N2、N1和N0处理模型剖面硝态氮含量模拟效果见表 2。统计分析结果表明模型可以模拟不同氮肥施用量条件下土壤剖面硝态氮含量。

图 3 N2处理土壤剖面硝态氮模拟值和实测值对比 Figure 3 Comparison of simulated and measured nitrate concentrations in N2 treatment

图选项

表 2 2009—2011年不同施氮量条件下硝态氮验证参数 Table 2 Nitrate concentrations in different nitrogen applications from 2009 to 2011

表选项

对土壤剖面含水量和硝态氮含量进行参数率定和验证后，就可以进行作物产量率定和验证。2009—2011年作物产量模型模拟值与实测值比较见图 4。各处理的产量模拟值均小于实测值，不同处理间产量的模拟值之间差异不明显，而且不同处理间实测值之间差异也不明显。这主要是因为试验开始前施氮量很大，导致土壤残留和积累的氮素已经能够满足作物生长需要。

图 4 2009—2011年不同施氮量条件下产量关系比较 Figure 4 Comparison of simulated and measured yields in different nitrogen applications from 2009 to 2011

图选项

对土壤剖面水分和氮素，以及作物产量等参数率定和验证后，模型主要参数就被确定下来，可以开展不同施氮量条件下土壤氨挥发和硝态氮淋溶的预测研究。通过模型对土壤氨挥发进行的预测表明，土壤氨挥发与施氮量呈正相关关系。轮作周期内施氮量分别为575、400、215、0 kg·hm^-2时，玉米季土壤氨挥发分别为109.21、70.68、32.04、0.64 kg·hm^-2。氨累计挥发损失约占全年施氮量的17.2%。预测结果稍高于其他文献资料，具体模拟过程还有待细化^{[26, 27]}。

土壤剖面1 m以下硝态氮淋溶预测结果见图 5。可以看出，2009—2011年土壤剖面硝态氮淋溶量随施氮量的减少而降低，在2010年不同施氮量处理硝态氮淋溶量对应比其他年份有所降低。这是由于2010年全年降雨量为432 mm，均小于其前后两年200 mm左右，同时也说明土壤硝态氮淋溶与降水量密切相关。轮作周期内施氮量分别为575、400、215、0 kg·hm^-2时，玉米季硝态氮淋溶量平均为79.71、51.04、26.33、3.45 kg·hm^-2。

图 5 2009—2011年不同施氮量条件下预测的1 m以下土层硝态氮淋溶情况 Figure 5 Predicted nitrate leaching below 1 meter in different nitrogen applications from 2009 to 2011

图选项

本文运用RZWQM模型对不同施氮水平下华北农田土壤水分氮素迁移转化以及氨挥发和淋溶进行了模拟和预测。通过依次对水分模块、养分模块和作物生长模块参数的率定和验证表明模拟结果在可接受范围内，再次印证RZWQM模型已经可以应用于平原地形氮素迁移转化规律的研究^[28]，也进一步证实RZWQM模型可以预测和评估平原地区施肥条件下土壤淋溶和氨挥发。当土壤类型为粘壤土、饱和导水率较低、pH值偏碱性时，模型参数率定过程中需要注意初期土壤含水量与实测值进行耦合调节，这是后期准确模拟剖面硝态氮含量的基础^{[29, 30]}。此外，模型对不同施氮量水平下土壤剖面硝态氮模拟效果有一定的偏差，RMSE值最大达到12.04 mg·kg-1，可能因为前期施肥量过高导致土壤剖面氮素含量累积，对此问题还需要更进一步提高模型模拟效果。总体来看，模型参数率定和验证的效果可以接受。施氮量显著影响土层中硝态氮的累积，在土壤结构和降水的作用下会明显发生淋溶^[31]，对环境及地下水造成污染。在封丘国家试验站有学者得出，在冬小麦-夏玉米轮作体系下，全年氮肥施用量最大为400 kg·hm^-2，本研究3年施肥情况显示，2 m内不同土层硝态氮可累计达到260~418 kg·hm^-2，如继续增加施氮量，将会造成严重的地下水源污染^[32]。在本文不同施氮量处理下，3年土壤平均残留硝氮量分别为397.54（N3）、225.76（N2）、109.03（N1）、6.88 kg·hm^-2（N0）。基本符合华北地区情况。

对硝态氮淋溶和氨挥发过程的模拟表明，玉米季土壤硝态氮淋溶量随施氮量下降而呈显著降低。轮作周期施氮量分别为575、400、215 kg·hm^-2时，2009—2011年硝态氮淋溶量均值分别占轮作周期总施氮量的13.9%、12.8%和12.2%，分别占玉米季施氮量的35.4%、34.0%和40.5%。表现出玉米季硝态氮淋溶较大，而相对前茬小麦作物当季玉米施氮量却较低。主要因为以往的研究多以轮作周期为对象，分析冬小麦-夏玉米轮作周期内氮素迁移转化，而忽略了不同作物间的硝态氮淋溶相对重要程度。本研究发现虽然前茬小麦施氮量大于玉米季，反而玉米季淋溶量较大。这是因为华北地区降水主要集中在7—9月，麦季土壤存留的大量硝态氮在玉米季随降水下渗和淋失。由此推测，虽然麦季施氮量过高但淋溶可能相对玉米季较小，而在玉米季发生了硝态氮的损失和污染。此结论还有待结合整个轮作周期的模型模拟进行深入分析。氨挥发量和施氮量也具有正相关关系，3年氨挥发均值分别占轮作周期总施氮量的18.9%、17.7%和15.0%。模型预测结果基本符合当地情况。

模型不仅可对不同施氮量土壤剖面水分和氮素迁移情况进行模拟，同时也可以评估硝态氮淋溶损失以及作物产量。有研究表明，施肥过量会增加土壤硝态氮含量，增加氮素淋溶的风险和通过挥发途径的氮素损失^[33]。模型模拟结果表明，2009年玉米季施氮225 kg·hm^-2作物产量达到最大，但是与其他处理间差异没有达到显著水平，而在模拟其他处理时存在一定偏差。这可能是因为农田养分投入和作物养分吸收间存在不平衡^[34]，并且土壤本底值偏高，在不施肥的情况下作物利用土壤中原有氮肥而表现出对施氮量反映不明显。2010年玉米季施氮量65 kg·hm^-2下，作物产量同样能保持在较高水平，全年总施氮量为215 kg·hm^-2左右就可以获得较高的玉米产量。从玉米产量、氮素的利用以及土壤中硝态氮积累综合考虑，夏玉米合理施氮量应该控制在195～260 kg·hm^-2范围内^{[35, 36]}。本试验结果进一步证明了该结论，在华北地区应进一步降低施氮量，同时也证实RZWQM模型在华北平原地区对农田氮素淋溶评估的可操作性。而RZWQM模型可以有效模拟作物氮素迁移、淋溶和氨挥发，在实测数据参数验证基础上，为提供合理有效控制氮素淋溶的施肥量提供基础。

通过华北平原的玉米季土壤剖面含水量和硝态氮含量以及作物产量等实测数据，对RZWQM模型参数进行了率定和验证，并分析了该模型用于评估硝态氮淋溶和氨挥发的可行性。结果表明，冬小麦-夏玉米轮作周期施氮量分别为575、400、215和0 kg·hm^-2时，2009—2011年土壤剖面硝态氮平均均方误差分别为6.56、7.30、3.64、1.53 mg·kg-1，平均相对误差分别为26.7%、17.7%、32.7%和-37.2%。硝态氮淋溶和氨挥发平均损失量占轮作周期施氮量的13.0%和17.2%，且随施氮量的减少而降低。本研究对大田条件下土壤氮素运移的模型模拟取得了较为理想的结果。在此结果上，推算轮作周期施氮量215 kg·hm^-2时可以有效控制硝态氮的淋溶损失。此外，进一步完善RZWQM模型小麦季的参数率定对于拓展模型的预测功能极为重要。