基于DRAINMOD模型的不同灌排模式稻田水氮运移模拟

谢阳村，徐 敏*，高世凯，沈连起

（1.生态环境部 环境规划院，北京 100000；2.华北水利水电大学 水利学院， 郑州 450046；3.寿光市水利局，山东 寿光 262700）

0 引 言

【研究意义】氮素是农作物生长发育的关键因素，在传统灌排模式下, 由于灌溉和降雨造成的排水较多，大量氮素直接进入水体，导致农业面源污染[1-4]。已有研究表明，单一水稻控制灌溉与控制排水技术可高效利用养分和水分，减轻环境负荷，实现节水高产、减排、控污的目标[5-6]。但将水稻控制灌溉与控制排水技术进行有机结合，形成水稻全生育期的控制灌排技术[4,7]，进行节水减排控污效果的研究却较少。此外，农业面源污染区域发生具有随意性、污染排放不确定性以及空间异质性等特点，对其影响过程、影响因素以及产生的机理研究造成难度[8-9]。【研究进展】辅助于完善农田排水管理以及氮素流失的农田排水模型已经引起国内外的广泛关注[10-12]，但对于稻田的排水和氮素运移模型较为少见。DRAINMOD 模型是美国农业部推荐的农田排水模型，模型不断发展更新，并逐渐得到推广应用[13-15]。部分学者利用DRAINMOD模型尝试应用于稻田，认为该模型能够较好地模拟灌区稻田径流排水过程和地下水位动态变化[16-18]，也有学者利用DRAINMOD 模型对暴雨排水过程氮素变化进行探究[19-20]。高学睿等[21]将DRAINMOD 模型应用到气候相对湿润的湖北漳河灌区，对灌区稻田的排水和氮素流失规律进行模拟，发现模型可有效地模拟预测南方典型灌区稻田排水特性及氮素运移特性，相对误差范围为2.8%～16.0%。Hashemi 等[22]发现DRAINMOD 模型可以较好地模拟不同排水深度和间距组合下稻田硝态氮损失和排水流量，相对误差范围为1.0%～35.0%。【切入点】综上，对于DRAINMOD模型在稻田的应用研究多基于传统灌排模式，但对于不同灌排模式下稻田的应用研究较少。【拟解决的关键问题】采用Morris 方法对DRAINMOD 模型的水氮运移参数进行灵敏度分析，利用DRAINMOD 模拟不同灌排模式稻田水氮动态变化规律，探讨DRAINMOD模型在水稻不同灌排模式下的适用性，以期为我国南方稻作区灌排实践提供科学依据。

1 材料与方法

1.1 试验概况及农艺措施

试验研究于2017年和2018年6―10月在涟水县水利科学研究站的试验大田完成，试验区位于江苏省淮安市涟水县朱码镇王二庄村境内，地理坐标东经119°16′，北纬33°50′，属于亚热带湿润性气候区，多年平均气温14.4 ℃，年降水979.1 mm，降水量时程分布不均，汛期降水量占全年降水量比重较大，年均蒸发1 385.4 mm，年均日照时间2 280 h，年均无霜期240 d。试验区0～100 cm 土层土壤为壤土，土壤体积质量为1.42 g/cm3，pH 值为6.82，田间持水率为27.9%（质量），有机质量为2.19%，全氮量、全磷量分别为0.98、1.12 g/kg，速效氮量为21.5 mg/kg。采用格田小区试验，格田规格为90 m×27 m，格田田埂高30 cm以上，底部宽40 cm，上部宽30 cm，除农沟一侧的田埂外，格田周边嵌入50 cm 的薄膜并覆盖至田梗，避免了各小区间水位影响。试验区布置见图1。

图1 试验区布置 Fig.1 The layout of the test area

供试水稻品种为两优9918，该品种株型松紧适中，群体整齐，茎秆坚韧，叶色淡绿，叶片内卷挺直，熟期落色好，抗病虫害性能较好，耐高温能力强。该品种有效穗数241.5 万穗/hm2，穗总粒数175.7粒，结实率85.5%，千粒质量27.5 g，株高可达122.7 cm 左右，全生育期平均为130.8 d，多年平均产量为597.1 kg。水稻农艺措施见表1，复合肥N、P、K 质量比为15∶15∶15，尿素含氮量为46.4%。

表1 2017―2018年水稻农艺措施 Table 1 Agronomic measures of rice in 2017―2018

1.2 试验设计

考虑到南方地区水稻生长期与汛期同季，控制灌排模式可以在保持较低灌水下限的同时，提高雨后蓄水上限，扩大稻田的储水库容，达到充分利用天然降水，减少灌排水频率、灌排定额和氮磷排放量，并依据前人连续几年的试验结果[7,23,24]及江苏省水稻节水灌溉技术规范，各处理水位调控方案见表2。

表2 各处理农田水位调控方案 Table 2 Water control program of each treatment

选取水稻分蘖期、拔节孕穗期、抽穗开花期和乳熟期4 个生育阶段进行控水试验，共3 个处理，每个处理设3 个重复。传统灌排（CK）农田水位达到灌水下限值-20 cm 立即灌水至6 cm，雨后水层深度超过6 cm 即排水。控制灌排分为轻旱控排（LCID）和重旱控排（HCID），LCID 达到灌水下限值-20 cm 立即灌水至6 cm，分蘖期雨后水层深度超过10 cm 即排水，其他生育期雨后水层深度超过20 cm 即排水；HCID 达到灌水下限值-50 cm 立即灌水至6 cm，雨后排水方式与LCID 相同。

1.3 试验观测指标及方法

1）气象资料。由田间试验小区附近的自动气象站收集，主要包括降雨时间、降雨历时及降雨量和大气温度、大气压强、大气相对湿度、风速、水面蒸发量、有效辐射量等。

2）农田水位。每天09:00 对农田水位进行观测。当农田地面有水层时，通过竖尺在固定观测点测量农田地面水层深度。无水层时，通过在试验田块中间等距离安装的3 个地下水位观测井记录各小区农田水位的埋深。

3）灌水量。通过在农渠上安装的三角量水堰测定各格田灌溉水量。

4）排水量。采用水位差法，记录排水前后格田水层深度，根据二者差值计算排泄水量，其中当降雨引起的径流时应加上排水时段内的降雨量，日内发生多次排水按1 次计算。

5）水样提取及分析。在格田小区选择4 个取样点，地表水用塑料注射器不扰动土层分别在水面、1/2水深、底部取水，混合装入采样瓶，一般5～7 d 取1 次水样，但施肥后以及暴雨后1 周内隔日取样；渗漏水通过地下水观测井，深度为1.5 m，井管底部密封，管壁打孔并用无纺布过滤，井口加盖防雨。取样时先用真空泵清空井内积水，1～2 h 后吸取井管中渗漏水，取样频率同地表水。排水时在排水出口取水，每30 min 取1 次，当产流均匀时可适当延长步长。所有水样采样后贴好标签，并进行冷藏（3 ℃）处理，在24 h 内进行水质分析。水样中铵态氮、硝态氮分别采用絮凝沉淀纳氏试剂光度法、紫外分光光度法，所用仪器主要为UV2800 岛津紫外分光光度仪。

2 模型原理及主要参数

2.1 模型原理

DRAINMOD 模型以简单的土壤剖面范围内水平衡原理为基础来描述农田水文变化过程，可以准确地预测出农田水位、地表径流和地下排水量以及作物产量等。在时间增量Δt内，土壤水量平衡方程可表示为：

式中：ΔV为土体水量的变化量（cm）；D为侧向排水量（cm）；ET为蒸发蒸腾量（cm）；DS为深层渗漏量（cm）。

在单元Δt时段内，地表水量平衡方程为：

式中：ΔS为地表蓄水量的变化量（cm）；P为降雨量（含灌水量）（cm）；F为入渗量（cm）；RO为地表径流量（cm）。

DRAINMOD 模型可以用于模拟湿润地区浅水位条件下农田非饱和区（一维，垂向）及饱和区（二维，垂向和侧向）的土壤水和氮素转化运移，模型采用多相一维对流弥散方程表达式来模拟计算农田中的氮素运移，并用一阶有限差分法求解。多相一维对流弥散方程式为：

式中：θa为土壤液相体积分数；θg为土壤气相体积分数；ρb为土壤固相干体积质量（M/L3）；Ca为某种形式氮素液相浓度（M/L3）；Cg为某种氮素形式气相浓度（M/L3）；Cs为某种形式氮素固相浓度（M/M）；Da为水动力扩散系数（L2/T）；dg为分子扩散系数；Va液相体积通量；S为源汇项；t为时间（T）；z为空间坐标（L）。

表3 土壤参数和排水参数 Table 3 Main input parameters of soil and drainage

表4 氮素参数 Table 4 Nitrogen parameters

2.2 模型参数

土壤参数通过入渗试验、野外钻孔法、环刀法等方法测定。由于原系统针对的是暗管排水，在应用于明沟排水时需要进行相应的参数调整，排水深度和排水间距试验测定；排水管有效半径和地面到不透水层深度通过模型率定；由于土壤质地和试验条件相似，排水系数和Kirkham积水深度参考文献[25]确定。氮素相关参数通过模型率定和参考文献[26-27]确定。作物参数有效根系深度根据2017年和2018年实测数据以及相关经验进行设定。模型初始参数见表3—表5。采用2017年试验观测数据进行参数率定、2018年数据对模型进行验证。

表5 作物参数 Table 5 Growth parameters

3 模型参数敏感性分析

3.1 模型参数选择及分布

敏感性分析分为局部敏感性及全局敏感性，其中全局敏感性考虑了参数之间的相互影响，更加适用于DRAINMOD 模型。将参数设定为均匀分布，其参数分布范围参考文献[25]，具体见表6。

表6 模型参数服从的分布 Table 6 Obedient probability distribution of model parameters

水分运移模拟中水位变化选择通过衡量模拟值与观测值的拟合程度进而评估参数敏感性，即Nash Sutcliffe 效率（Nash）、相对误差总量RE以及相关系数CC，计算式为：

式中：Qi为农田水位实测值（cm）；Mi为农田水位模拟量（cm）；n为模拟次数；为农田水位实测平均值（cm）；为农田水位模拟平均值（cm）。

其他参数采用定性分析的Morris 法[28]进行敏感性分析，该方法通过参数变化效应的平均值μ和标准差σ来判断参数的敏感性，其中μ值反映参数敏感性大小，μ值越大，代表参数对模型输出结果的影响越大；σ值反映参数间交互作用的大小或者表示参数为非线性影响，σ值越大，代表参数交互作用越大。

3.2 水分运移模拟敏感性分析

对于选取的水分运移参数，依据Morris 方法共进行20 次抽样，得到100 组模型参数的抽样样本，运行模型得出各个参数的灵敏度，结果如图2 所示。

由目标函数均值μ可以得出，I.L对模型农田水位输出结果影响最小，De影响最大，三层侧向饱和导水率在不同目标函数上有所区别，但Ksat2影响最大，尤其在Nash和CC上响应较明显；对于地表排水量和地下渗漏量，Ksat2与De的影响最大，I.L影响最小。

图2 水文参数均值与标准差散点 Fig.2 Mean value and standard value of hydrological parameters

由目标函数标准差σ可以得出，整体呈现出“山谷”型，De受其他参数的影响（或非线性作用）最为明显，而Ksat2基本不受其他参数的影响；对于地表排水量和地下渗漏量，I.L处于较高水平，De和Ksat2受其他因素影响较小。

3.3 氮素运移模拟敏感性分析

对于选取的氮素运移参数，依据Morris方法共进行20次抽样，采样220组数据得到220组模型参数的抽样样本，运行模型得出各个参数的灵敏度，如图3所示。

由目标函数均值μ可以得出，对于地表铵态氮径流损失，λ对输出结果的影响最明显，因为扩散是铵态氮运移的主要机理，当铵态氮垂向运移距离较短时会导致较多的铵态氮通过径流方式排出，其次为Toptdec，因为当土壤温度与分解温度相近时，有利于有机质分解为铵态氮；对于地表硝态氮径流损失，反硝化反应起主导作用，反硝化参数Vmaxden、Kmden、Toptden对输出结果影响最大，有机质分解参数Toptdec、Kdeca、Kdecs影响最小；对于地下铵态氮渗漏损失，输出结果基本只受到有机质分解和硝化反应参数影响，依次为Kmnit、Toptdec、Vmaxnit、Kdeca、Topt-nit、Kdecs；对于地下硝态氮渗漏损失，输出结果基本只受到硝化反应及反硝化反应参数影响，依次为Vmaxden＞Kmnit＞Vmaxnit＞Kmden＞Toptden＞Toptnit。

由目标函数标准差σ可以得出，整体呈“山峰”型趋势，地表铵态氮径流损失参数Toptdec、Toptnit、Kdeca、λ相互影响较大；地表硝态氮径流损失参数中仅Kmden略受其他参数影响；地下铵态氮渗漏损失Toptnit和Kdeca相互影响较大；地下硝态氮渗漏损失Kmnit受其他参数影响大。

图3 氮素参数均值与标准差散点 Fig.3 Mean and standard values of nitrogen

图4 农田水位变化模拟值和实测值对比 Fig.4 Simulated and measured values of farmland water level

表7 农田水位拟合性能 Table 7 Fitting performance of farmland water level

4 模拟结果与分析

4.1 不同灌排模式下稻田水分动态

各处理农田水位变化模拟值和实测值对比如图4所示。从图4 可以看出，实测的农田水位与模型模拟结果之间吻合较好，模拟结果能较好地反映田间农田水位的变化情况。其中，正水位模拟程度较好，但负水位（即水位低于土壤表面）则模拟结果存在较大差别，但属于可接受范围。由表7 可以得出，农田水位拟合性能参数Nash在0.81～0.97 之间，均大于0.8，模拟效果较好；Re绝对值在0.01～0.53 之间，且主要为地下水位模拟误差；CC在0.86～0.98 之间，模拟值与实测值变化趋势较为一致。

不同灌排模式灌排水量模拟值和实测值见表7。各处理灌水量和排水量模拟值与实际值差别小于10%，模拟效果较好。与CK 相比，2017年LCID 和HCID 灌水量分别减少16.1%、37.1%，排水量分别减少68.6%、72.6%；2018年LCID 和HCID 灌水量分别减少9.7%、26.2%，排水量分别减少33.0%、41.9%。说明控制灌排模式通过合理调蓄雨水资源，减少稻田排水次数和排水量，对于节约灌溉用水有重要意义。而HCID 与LCID 相比，灌水量更少，说明节水省工效果更加显著。

4.2 不同灌排模式下稻田氮素动态

各处理氮素径流损失累计量如图5 所示。从整体趋势看，各处理模拟值与实测值在前期径流损失量较大，在移栽后几天发生较大变化是由于生育前期连续遇雨，而返青期和分蘖初期需要保持浅水层，因此地表排水量较大，且此时基肥的施入引起地表水中氮素浓度较高。其中，CK 径流损失量最大，主要是由于前期排水量较多，而控制灌排拦蓄雨水，减少了氮素的排放。与实测值相比，铵态氮径流损失模拟较好，虽然在2018年存在些许差异，主要是由于排水影响，而非氮素浓度导致，硝态氮径流损失模拟存在较大差别。

图5 氮素径流损失累计量 Fig.5 Cumulative nitrogen runoff loss

氮素径流损失累计总量模拟值与实测值对比见表8。各处理氮素径流损失累计总量模拟值与实际值差别小于20%，模拟效果较好。从铵态氮径流损失累计总量看，铵态氮模拟值与实际值相差不大，且与地表排水量模拟值变化情况保持一致；硝态氮径流损失累计总量在2017年各处理均有一定程度的升高，而2018年CK 模式下有所降低，这主要是模型对排水中硝态氮质量浓度模拟存在偏差所致。与CK 相比，2017年LCID 和HCID 铵态氮径流损失累计总量分别减少41.2%、45.0%，硝态氮径流损失累计总量分别减少49.6%、53.1%；2018年LCID 和HCID 铵态氮径流损失累计总量分别减少43.6%、47.4%，硝态氮径流损失累计总量分别减少29.8%、32.8%。说明控制灌排模式约束了地表水高质量浓度氮素时期的稻田排水，实现了减排控污的目标。

表8 不同处理氮素径流损失累计总量 Table 8 Fitting performance of nitrogen runoff loss

5 讨 论

模型参数的合理选取对于模拟结果的精确性起到关键作用，而大多数参数是依靠试验数据率定间接取得的，采用灵敏度分析来评价各个参数的不确定性对模拟结果的影响，确定出比较敏感的参数，有针对性地对这些参数进行测量和调节，从而选取合理的模型参数[29]。通过DRAINMOD 模型进行模拟时，20～40 cm 土层侧向饱和导水率Ksat2的变化对稻田水分运移结果影响较大，但与其他参数相互作用不明显，进行模拟输入时，可配合其他参数在适当范围内进行微调，来保证模拟精度。不透水层深度I.L和排水管有效管径De极易受其他参数影响，故其初始输入值的精度直接影响模型的模拟效果。铵态氮模拟的准确性取决于弥散系数、硝化反应参数、有机质适宜分解温度，硝态氮模拟的准确性取决于反硝化参数和硝化参数[19,21]，应该加强试验手段保证现场测量的精度。

DRAINMOD 模型进行模拟计算时，当农田水位为负值时，模拟值与实测值之间存在较大偏差，这可能是因为地下水位较浅时，土壤水力学特性对地下水位的预测比较敏感[14]，因此模拟效果存在偏差，但并未显著影响到灌水及灌水后地表水层深度。与采用 HYDRUS-1D 模型模拟稻田中水分运移的研究相比较[24]，DRAINMOD 模型表现良好。铵态氮模拟值与实际值相差不大，且与地表排水量模拟值变化情况保持一致，说明铵态氮径流负荷损失模拟偏差主要是由排水量所引起的。而DRAINMOD 模型低估了生育前期排水中硝态氮浓度，高估了中后期硝态氮浓度，使得曲线呈现差别，还需进一步研究。

本研究对不同灌排模式下灌排水量和氮素径流负荷进行模拟，结果发现控制灌排模式将水稻控制灌溉技术与控制排水相结合，较传统灌排模式，可高效利用养分和水分，充分发挥稻田的湿地效应，减少灌排水量和稻田氮磷污染物负荷，实现节水减排的目标，与前人研究结果[28]一致。同时，已有研究发现控制灌排通过调控稻田水分状况，可以实现水稻高产，与传统灌排相比，减产少于10%，且差异不显著[7,30-31]。但是，对水稻田来说，由于田间环境有很大差异，稻作区的沟、田水体具有相互影响、协同控制的特点，DRAINMOD 模型对于不同灌排模式的适用性需要进一步在不同田间尺度、不同水文年型和不同地理分区中进行检验与推广。

6 结 论

1）20～40 cm 土层侧向饱和导水率Ksat2对稻田水分运移模拟结果影响最为明显，同时各个参数间存在非线性作用，以不透水层深度I.L和排水管有效管径De受其他参数影响最为明显；弥散系数λ、硝化反应参数（Kmnit、Vmaxnit、Toptnit）、有机质适宜分解温度Toptdec对稻田铵态氮运移模拟结果影响较大；反硝化参数（Vmaxden、Kmden、Toptden）和硝化参数（Kmnit、Vmaxnit、Toptnit）对稻田硝态氮运移模拟结果影响较大。

2）DRAINMOD 模型能够较好模拟不同灌排模式的稻田水位，拟合度在正水位时较高，但在负水位时较低；DRAINMOD 模型能够较好地模拟不同灌排模式氮素运移，其中铵态氮的拟合度较好。控制灌排模式通过合理调蓄雨水资源，约束地表水高浓度氮素时期的稻田排水，可以实现节水减排的目的。