硝酸盐是无机氮肥在土壤中的主要形式之一, 对作物的生长发育起着重要作用^[1], 硝酸盐极易溶于水, 当施肥量超过植物需要及土壤反硝化能力时, 氮素就会在降雨或灌溉作用下随土壤溶液运移, 通过淋溶进入含水层污染地下水^{[2, 3]}.氮素的淋失不仅造成农业水资源严重浪费, 化肥和农药等的大量流失, 也导致了地下水水质的退化, 成为引起地下水污染最常见、分布最广泛的农业非点源污染物之一^{[4, 5]}.如何在保护地下水水质和防控非点源污染的同时合理进行农业施肥成为国内外学者研究的热点问题^[6~8].

太行山低山丘陵区是华北平原的生态屏障也是华北平原地下水补给区, 属于典型的农林交错地带, 人类活动频繁、土地利用类型多样.该地区土壤结构不良, 石砾含量大, 且多属粗骨土, 易侵蚀.土壤风化层薄, 土层厚度在20~50 cm之间, 土壤保肥保水能力差, 具有较高的渗透性能^[9].近年来低山丘陵区农业活动如经济林和农田面积大幅度增加, 氮肥过度施用以及人类污染物的排放造成该区域地下水硝酸盐浓度逐年升高^[10].因此研究太行山低山丘陵区典型农田土壤氮淋溶过程对保护补给区地下水具有重要意义^[11].

目前, 计算和评估农田氮淋溶的方法包括直接测定法和其它模型法：前者包括淋溶盘或土水势法, 由于其仅能测定点上数据, 区域代表性差; 后者在国内正在推广和改进, 同时在应用时需要针对国内某些特殊的区域环境进行耦合调整来提高模拟的精度^{[12, 13]}.模型模拟法因具有操作简单、模拟准确和代表性广等优点而成为评估和预测淋溶的核心方法, 目前应用较多的模型有GLEAMS(groundwater loading effects of agricultural management systems)^{[14, 15]}、SWAT(soil water assessment tool)^{[16, 17]}等, 其在综合分析作物生产和土壤氮素运移方面还有待加强^{[18, 19]}.而根区水质模型(root zone water quality model, RZWQM)耦合了农业生产管理以及环境影响的模块, 成为预测和评估农田氮素淋溶较为普遍的工具^{[20, 21]}.

因此, 本文利用RZWQM模型模拟太行山低山丘陵区冬小麦-夏玉米农田轮作系统土壤硝态氮和水分运移状况, 并且通过模拟3种不同施氮量以及不同降水年型的情景下, 分析土壤硝态氮淋溶状况.本研究对于预测山区农田氮淋溶量对地下水的影响以及合理控制山区农田施肥保护地下水环境等方面具有一定意义.

冬小麦-夏玉米试验在中国科学院太行山山地生态试验站(114°15′50″E, 37°52′44″N)进行.该地区属于暖温带半干旱大陆性气候, 年平均气温13℃(-1.6℃~26.3℃), 年平均降水量513.2 mm, 其中春季降水量较少, 6~9月降水量占全年降水量的75%, 年平均蒸发量为1 598 mm, 海拔为49~812 m^[22]. 图 1为2015~2016太行山降水量变化特征, 其中2015年降水量为469.4 mm, 2016年降水量为927.1 mm, 2016年7月19~20日降水量为504.2 mm.利用太行山站过去30年的降水数据进行统计分析, 表明2015年为平水年, 2016年为丰水年, 同时根据降水划分标准判定2016年7月19~20日降水属于大暴雨即极端降水^{[23, 24]}.

图 1

Fig. 1

图 1 2015~2016年太行山站降水量变化特征 Fig. 1 Precipitation at Taihang Station during 2015 and 2016

该地区上游为石质山区地质结构为前寒武系滹沱群变质岩组, 出露的岩石主要为花岗片麻岩, 在其表层覆盖有0.5~10 m厚的风化层.土壤主要为褐土, 土层薄, 角砾石含量大, 易侵蚀, 具有较高的渗透性能.该试验小区自1986年建站以来施肥灌溉措施与山区农民施肥灌溉相一致, 小区土壤理化性质基本稳定.土壤剖面基本理化性质和土壤机械组成见表 1.

表 1 (Table 1)

表 1 研究点土壤基本理化性状和机械组成特征1) Table 1 Basic physical and chemical properties and mechanical characteristics of the soil in the study 土壤深度/cm 土壤理化性质 机械组成/% 土壤容重 /g·cm-3 饱和导水率 /cm·h-1 田间持水率 /cm3·cm-3 土壤温度 /℃ 初始硝态氮 /μg·g-1 黏粒 粉粒 砂粒 0~10 1.25±0.02 4.80±0.2 25.72±1.2 20.5±0.5 9.16±1.0 2.3±0.1 37.1±2.1 60.6±3.8 10~20 1.27±0.03 4.52±0.1 23.39±1.5 19.1±0.3 11.49±1.3 2.0±0.2 36.8±1.8 61.2±3.2 20~30 1.30±0.05 4.25±0.4 21.52±1.7 18.0±0.2 4.69±0.9 2.3±0.1 34.7±2.2 63.0±2.8 30~50 1.32±0.02 4.13±0.3 17.62±1.1 17.4±0.4 8.32±0.7 2.1±0.3 36.9±2.3 61.0±3.1 50~70 1.50±0.04 3.88±0.2 18.11±1.9 17.0±0.5 16.22±1.8 2.4±0.2 37.4±1.9 60.2±3.5 70~100 1.56±0.03 2.18±0.2 14.75±1.6 15.1±0.3 20.01±2.1 2.4±0.1 39.5±2.4 58.1±2.9 1)表内数值表示平均值±标准差(n=3)

表 1 研究点土壤基本理化性状和机械组成特征1) Table 1 Basic physical and chemical properties and mechanical characteristics of the soil in the study

冬小麦-夏玉米耕作方式为常规翻耕.冬小麦-夏玉米轮作的施氮量如下：冬小麦施氮量(以N计)为300 kg·hm^-2, 小区面积为310 m²(长22.6 m, 宽13.7 m), 同时在小区四周设置宽度为1m的保护行, 氮肥于小麦播种前和返青期分2次施用, 底追各半, 玉米季不施肥.冬小麦灌溉2次分别在种植前和返青期, 灌溉方式采用当地农民大水漫灌的灌溉方式, 玉米季无灌溉.

土壤水动态监测利用含水量反射传感器(Water Content Reflectometer, Model CS615, Campbell Scientific, Inc., Logan, UT, USA)监测, 每个传感器由两根20 cm的波导管(waveguides)组成.埋设深度分别为10、20、30、50、70和100 cm, 传感器探针水平放置, 每30 min监测一次, 所得土壤含水量为体积含水量.

土壤硝态氮含量分别在作物生长期(冬小麦生长初期、返青期和收获期; 夏玉米生长初期和收获期)利用剖面采样和负压法抽提土壤水测定.剖面采样测土壤硝态氮含量利用离心萃取法, 即称取风干土样品10 g, 用超纯水按1:5浸提, 测定上清液的EC, 用0.25 μm的醋酸纤维素滤膜过滤, 采用离子色谱仪(ICS-2100, Dionex, 美国)测定NO₃^-含量.

RZWQM模型是美国农业部大平原系统研究所(USDA-ARS, Great Plain System Research Unit)于1992年推出的农业系统作物和环境管理模型.该模型整合了作物根区所有影响作物生长的物理、生物和化学过程, 主要应用到我国平原地区, 能够准确模拟土壤剖面、土壤表层, 浅层或者深层渗漏过程中的硝酸盐含量^{[25, 26]}. RZWQM主要由物理、化学、养分、杀虫、作物和管理这6个模块构成, 其中包括“日”和“时”两个时间尺度, 以“日”为尺度计算离子、肥料、灌溉水、耕地措施, 以及潜在蒸散发和蒸腾量^[27].水分迁移及营养化学过程以“时”尺度计算, 包括土壤水分的再分配、营养盐的迁移转化、渗滤等.

模型初始输入数据包括当地全年的气象数据和土壤基本数据(如表 1).其中气象数据包括模拟期间日降水、日最高温、日最低温、风速、相对空气湿度等; 土壤基本数据包括土层容重、田间持水量、土壤pH、土壤含水率、剖面硝态氮和铵态氮初始含量.

本文以2015年田间种植的冬小麦实测数据输入RZWQM模型进行模型参数率定, 以2016年夏玉米模型输出数据与实测值进行比较, 然后进行RZWQM模型参数验证.模型参数的率定包括水分、养分和作物模块.

模型率定效果的评价是判定参数优化的关键, 不同统计指标各有优缺点^[28], 本文在模型参数率定过程中采用归一化误差(均方根差/平均值)进行参数优化.公式如下：

(1)

(2)

(3)

式中, RMSE表示均方误差; RE_i表示第i个实测值的相对误差; MRE表示平均相对误差; N表示物理量实测总个数; P_i表示模拟值; Q_i表示实测值.率定结果标准主要通过RMSE和MRE表示, 表现为RMSE达到最小值, 而MRE趋近于0.

利用2015年冬小麦整个生育期土壤剖面实测水分数据进行模型参数的率定, 以夏玉米剖面实测数据进行模型的参数验证.将实测的土壤饱和含水率、残留含水率、田间持水率、饱和导水率、水分特征曲线等土壤水力参数作为初始输入条件, 然后根据实测值与模拟值的比较分析, 相应地调整各土壤剖面水力参数的取值, 主要是饱和含水率、田间持水率、饱和导水率, 从而保证模型结果符合率定准则.小区各层土壤水力参数实测值与调整值见表 2.

表 2 (Table 2)

表 2 土壤水力参数实测值与率定后调整值 Table 2 Observed and simulated values of the soil hydraulic parameters 土壤深度/cm 饱和含水率/cm3·cm-3 田间持水率/cm3·cm-3 饱和导水率/cm·h-1 实测值 模拟值 实测值 模拟值 实测值 模拟值 0~10 0.36 0.357 24.59 25.723 4.85 4.796 10~20 0.39 0.382 21.92 23.395 4.59 4.516 20~30 0.43 0.419 20.81 21.521 4.20 4.253 30~50 0.44 0.432 17.01 17.626 4.05 4.125 50~70 0.49 0.514 17.93 18.114 3.58 3.881 70~100 0.46 0.475 13.82 14.755 2.02 2.175

表 2 土壤水力参数实测值与率定后调整值 Table 2 Observed and simulated values of the soil hydraulic parameters

从图 2冬小麦土壤剖面含水量模型率定实测值与模拟值的比较可以发现, 实测值与模拟值变化趋势一致.其中, 土壤剖面表层土壤含水率受降水或灌溉影响明显, 而模拟结果都大致呈现出表层含水率较大的波动过程.下层70 cm土壤含水率波动过程明显不如表层强烈, 底层100 cm土壤含水率较为稳定, 模拟结果体现了不同土层含水率的变化特征. 2015年冬小麦在0~10、10~20、20~30、30~50、50~70和70~100 cm土层的含水量率定结果RMSE分别为0.017、0.014、0.015、0.010、0.011和0.004 cm³·cm^-3, MRE分别为8.81%、8.50%、6.43%、8.26%、6.42%和2.64%.

图 2

Fig. 2

图 2 冬小麦土壤剖面含水率模拟值与实测值对比 Fig. 2 Comparison of the simulated and measured soil water contents in the winter-wheat season

利用率定调整后的参数, 验证2016年夏玉米土壤含水量及硝酸盐分布.夏玉米土壤含水量验证结果见图 3.从中可以看出, 模拟值反映了实测值的变化过程且方差分析结果表明, 实测值与模拟值之间差异不显著, 对土壤剖面土壤含水量模拟效果可以接受. 2016年夏玉米在0~10、10~20、20~30、30~50、50~70和70~100 cm土层的含水量验证结果RMSE分别为0.031、0.022、0.020、0.025、0.019和0.013 cm³·cm^-3, MRE分别为7.3%、6.1%、5.6%、7.3%、4.9%和5.36%.方差分析结果表明, 下部土层含水率实测值和模拟值能够较好地拟合, 效果优于上部土层, 这主要是由于上层水分受到地表水分蒸散量与水分下渗相互作用的影响变化较大.但整体上实测值与模拟值之间差异不显著, 模型水分模块可以较好地反映土壤实际含水率.

图 3

Fig. 3

图 3 夏玉米土壤剖面含水率模拟值与实测值对比(2016年) Fig. 3 Comparison of the simulated and measured soil water contents in the summer-maize season (2016)

图 4为2015年冬小麦不同生长期土壤硝态氮模拟值与实测值拟合结果.从图 2中可以看出, 模拟值均反映了硝态氮随小麦生长期的变化趋势, 生长前期上层硝态氮含量低, 下层高, 收获前期体现出表层含量较高, 下层土壤硝态氮含量逐步降低的特点, 从时间上也大致体现出硝态氮先升高后降低的趋势.另外, 总体上, 下层土壤硝态氮模拟结果优于土壤表层, 这也是由于上层土壤受施肥、灌水、降雨等因素影响较大的原因. 2015年冬小麦在0~10、10~20、20~30、30~50、50~70和70~100 cm土层不同时期的土壤硝态氮率定结果平均RMSE值分别为3.60、4.48、5.75、3.14、3.39和2.68 mg·kg^-1, MRE平均值分别为8.9%、4.5%、2.2%、1.9%、3.5%和0.9%.

图 4

Fig. 4

图 4 2015年冬小麦土壤硝态氮率定 Fig. 4 Calibration of the soil nitrate nitrogen content in the winter-wheat season of 2015

利用2016年夏玉米土壤硝态氮实测值对模型参数进行了验证.夏玉米土壤硝态氮模拟值与实测值变化趋势相似(图 5), 2016年夏玉米在0~10、10~20、20~30、30~50、50~70和70~100 cm土层不同时期的土壤硝态氮验证结果平均RMSE值分别为5.64、5.75、6.54、6.59、1.77和4.47 mg·kg^-1, MRE平均值分别为6.3%、2.6%、3.6%、3.3%、2.3%和0.8%.且方差分析结果表明, 实测值与模拟值之间差异不显著, 对土壤剖面硝态氮含量模拟效果可以接受.统计分析结果表明, 模型可以模拟不同氮肥施用量条件下土壤剖面硝态氮含量以及硝态氮淋滤量.

图 5

Fig. 5

图 5 夏玉米土壤剖面硝态氮模拟值与实测值对比 Fig. 5 Comparison of the simulated and measured soil nitrate contents in the summer-maize season

图 6显示了2015年冬小麦生育期内土壤剖面硝态氮含量变化.小麦季0~30 cm土层内土壤硝态氮在生育期内存在显著变化, 总体上0~30 cm土层硝态氮含量明显高于当季冬小麦初始硝态氮含量且该土层内出现硝态氮浓度峰值, 其峰值大小随冬小麦生育期的推进而增加(即峰值随着氮肥的施入次数增加而增加).冬小麦返青灌溉前(2016年4月16日)30 cm以下土层土壤硝态氮浓度与当季冬小麦初始值相近, 说明小麦季返青前30 cm土层以下土壤硝态氮含量基本维持稳定状态.这表明冬小麦返青期之前降水、施氮肥和作物吸氮等主要影响0~30 cm土层内土壤硝态氮浓度, 未影响到30 cm以下土层土壤硝态氮含量.冬小麦返青期进行灌溉施肥措施后, 30 cm以上土层土壤硝态氮含量会继续升高, 而30 cm以下土壤硝态氮含量会急剧降低.研究发现土壤硝态氮的淋溶是随着土壤水分的排出而出现^[29], 本研究中由于返青期灌溉水以及小麦后期(5~6月)的几次降水入渗补给到30 cm以下土层(图 2), 因此30 cm以下土层部分土壤硝酸盐会随着降水的入渗而淋失进入地下水.

图 6

Fig. 6

图 6 冬小麦生育期内土壤剖面硝态氮分布 Fig. 6 Concentrations of nitrate nitrogen in the soil profile in the winter wheat season

图 7显示了2016年夏玉米生育期内土壤硝态氮剖面分布.从中可以看出夏玉米生育期内土壤硝态氮分布变化规律与冬小麦不同(图 6), 冬小麦生育期内0~30 cm土层土壤硝态氮发生变化, 30 cm以下在小麦返青期之前没有变化, 小麦返青期后随着降雨量增加, 30 cm以下有部分硝态氮淋滤出100 cm土层; 而夏玉米季30 cm以上土层和30 cm以下土层土壤硝态氮均显著变化, 30 cm以上土层土壤硝态氮先随着施氮肥增加而增加, 而后由于大暴雨引起的淋滤而降低, 30 cm以下土层土壤硝酸盐总体上在玉米季初期保持不变, 而后受上层硝态氮随降水淋滤作用而导致下层增加, 最后随着降水淋滤到地下水而导致土壤硝态氮含量降低.

图 7

Fig. 7

图 7 夏玉米生育期内土壤剖面硝态氮分布(2016年) Fig. 7 Concentrations of nitrate nitrogen in the soil profile in the summer maize season(2016)

分别通过模型计算2015~2016年冬小麦和夏玉米季土壤硝态氮淋滤量, 本研究发现小麦季施氮量为300 kg·hm^-2时, 小麦季土壤硝态氮淋滤量仅为2.12 kg·hm^-2, 而玉米季土壤硝态氮淋滤量却达到了59.9 kg·hm^-2.因此在太行山低山丘陵区农田土壤硝态氮淋滤主要发生在玉米季(雨季), 小麦季硝态氮淋滤量较少.

对土壤剖面水分和硝态氮等参数率定和验证后, 利用RZWQM模型开展不同灌溉量或降水量、不同施氮量水平下土壤硝态氮淋滤的预测研究.通过以上模拟发现, 2015~2016冬小麦-夏玉米季土壤硝态氮淋滤主要发生在夏玉米季(雨季), 在冬小麦季淋滤量仅为2.12 kg·hm^-2, 所以本研究只模拟预测玉米季土壤硝态氮淋溶.通过太行山低山丘陵区农田施氮量调查可知, 当地农田施氮量主要有3种情况：全年不施肥(0 kg·hm^-2), 当地冬小麦平均施氮量(300 kg·hm^-2), 当地冬小麦较高施氮量(450 kg·hm^-2).另外由于该区域农田灌溉较少, 仅在小麦种植前和返青期灌溉, 而玉米季无灌溉.因此本研究只模拟不同施氮量和不同降水年型对土壤硝态氮淋溶量的影响.通过模拟太行山区农田3种不同施氮和不同降水年型情景下土壤硝态氮淋溶量, 以期选择更加合理的山区农田施氮量.轮作周期内全年施氮量分别为0、300及450 kg·hm^-2和不同降水年情景时, 玉米季土壤剖面1 m以下硝态氮淋溶预测结果见图 8.

图 8

Fig. 8

图 8 RZWQM模拟不同降水年和不同施氮量条件下1 m以下土层硝态氮淋溶情况 Fig. 8 Predicted nitrate nitrogen leaching below 1 m in different nitrogen applications in 2013, 2015, and 2016 using RZWQM

可以看出, 2013(枯水年)、2015(平水年)和2016 (丰水年)年土壤剖面硝态氮淋滤量随施氮量的增加而增加. 2013年不同施氮量处理土壤硝态氮淋滤量最低, 这是由于2013年全年降水量为382mm, 玉米生长期降水量为336mm, 均小于2015和2016年降水量, 同时也说明土壤硝态氮淋滤与降水量密切相关.

通过对土壤硝态氮淋滤的模拟表明, 玉米季土壤硝态氮淋滤量随施氮量下降而呈显著降低.通过对比太行山低山丘陵区3种施氮水平, 当施氮量低于300 kg·hm^-2时可显著减少土壤硝态氮的淋溶损失, 而施氮量较高时其硝态氮淋滤量也显著增加, 尤其是在丰水年特别是受强降水的影响淋滤量增加更显著.如模拟结果显示施氮量为450 kg·hm^-2时, 在2013和2015年玉米季平均硝态氮淋滤量为63.35 kg·hm^-2, 但是在2016年夏玉米季发生强降水(7月19日247 mm, 7月20日256 mm)导致其硝态氮淋滤量增加到了136.5 kg·hm^-2, Asadi等^[30]和夏梦洁等^[31]同样发现较强的降水导致较多的硝态氮淋滤量.当轮作周期施氮量分别为300 kg·hm^-2和450 kg·hm^-2时, 2013(枯水年)、2015(平水年)和2016(丰水年)年玉米季硝态氮淋溶量均值分别占轮作周期总施氮量的13.20%~19.63%和14.58%~21.96%.表现出玉米季硝态氮淋滤量较大, 而相对前茬小麦作物当季玉米施氮量却较低或者不施氮.本研究发现虽然前茬小麦施氮量大于玉米季, 反而玉米季淋溶量较大.这是因为华北太行山低山丘陵区降水主要集中在6~9月, 麦季土壤存留的大量硝态氮在玉米季随降水下渗和淋失.

由以上结果可见, 在太行山低山丘陵区应进一步降低施氮量和灌溉量, 可以减少土壤硝态氮的淋溶损失, 而RZWQM模型可以有效模拟作物氮素迁移和淋溶, 在实测数据参数验证基础上, 为太行山低山丘陵区合理有效控制硝态氮淋滤提供依据.

针对太行山低山丘陵区农田氮淋失引起的地下水硝酸盐浓度上升问题, 本研究通过典型小麦-玉米轮作农田土壤水分、土壤硝态氮观测, 利用RZWQM模型模拟预测了不同施氮水平在不同降水年型下1m土壤剖面硝态氮淋失量.取得的主要结论有：土壤硝态氮淋溶主要发生在夏玉米季(雨季), 土壤硝态氮淋失量与施氮量和降水量有显著关系, 平水年和枯水年在当地中等施氮水平下(300 kg·hm^-2), 氮的淋失量无显著变化, 高施氮水平(450 kg·hm^-2)受降水的影响氮的淋失量显著升高; 受特大暴雨的影响, 使得土壤氮的淋失量呈非线性增长的趋势, 表明大降水事件不仅对地下水形成大量补给, 很大程度上也增加了累积在农田土壤中的硝态氮的淋溶损失, 增加了对区域地下水硝态氮潜在污染威胁.