不同氮肥液浓度下微孔陶瓷根灌土壤水氮运移与流失特性的数值模拟研究

杨 雪，刘旭飞，韩梦雪，张 林

（1.中国科学院教育部水土保持与生态环境研究中心，陕西 杨凌 712100；2.中国科学院水利部水土保持研究所，陕西 杨凌 712100；3.中国科学院大学，北京 100049；4.西北农林科技大学水利与建筑工程学院，陕西 杨凌 712100）

0 引 言

施用氮肥是提高农作物产量最简单、有效的方式之一，已被作为关键影响因子用来调控作物生长发育和最终产量[1]。但是，在实际的农业生产中，农民经常为提高作物产量而盲目施肥。过高的氮肥投入不仅降低其利用效率，还会造成环境问题。我国作物的氮肥利用率普遍小于50%，远低于发达国家[2]。而未被利用的氮肥随着土壤水分的深层渗漏而淋失，加重地下水体的非点源污染，Niu 等人的研究表明地下肥水可能主要是由历史上有机肥等不断下渗的结果，导致深层包气带NO3-大量累积，加剧了地下水NO3-污染风险[3]。为解决氮肥利用率低的问题，水肥一体化技术应运而生。这一技术将灌水和施肥融为一体，借助压力灌溉系统，将水肥液均匀、准确地输送到作物根部，是氮肥增效的关键措施之一[4]。苏明辰等[4]调研发现，与传统施肥方法相比，水肥一体化技术将水分利用率提高了40%～60%，肥料利用率提高了30%～50%。袁念念等[5]研究表示水肥一体化技术将肥料集中输送到作物根部，增加了根系密集区域的氮素含量，减弱了硝态氮深层淋溶风险。然而，目前水肥一体化技术是周期性地将氮肥液施到作物根部，而氮素一旦进入土壤会被脲酶迅速催化水解，这种短时间内大量水肥供给的方式会增加氮肥无效损失风险，造成作物后期氮肥供应不足，导致氮肥的过量施用[6-8]。因此，推进高效水肥供施技术发展，对提高水肥利用效率具有重要意义。

近年来，西北农林科技大学自主研发了一种微孔陶瓷根灌技术（SICE）。相对于地下滴灌，SICE在大幅度降低工作压力（＜1.0 m）的同时，通过微孔陶瓷灌水器内部微米级的孔隙实现了小流量出流。在灌水量相同的前提下，大大延长了灌水时间，减缓了氮素进入土壤的速率，使氮素的肥效可以缓慢释放，较传统灌溉方式有效养分供应时间长，整个生育期内土壤水氮环境稳定，氮肥利用效率高。目前，已有研究论证了SICE运行可靠性和田间实际应用效果[9-11]。然而，以往研究多聚焦在土壤水分，目前利用SICE 进行氮肥液持续微量供给时的土壤水氮运移及其流失特性尚不明晰，限制了该技术在实际中的推广应用。

近年来，众多学者建立了大量土壤水氮联合运移模型。其中，由美国盐土实验室（US Salinity laboratory）开发的HYDRUS 模型最为广泛应用。Doltra 等对滴灌施肥条件下种植了不同作物的砂质壤土进行建模，结果表明该模型对0～90 cm 土层内铵态氮（NH4+-N）和硝态氮（NO3--N）含量的预测十分精确，平均绝对误差值低于观测值的平均标准差[12]。聂坤堃等利用HYDRUS-2D 反演所得参数对沟灌肥液入渗过程进行了模拟，结果显示土壤水分运移模拟值与实测值相对误差均小于7.0%，溶质运移相对误差均小于11.5%[13]。上述研究均获得了较高的拟合效果，由此可见，利用HYDRUS 模型对SICE 土壤水氮运移规律及流失特性进行预测是可行的。

为此，本研究采用HYDRUS-2D 建立了SICE 持续微量供给下的土壤水氮联合运移模型，并于中国青海开展了微孔陶瓷根灌枸杞田间试验，用以模型校准和准确性验证，最后以地下滴灌（CK）为对照，利用模型研究了3 种氮肥液浓度（211.54 mg/L、169.23 mg/L、126.92 mg/L）下SICE 土壤水氮运移与流失特性。研究结果可为阐明SICE 土壤水氮运移特性和水肥高效利用内在影响机制提供理论依据。

1 材料与方法

1.1 试验区概况

研究地点位于青海省海西州德令哈市怀头他拉镇（东经96°44'，北纬37°21'）（见图1），地处柴达木盆地东北边缘，平均海拔2 874 m。属大陆性高原气候。年平均气温2.4 ℃，昼夜温差大，最高平均温差可达15.8 ℃。年平均降水量110.3 mm，集中在5-8 月，无霜期97 d 左右，年平均蒸发量2 000～2 500 mm。试验开始前，取试验地0～100 cm 土壤进行检测，按照USDA 分类标准为壤质砂土。该地土壤保水能力差、养分含量低，土壤有机质含量为5.0 g/kg，速效氮含量108.6 mg/kg、速效磷含量为3.6 mg/kg，速效钾含量为2.9 mg/kg。其基本物理参数见表1。

图1 研究区位置及模拟的概念几何（单位：cm）Fig.1 Study area location and conceptual geometry of simulation

表1 土壤物理特性Tab.1 Soil physical characteristics

1.2 试验设计与方法

选取试验地的枸杞园区进行试验，品种为“宁枸七号”，其生育期共140 d，试验期为2019 年。种植时株行距为0.8 m×1.2 m。以SICE 实现水分和氮素的持续供给，陶瓷灌水器在西北农林科技大学旱区节水农业研究院中制备完成，尺寸为6.0 cm×4.0 cm×2.0 cm(长×内径×外径)，设计流量为0.25 L/h(工作压力为0.2 m时空气中的出流量)。按照一管一行的布置方式将陶瓷灌水器埋设在每棵枸杞根层土壤，埋深为30 cm（枸杞根系密度最大位置），距离树干20 cm，布设长度40 m。供水装置为体积1.1 m3的水桶，此外，该系统还由供水水箱、浮球阀、过滤器、施肥器、毛管、阀门等组成，其详细结构见文献[14]中的描述。以地下滴灌（CK）为对照，选用耐特菲姆超级台风贴片式滴灌带，设计流量为3.75 L/h。灌溉定额1 820 m3/hm2，全生育期共灌水7 次，每次间隔15～20 d，灌水定额260 m3/hm2，每次灌水6.67 h，每棵树灌水24.96 L。各处理行间设置0.2 m 缓冲区，并深埋塑料薄膜予以隔离，防止处理间土壤水分的横向交换。施用氮肥选用尿素，含氮量为46.7%。从高到低设置3 种施肥浓度，采用不完全组合试验，共包括4个处理（表2）。特别地，由于CK 属于间歇灌溉，土壤水氮情况呈现周期性的变化规律，再次灌溉会重复上一过程。此外，模拟时间尺度过长会忽略单次灌水后枸杞根系层土壤水氮变化细节，故本研究选择一个灌水周期（20 d）进行模拟。

表2 试验处理Tab.2 Experimental treatment

2 模型构建与验证

2.1 模型基本方程

2.1.1 土壤水分运动基本方程

微孔陶瓷根灌与地下滴灌的土壤水分运动均为三维流动问题，但该系统下水分入渗具有轴对称特点，即水分以陶瓷灌水器中心轴为对称轴均匀的向四周入渗。进一步地，假设模拟区域内的土壤性质各向同性且均匀；入渗水流为连续介质且不可压缩，颗粒骨架保持不变；忽略温度、水分滞后作用等因素以及溶质浓度变化产生的势梯度对水分运动及分布的影响[15]。故此时水分入渗可简化为二维平面问题。采用Richards方程描述土壤剖面水流运动：

式中：θ为土壤体积含水率，cm3/cm3；h为土壤负压力水头，cm；K(h)为非饱和导水率，cm/d；t为模拟时间，d；x为径向坐标，cm；z为垂直坐标，设向上为正，cm。

2.1.2 土壤氮素运移基本方程

土壤中氮素的迁移转化受多因子的综合影响。本研究选用尿素作为肥料，尿素分子进入土壤后在适宜的水分和酶的影响下，发生一系列生物化学反应。本研究忽略铵态氮挥发、以及土壤中气相和温度的影响，优先考虑水解、吸附、硝化、反硝化和矿化作用。进一步假设溶质在土壤溶液中是通过对流—弥散作用运移的。其中，矿化作用用零级动力学反应方程描述，水解、硝化与反硝化用一级动力学反应方程描述[15]。氮素运移基本方程详见公式（2）～（4）[16]。

尿素态氮运移方程：

铵态氮（NH4+-N）运移方程：

硝态氮（NO3--N）运移方程：

式中：C1、C2、C3分别为尿素态氮、铵态氮、硝态氮质量浓度，mg/cm3；k0，k1，k2，k3分别为有机质矿化、尿素水解、铵态氮硝化和硝态氮反硝化速率常数，d-1；kd为铵态氮吸附速率，cm3/g；Dxx和Dzz为水动力弥散系数张量的分量，cm2/d；包含纵向弥散系数DL、横向弥散系数DT，cm；和自由水中的分子扩散系数DW，cm2/d[17]；qx和qz分别为x、z方向的土壤水通量，cm3/d；s为固相铵态氮吸附量，s=kdC2；ρ为土壤干容重，g/cm3。

2.2 初始条件和边界条件

本研究模拟采用的概念几何和边界条件如图1所示，其为宽度120 cm，厚度100 cm 的土层。初始土壤溶质含量设为零，初始土壤含水量设定与试验地实际条件相同。边界条件分为水流边界和溶质运移边界。水流上边界设置为大气边界，输入通量包括降水、作物潜在蒸腾和棵间蒸发量，地下滴灌带和陶瓷灌水器均设置为变流量边界，输入通量包括灌溉和溶质通量。左右两侧设置为无通量边界。下边界假定为自由排水边界。对于溶质运移边界，由于滴灌带和陶瓷灌水器周围未形成积水，故此时边界条件采用第三类边界条件。同理，溶质运移的上边界和下边界同样选择第三类边界条件。

2.3 模型参数及验证

2.3.1 土壤水力参数

模拟过程中土壤含水率θ以及非饱和导水率K(h)等水力特征参数均采用Van Genuchten模型进行拟合[18]：

式中：θr和θS分别为土壤残余含水率和土壤饱和含水率，cm3/cm3；KS为土壤饱和导水率，cm/d；Se为相对饱和度，cm3/cm3；α为土壤进气的倒数；n为孔隙分布指数；l为空气联通系数，作为许多土壤的平均值，被估计为0.5[19]。

拟合率定后得到土壤水分特征参数如表3所示。

表3 土壤水分特征参数Tab.3 Soil hydraulic properties

2.3.2 土壤溶质参数

本氮素运移模型中忽略土壤中气相和温度的影响，氮素的运移参数包括纵向弥散系数DL、横向弥散系数DT和自由水中的分子扩散系数Dw；氮素转化包括矿化、水解、吸附、硝化、反硝化5 个参数[15,16]，并根据实测数据进行率定，见表4和表5。

表5 氮素转化参数Tab.5 Nitrogen transformation parameters

2.3.3 模型验证

本研究使用决定系数（R2）和相对均方根误差（NRMSE）评估土壤水氮运移模型计算准确度。R2越趋向于1表示模拟值与实测值越接近[20]，采用公式（9）进行计算；NRMSE的取值小于0.2 则表示模拟值与实测值具有高度匹配性，其计算方法如公式（10）所示：

式中：n为对比次数；Pi为模型值；Oi为实测值；Oˉ为平均实测值。

图2为土壤含水率和土壤NO3--N含量计算值与实测值对比结果。由图2可知，R2的拟合结果分别为0.891 1和0.854 9，表示该水氮联合运移模型分别解释了89.11%和85.49%的不确定性，说明实测土壤含水率和硝态氮含量中约有89.11%和85.49%可由对应的模拟值来说明或决定。此外，NRMSE的值均小于0.2。说明构建的模型能较准确的反映土壤水氮时空动态变化规律。在此基础上，可进一步研究不同处理对土壤NH4+-N 时空分布的影响，并对根层土壤水氮流失进行预测。

图2 土壤含水率和土壤NO3--N含量的模拟值与实测值的相关性Fig.2 The correlation between the simulated and measured values of soil water content and soil NO3--N content

3 模型应用

3.1 土壤水氮时间变化特征

图3为一个灌水周期结束时根系层土壤平均含水率、土壤平均NH4+-N和NO3--N含量在不同处理下的时间变化特征。

图3 不同处理土壤水氮随时间变化图Fig.3 The variation of average soil water and nitrogen content with time in different treatments

由图3 可知，对于土壤含水率的变化，CK 处理前期快速增加至0.1 cm3/cm3以上，紧接着快速降低至0.07 cm3/cm3左右，而T1、T2、T3 处理的土壤含水率始终维持在0.09 cm3/cm3左右，相比与CK 形成的波动的土壤水分环境，SICE 更能保证土壤水分环境的稳定。

对于土壤平均NH4+-N含量的变化，CK处理先快速增加到16.70 mg/kg，随后逐渐降低，在第13 天左右接近完全消散。而T1 处理土壤平均NH+4-N 含量增加到13.33 mg/kg 后下降，后续稳定在3.11 mg/kg。相比于CK，T1 处理土壤NH4+-N 含量下降得更慢，在根系层的留存时间更长。此外，T1 处理中土壤NH4+-N 含量峰值较CK 降低，这可能是由于前期CK 灌水量大，充足的水分有利用于尿素氮向铵态氮的转化，而T1 的出流量小，前期土壤水分不足，抑制了NH4+-N的产生。

对于土壤平均NO3--N含量的变化，呈现和土壤NH4+-N含量类似的变化规律，4 种处理土壤平均NO3--N 含量皆在经历一个峰值后逐渐降低。对比CK 和T1 发现，随着时间的增加，T1 较CK 处理土壤NO3--N 含量下降幅度小、在土壤中的浓度高，说明相较于CK，SICE 在整个周期内NO3--N 供应更加充足、稳定。此外，通过T1、T2、T3 对比发现，随着氮肥液浓度的降低，土壤NO3--N 含量的峰值降低，曲线斜率变缓，波动范围减小。而且，相较于土壤NH4+-N含量，土壤NO3--N含量达到峰值的时间相对滞后（从3 d 增加到7 d），这是因为施用的肥料是尿素的关系，会有一个尿素氮向铵态氮，再向硝态氮的转化过程。

综上所述，相较于CK，SICE 土壤平均水氮含量的变化速率明显变小，有效水分和养分供应时间长，保证了整个灌水周期水氮的充足，使得土壤含水量和含氮量更加稳定。

3.2 土壤水氮空间变化特征

图4 为一个灌水周期结束时（第20 d）根系层土壤水分、土壤NH4+-N和NO3--N在不同处理下的空间分布特征。

图4 第20天土壤含水率，土壤NH4+-N含量和土壤NO3--N含量的空间分布Fig.4 The spatial distribution of soil water content,soil NH4+-N content and soil NO3--N content on the 20th day

由图4 可以看出，一个灌水周期过后，CK 处理的根系层水分含量明显小于T1，T1 处理的湿润体基本保持不变，最大含水率的土层位于枸杞根系密度最大处（-30 cm），这有利于将氮肥液保存在枸杞根系层中被作物吸收。此外，通过对比T1、T2、T3 处理还发现氮肥液浓度的增加对土壤水分的分布无显著影响。

对于根系层NH4+-N的分布，CK处理的土层在一个灌水周期后NH4+-N 的含量接近于0，而T1 处理还在持续出流，使得枸杞根系密度最大处的NH4+-N含量稳定在3 mg/kg以上。对比T1、T2、T3，随着氮肥液浓度的增加，土壤各剖面NH4+-N 的含量增大，并对其在土层中的分布范围有轻微促进作用。并且，NH4+-N呈现聚集在灌水器周围的现象。

对于根系层土壤NO3--N 的分布，一个灌水周期结束后，CK 处理的NO3--N 含量降低至8.11 mg/kg，T1 处理在14.87 mg/kg 以上。对比T1、T2、T3 发现，氮肥液浓度的增加对土壤中NO3--N含量影响较大，而对其在剖面分布无明显影响。此外，与NH4+-N集中吸附在灌水器周围相比，土层NO3--N分布范围更广，并趋于向下层土壤运移，但始终集中在枸杞根系土层（0～-60 cm）。

上述研究成果均表明，和CK 相比，SICE 减小了根系层土壤水肥环境的波动，有效保证了土层根、水、肥的高度匹配。

3.3 土壤水氮流失分析

选择一个灌水周期进行模拟，计算获得的根系层土壤水氮流失结果如图5所示。

图5 根层土壤水氮流失特性Fig.5 Characteristics of water and nitrogen loss in root layer soil

对于根系层土壤水分流失。CK 和SICE 的土壤水分渗漏量分别为14.94、14.16 mm。一个灌水周期内两者相差较小，前期土壤水分渗漏曲线高度重合，直至第10 d 开始出现差异，10 d 后CK 处理的渗漏量持续保持最大。这可能是因为CK 短时间内大量供水，而T1、T2、T3 持续微量出流，导致CK 土壤水分损失曲线的斜率大于T1、T2、T3 处理，随着时间增加，差异累积，直至在第10 d 后显现出来。这里T1、T2、T3处理的根层水分损失高度重合。造成此现象的原因是模型中计算水分损失取决于输入的土壤特征参数，此设定相同，本模拟只改变了肥液浓度，忽略了实际土壤中各种因素的相互作用。

对于根系层土壤NH4+-N的流失。土壤NH4+-N淋失量在第5 d 开始出现差异，一个灌水周期内土壤NH4+-N 淋失量为CK（0.094 kg/hm2）＞T1（0.081 kg/hm2）＞T2（0.080 kg/hm2）＞T3（0.079 kg/hm2），随氮肥液浓度的增加淋失量增大，但各处理间相差较小，最大差值为0.015 kg/hm2。这是因为NH4+-N本身易与土壤胶粒结合，容易被大量吸附停留在灌水器附近的土层，向土壤下层运移的NH4+-N 含量很少，而且土层中的NH4+-N 会继续向NO3--N 转化，导致NH4+-N 淋失量差异不明显。

同理，根系层土壤NO3--N 变化趋势类似，只是时间相对滞后。不同处理土壤NO3--N 的淋失量结果相差较大，分别为CK（1.49 kg/hm2）＞T1（0.73 kg/hm2）＞T2（0.63 kg/hm2）＞T3（0.52 kg/hm2），与CK 相比，一个灌水周期内T1 处理每公顷减小了51.01%的NO3--N 淋失量。其中最大差值为0.97 kg/hm2。这是因为NO3--N 与土壤胶粒性质相似，不易被土壤吸附，灌水后会随着水分在土壤剖面内不断迁移，导致NO3--N 更容易淋失。

总体来看，SICE 处理土层水分渗漏量和氮素淋失量均小于CK 处理。其中，根系层土壤水分和NH4+-N 的淋失量差值较小，而NO3--N 淋失量则差异较大。对于根层水氮渗漏，SICE 表现明显优于CK 处理，在确保根系层土壤水氮环境稳定的同时渗漏量减小，降低了硝酸盐污染物向地下水迁移的风险。

4 讨 论

本研究采用HYDRUS-2D 建立了SICE 持续微量供给下的土壤水氮联合运移模型，用以研究SICE 土壤水氮运移与流失特性。SICE 相较于传统灌溉方式最大的不同在于灌溉持续时间和出流特性。当前以间歇灌溉为主的水氮供施方式，存在的最大问题即土壤水分和养分的剧烈波动。与之相反，SICE属于连续灌溉。在灌水量相同的前提下，可利用陶瓷灌水器内部微米级孔隙实现水肥液的持续微量供给[21]，减缓水分和氮肥进入土壤的速率，延长灌溉时间，使得肥效可以缓慢释放，保证了整个生育期内土壤水氮环境的稳定。

地下滴灌一次的灌水时间仅为数小时，其初期出流量大，由于土壤的各向同性此时重力势占主导地位，短时间内大量水肥供给的方式可能会导致优先流的产生，造成水氮的深层渗漏，致使根区氮素不足，无法满足作物生长后期的需求。相反地，SICE 持续微量的出流方式，更倾向于形成基质流，水分和养分均匀下渗，有效的将土壤NH4+-N和NO3--N保留在作物根系层。例如本文3.2 的研究中，相较于CK，SICE 的0～60 cm根系土层中明显有更多的氮素留存（图4）。

此外，根层氮素含量的增加，也可能是因为水的可用性增加，理想的土壤水分含量更利于尿素的转化，使作物更好地吸收养分。对比CK 和T1，在氮素运移和反应参数设定相同，且肥液浓度一致的情况下，T1 处理土壤NO3--N 含量峰值大于CK 处理（图3）。这是因为不同形态的氮在土壤中会相互转化，当尿素随水施入土壤后，先在脲酶的作用下生成碳酸铵[22]。碳酸铵进一步与水发生反应转化为NH4+-N，NH4+-N 再氧化为NO3--N。由于前期CK灌水量大，形成了饱和土壤水环境，加快了生成NO3--N 这一反应进程。而硝酸根带负电荷，不易被带负电荷为主的土壤胶体吸附，其主要存在于土壤溶液中，容易随水移动而引起硝酸根的淋失[23]。因此，CK 将较于SICE 处理的NO3--N 淋失量更大（图5）。而且，CK 处理前期土壤湿度过大，在通气不良情况下，根层中的NO3--N 可经反硝化作用还原成不能被作物利用的氧化亚氮，氧化氮和氮气，造成脱氮损失[24,25]。

研究还发现，土壤NH4+-N和NO3--N空间分布存在较大差异，土壤NH4+-N 的分布半径明显小于NO3--N。一方面是因为NO3--N 和生物作用使得NH4+-N 在土层中短暂存在后就转化为NO3--N。另一方面是因为区别于NO3--N 的化学特性，铵离子带正电荷，容易与土壤胶粒上已有的阳离子进行交换，并吸附于土壤胶粒，不仅吸附在土壤表面，还可进入粘土矿物的晶体中，成为固定态铵离子[26]，移动性较小，比较容易被土壤“包存”[27]。也容易成为交换态养分，逐步供给作物吸收利用[28,29]。因此相较于NO3--N 较广的分布范围，NH4+-N 集中吸附在灌水器周围。

此外，本研究并未发现氮肥浓度对土壤水分的运移有明显影响，与此相反，有研究表示土壤NO3--N 浓度会影响土壤水分分布特性，从而导致NO3--N 浓度高的土层吸水量增加[30]。因此，未来有必要进一步考虑灌溉水和氮肥之间的相互影响，以更好地优化SICE的灌水施肥策略，推进其工程化应用。

5 结 论

综上，本文利用HYDRUS-2D 构建了SICE 持续微量供给下的土壤水氮联合运移模型，并于中国青海开展了微孔陶瓷根灌枸杞田间试验，用以模型校准和准确性验证，最后以地下滴灌（CK）为对照，利用模型研究了3 种氮肥液浓度（211.54 mg/L、169.23 mg/L、126.92 mg/L）下SICE 土壤水氮时空分布特征，并对土壤水分渗漏及氮素淋失动态进行了模拟分析，得出以下结论：

（1）校验后的HYDRUS-2D 模型在预测土壤含水量、NO3--N 运移过程中表现良好，进一步地，可以对SICE 土壤NH4+-N 的时空分布规律，以及土壤水分深层渗漏和氮素淋失进行评估。

（2）与CK 相比，SICE 减缓了氮肥液进入土壤的速率，延长了有效水分和养分的供应时间（从6.67 h 延长至480.00 h），降低了根系层土壤水分、NH4+-N 和NO3--N 的时间波动性，保证土壤水氮相对充足的同时还维持了土壤水氮环境的相对稳定；其次，氮肥液持续微量供给期内，SICE 将土壤水分、NH4+-N 和NO3--N 维持在枸杞根系密度最大处的土层深度，充分保证了土层根、水、肥的高度匹配，为枸杞水氮高效利用提供了有效保障。

（3）对于根层土壤水分渗漏量、土壤NH4+-N和NO3--N淋失量，三者的处理结果均为CK＞T1＞T2＞T3，与CK 相比，T1处理每公顷减小了51.01%的NO3--N淋失量。SICE在降低水分深层渗漏和氮素淋失方面明显优于CK，降低了硝酸盐污染物向地下水迁移的风险。