滴灌施肥条件下土壤水氮运移数值模拟

黎会仙，王文娥，胡笑涛

(西北农林科技大学旱区农业水土工程教育部重点实验室,陕西 杨凌 712100)

水、肥是影响作物生长的2个重要因子，在实际生产中如果灌溉和施肥方式不合理，将造成水、肥浪费及地下水污染。滴灌施肥是一种节水节肥的先进技术[1]，目前滴灌随水施肥过程中土壤水分、养分的运移分布规律大部分是基于试验获取，而试验方法费时费力。计算机数值模拟具有简便、灵活等特点，能实现对各种灌溉情况的预测[2]，建立适当的数学模型可以更高效准确地研究滴灌施肥条件下水分、养分的运移转化过程。近年来，相关学者建立了不同滴灌条件下土壤水分的运移模型[3-4]，并且滴灌条件下土壤水和溶质运移模型都已取得发展[5-6]。脱云飞等[7]建立了膜孔肥液自由入渗的土壤水氮运移模型，较好地反映了土壤水氮运移转化规律。有学者利用 HYDRUS-1D 软件对不同情景下土壤水氮的迁移变化进行动态模拟，均获得了拟合度较高的水氮模拟模型[8-9]。EU-Rotate_N和HYDRUS-2D 模型均能模拟氮的淋溶损失，较精确地预测氮含量，在模拟氮的吸收及运移方面各自具有优势[10]。Tournebize等[11]应用SWMS模型模拟了植被覆盖条件下氮的运移。Gärdenäs等[12]模拟了不同微灌系统、灌溉施肥措施和土壤类型对硝酸盐淋失的影响。热合木等[13]模拟了膜下滴灌土壤水、盐分运移，获得较好的效果。

以上研究者应用不同的模型对土壤水、溶质在不同灌溉施肥条件下的运移进行研究，均获得较好的效果，但现有的土壤水、溶质运移研究大多仅考虑单一溶质的运移。滴灌施肥条件下土壤氮的动力学过程非常复杂，与土壤温度、pH、含水率、微生物、土壤类型等多种因素有关。

尿素是广泛使用的肥料，且在土壤中易分解转化，由于尿素态氮在土壤中运移转化的复杂性，目前对滴灌施肥条件下土壤尿素态氮的运移转化模拟较少，因此，有必要对尿素态氮在土壤中的运移变化进行模拟计算，以便对滴灌系统水肥的运移进行预测。本研究基于试验基础，利用HYDRUS软件对滴灌施肥条件下土壤水氮(尿素)运移及再分布过程进行模拟分析，并利用滴灌水肥入渗试验对模型进行验证，以确定模型对滴灌施肥条件下土壤水氮运移分布模拟的精度，为滴灌施肥系统的优化设计提供参考依据，提高水肥利用效率，促进增产增收。

1 试验设计与方法

1.1 土壤物理性质

供试土壤取自西北农林科技大学试验田0～40 cm土层，土壤基本物理参数见表1。

1.2 试验装置与方法

试验装置包括马氏瓶、圆土柱(高50 cm，半径35 cm)、滴头、EM50传感器及数据采集器，如图1所示。试验所用肥料为普通尿素(CON2H4)。

试验设置滴头流量分别为 2、3、4 L·h-1，灌水量为21 L，施肥浓度为300 mg·L-1。试验前将风干土过2 mm筛，土壤初始体积含水率为0.057 cm3·cm-3，初始土壤铵态氮、硝态氮含量分别为6.5、25.05 mg·kg-1。试验土柱分层装土，装土时进行搅拌，各层之间进行打毛、平整，以保证土柱内土壤分布均匀，各向同性均质。装土过程中埋设传感器，每个土柱埋设12个EM50传感器，分层布置在土体中，根据滴灌土壤水分入渗特点[14-15]及传感器测试精度和范围，传感器的布置方式见图1。装土完毕后静置24 h，使土体初始含水率均匀稳定。灌水前在土柱表面覆盖一层薄膜。试验前对所有仪器进行标定。

表1 土壤基本物理特性

注:图中数字所示为传感器分布位置。Note: The numbers in the figure show the distributed location of sensor.图1 试验装置及传感器布置Fig.1 Experimental equipment and sensor placement

试验开始后通过数据采集器监测土壤含水率随时间的变化情况，滴灌过程中每5 min监测一次，滴灌结束后每10 min监测一次。在滴灌结束24、48、72、96、120 h分别在5个纵向半剖面上用土钻按传感器布置的位置点取土，取回的土样用紫外分光仪测定土壤硝态氮、铵态氮含量。

2 滴灌施肥土壤水氮运移基本方程及模型

2.1 土壤水分运动基本方程

根据达西定律和质量守恒定律[16]及土壤水动力学原理，假定土壤均质和各向同性，忽略空气、温度和滞后效应对水分运动的影响，滴灌三维土壤水分入渗过程可用Richard方程来描述：

(1)

式中,φm为基质势(cm)；k(θ)为非饱和土壤导水率(cm·h-1)；θ为土壤体积含水率(cm3·cm-3)；t为时间(h)。

2.2 土壤氮素运移基本方程

尿素态氮运移方程

(2)

铵态氮运移方程

+k0ρ+k1θC1-k2θC2-k2ρs

(3)

硝态氮运移方程

+k2θC2+k2ρs-k3θC3

(4)

式中,x为空间坐标(角标i=1，2，3，j=1，2，3，表示空间3个方向)；C1、C2、C3分别为尿素态氮、铵态氮、硝态氮的质量浓度(mg·cm-3)；s为固相氨态氮吸附量，s=kdC2，kd为铵态氮吸附速率(cm3·g-1)；qi为土壤水通量(cm·h-1)；k0为有机质矿化速率(h-1)；k1为尿素水解速率(h-1)；k2为铵态氮硝化速率(h-1)；k3为硝态氮反硝化速率(h-1)；ρ为土壤干容重(mg·cm-3)；Dij为水动力弥散系数张量的分量(cm2·h-1)，包含纵向、横向弥散系数DL和Dr(cm)，以及自由水中的分子扩散系数DW(cm2·h-1)。运移及转化参数取值见表2。基于以上土壤水、溶质动力学原理，利用HYDRUS软件建立几何模型并进行数值计算。

2.3 土壤水力参数

本文土壤水分特征曲线θ(h)利用高速离心机测定，土壤饱和导水率K(h)用定水头法进行测定，土壤水分特征曲线采用Van Genuchten模型拟合[17-18]。

(5)

(6)

其中

Se=(θ-θr)/(θs-θr)

(7)

m=1-1/n,n>1

(8)

式中,θs为土壤饱和导水率(cm3·cm-3)；θr为土壤残余含水率(cm3·cm-3)；α,n,m为土壤水分特征曲线实测值(土壤含水率、土壤水吸力)拟合得到的参数(cm-1)；h为土壤负压水头(cm)；Ks为渗透系数(cm·min-1)；Se为土壤有效含水率(cm3·cm-3)；l为空隙连通性参数，一般情况取为0.5。拟合得到土壤水分特性的参数值为θr=0.03 cm3·cm-3；θs=0.48 cm3·cm-3；α=0.036 cm-1；n=1.56；Ks=0.125 cm·min-1。

2.4 土壤溶质运移参数

模型中氮素运移转化参数的确定根据国内外文献[19-23]，并根据实测试验数据以流量分别为2、3 L·h-1对各参数进行率定确定，模型中氮素转化参数的取值见表2。

2.5 定解条件

2.5.1 初始条件 假定土壤初始含水率和氮浓度在研究区域内均匀分布，初始条件可表示为：

θ(x,y,z,t)=θ0
(0≤x≤X，0≤y≤Y，0≤z≤Z，t=0)

(9)

C(x,y,z,t)=C0
(0≤x≤X，0≤y≤Y，0≤z≤Z，t=0)

(10)

式中,θ和C分别为土壤含水率，土壤铵态氮或硝态氮质量浓度(mg·cm-3)；θ0和C0分别为初始土壤含水率,初始土壤铵态氮或硝态氮质量浓度(mg·cm-3)；t为时间(h)，初始时刻时间为0 h。X，Y和Z为模拟区域边界(装置物理边界)的坐标。

2.5.2 边界条件 根据试验设计模拟研究区域为圆柱体(高50 cm，半径35 cm)，滴头在圆柱轴心的正上方，取圆柱体的一半为模拟研究区域。根据试验监测，滴灌开始时，土壤饱和区发展很快，随着时间的推移，饱和区扩展的速度逐渐减小，一定时间后饱和区半径(Rs)趋于定值为7.3 cm，所以模拟时忽略饱和区半径随时间的变化。饱和区在滴灌过程中为不随时间变化的定流量边界，不灌水时为零通量边界。上边界其余覆膜区域为零通量边界。故模拟区域内将上边界分为两种情况：

(1)饱和区流量和溶质浓度边界条件为：

(11)

(12)

(2)非饱和区流量和溶质浓度边界条件为：

(13)

(14)

式中,Ca为肥料溶液尿素态氮的质量浓度(mg·cm-3)；T为灌水历时(h)；C1，xi，Dij,qi与上式(4)中的含义相同；K(θ)为非饱和土壤导水率(cm·h-1)；σ(t)为灌水期间滴头处定流量边界通量(cm·h-1)；t>T时，σ(t)=0，z为竖轴方向的坐标。

模拟计算区域边壁为零通量边界；下边界为自由排水边界。

2.6 数值求解

利用HYDRUS软件对上述模型进行求解，用Galerkin有限单元法对数学模型进行数值计算[24]；时间采用隐式差分法离散[25]，空间离散采用伽辽金有限单元法。计算区内网格划分采用三棱柱单元，滴头饱和区半径内适当加密网格，以提高计算精度。根据试验设置，在模拟区设置与试验布置相对应的12个观测点，以便于对模拟结果进行验证，以滴头正下方土表为坐标原点，观测点坐标分别为(水平，纵向，垂向)：点1(5,0,-5)，点2(15,0,-5)，点3(25,0,-5)，点4(35,0,-5)，点5(5,0,-15)，点6(15,0,-15)，点7(25,0,-15)，点8(5,0,-25)，点9(15,0,-25)，点10(25,0,-25)，点11(5,0,-35)，点12(15,0,-35)。模拟的网格划分及观测点设置如图2所示半圆柱体。

表2 土壤尿素态氮运移转化参数

3 模型验证及结果分析

3.1 土壤含水率分布规律

土壤含水率的分布对作物的生长有直接影响，并影响土壤溶质等因素的变化，基于试验研究，用数值模拟的方法对滴灌及其再分布过程土壤含水率分布进行分析，以流量4 L·h-1土体剖面内各观测点土壤含水率对模型进行验证，对比了滴头流量为4 L·h-1时各观测点的模拟值与实测值。模拟与实测的土壤含水率分布用Kriging插值法进行插值计算，在Surfer软件中绘制。由图3可知不同时刻模拟与实测值的土壤含水率剖面分布基本一致。滴灌结束时，土壤含水率实测均值与模拟值在土层5 cm相对误差为4%，土层15 cm相对误差为10%，土层25 cm相对误差为7%。再分布过程，120 h部分观测点的相对误差较大。

注:图中数字所示为观测点的分布位置。Note: The numbers in the figure show the distribution of observation points.图2 模拟研究区域网格划分及观测点布置Fig.2 Flow domain gridding and layout of observation points in simulation region

注：图中土壤含水率单位为cm3·cm-3。Note: The unit of soil moisture in the figure is cm3·cm-3.图3 土壤含水率分布剖面模拟与实测对比Fig.3 Comparison of simulation and measured values of soil water content

在湿润土体内，滴灌结束时土壤含水率的分布随着离滴头距离的增加而减小；经再分布72 h，土体内水分向下层运移，土层25～30 cm的土壤含水率增大到0.2 cm3·cm-3。120 h土体剖面内的土壤含水率均有所下降，较滴灌结束时下降了18%，且土壤剖面内的土壤含水率分布趋于均匀稳定。对滴灌后土壤含水率的再分布过程模拟与实测值均能较好地吻合，各观测点相对误差的绝对值均在10%以内，模拟能较好地反映土壤含水率的分布规律，可用数值模拟的方法对不同滴灌系统条件下的土壤含水率分布情况进行预测。

3.2 土壤氮分布特征

氮肥是植物生长所必需的营养元素，有机态氮不能被植物直接吸收利用，铵态氮、硝态氮是植物能直接吸收利用的营养离子[27-28];尿素是目前广泛使用的氮肥，在土壤中能直接水解转化为铵态氮、硝态氮，研究滴灌施肥后土壤硝态氮、铵态氮量的时空分布有助于确定不同作物氮肥的次施用量及浓度。

3.2.1 土壤铵态氮时空分布特征 土壤铵态氮在土体中易被吸附，在土壤剖面中的分布决定于作物生长、灌溉方式、土壤性质等因素[29]。基于试验与模拟，以滴头流量4 L·h-1的土壤铵态氮动态分布进行模拟值与实测值的对比分析，模型在与试验布置相同位置坐标点处分别设置12个观测点，根据观测点的布置选取了有代表性的3个观测点(点1、6、8)进行分析(观测点位置见图2)。

由图4可知不同观测点由于硝化作用土壤铵态氮质量浓度均随时间的推移逐渐减小，且模拟值与实测值较接近。对滴灌结束后24～120 h的土壤铵态氮质量浓度模拟值与实测值进行对比分析，距滴头较近的观测点1(5，-5)土壤铵态氮质量浓度在土体中较大，72 h之前随时间的变化幅度较大，随后变化较小，其24、48、72、96、120 h各模拟值与实测值的相对误差分别为17%、20%、20%、17%、14%；观测点6(15,-15)、点8(5,-25)的铵态氮质量浓度在土体中随时间的变化幅度较小，24～120 h内模拟值与实测值的相对误差变化范围分别为14%～20%，12%～20%。观测点5(5,-15)的相对误差变化范围为1%～14%，对其它各观测点进行统计分析其相对误差最小值为2%，均在20%以内。由于实际土壤中微生物的作用不同，且土壤环境复杂，同一观测点处不同时间土壤铵态氮质量浓度模拟值与实测值相对误差存在差异。

图5见16页中土壤铵态氮质量浓度空间分布表明：滴灌后土壤铵态氮质量浓度主要分布于距滴头20 cm的范围内，随着离滴头距离的增加而逐渐减小，24 h土壤铵态氮质量浓度最大，到120 h时减少了40%。滴灌后不同观测点的铵态氮浓度不同，且硝化作用存在差异，120 h部分观测点的土壤铵态氮质量浓度低于本底值。模拟值与试验值的相对误差在合理范围内，土壤铵态氮分布结果与Breve[23]模拟的结果相近，可用该模型模拟土壤铵态氮的动态分布。

图4 不同观测点滴灌后24～120 h土壤铵态氮质量浓度变化Fig.4 The change process of concentration at different observation points 24～120 h after drip irrigation

3.2.2 土壤硝态氮时空分布特征 滴灌施肥后土壤硝态氮分布因灌溉方式、施肥种类等多种因素而存在差异。结合试验，利用数值模型对滴灌施尿素后土壤硝态氮的时空分布进行模拟，并对比分析实测值与模拟值。由图6可知，不同观测点土壤硝态氮质量浓度均随着时间的推移而逐渐增大，且模拟值与实测值较吻合。观测点1(5，-5)土壤硝态氮质量浓度随时间的增长较快，72 h之后变化较小，24、48、72、96 h和120 h内土壤硝态氮质量浓度模拟值与实测值的最大相对误差分别为18%、19%、16%、19%、20%；观测点6(15,-15)、点8(5,-25)随时间的变化幅度较小，24～120 h土壤硝态氮模拟值与实测值的相对误差变化范围分别为1%～20%，19%～20%。观测点5(5,-15)模拟值与实测值的相对误差变化范围为14%～20%，对其它各观测点进行统计分析其相对误差最小值为13%，均在20%以内。

24 h时观测点1，3，6，8，5的土壤硝态氮质量浓度分别为：0.437，0.467，0.451，0.482 mg·cm-3和0.447 mg·cm-3，相互间的差值小于0.05 mg·cm-3，各观测点在空间分布上差异不大，48 h之后土体中观测点1(5，-5)的土壤硝态氮质量浓度较大，120 h时为1.2 mg·cm-3。土壤铵态氮在硝化作用下逐渐转化为土壤硝态氮，由于滴头附近的土壤铵态氮质量浓度较大，土壤硝态氮质量浓度在滴头附近增加较多。由图6和图7可知：滴灌施肥后土壤硝态氮在整个土体内均有分布，48 h后土壤硝态氮质量浓度在离滴头20 cm的范围值较大，且空间分布上随距滴头距离的增加而减小，垂直方向上从距滴头5 cm的观测点1(5，-5)到距滴头25 cm的观测点8(5，-25)减少了53%。该空间分布特征与Breve[23]模拟的结果相近。土壤硝态氮的模拟值与实测值不同观测点、不同时间其相对误差存在差异，但其平均相对误差在18%以内，基本反映出土壤硝态氮分布的动态变化过程，可用该模型模拟滴灌施肥土壤硝态氮的动态分布。

图6 不同观测点滴灌后24～120 h土壤硝态氮质量浓度变化Fig.6 The change process of soil nitrate nitrogen concentration at different observation points 24～120 h after drip irrigation

图5 滴灌结束后72 h、120 h土壤铵态氮质量浓度分布Fig.5 The concentration distribution of in soil 72 h and 120 h after drip irrigation

图7 滴灌结束后72 h、120 h土壤硝态氮质量浓度分布Fig.7 The concentration distribution of nitrate nitrogen in soil 72 h and 120 h after drip irrigation

4 结论与讨论

尿素施入土壤后，会以分子形态被土壤胶体吸附, 但在合适的水热条件和丰富的有机物料供应情况下，会很快水解[30]。李世清[31]研究表明温度、土壤含水率、尿素起始浓度、土壤质地、外加脲酶对尿素水解速率均有影响；汪建飞[32]研究了培养时间、尿素浓度、土壤温度、土壤含水量对土壤外源尿素(施入土壤中的尿素)水解速率的影响，其影响大小顺序为：培养时间>尿素浓度>土壤温度>土壤含水量。樊庆鲁[33]指出在pH>5.6 的土壤中, 温度为 25℃时，尿素3 d可完全分解。王旭洋[34]研究了滴灌施肥(尿素)条件下速效氮在土壤中的运移，表明滴灌结束1 d，尿素水解，铵态氮含量达到最大值，随后随时间的推移逐渐减小；候红雨[35]研究表明尿素施入土壤6～24 h内，剖面中尿素的含量普遍降低了约66.7%。以上研究表明施入土壤中的尿素在不同的条件下，水解转化为铵态氮和硝态氮的速率不同。

本文采用试验与数值模拟的方法对滴灌施肥后土壤水肥的运移变化进行研究，得出滴灌结束24 h，尿素水解，铵态氮含量达到最大值。数值模拟中根据试验设计设置参数，并对参数值率定，对比分析了模拟与试验结果，模型能够反映土壤水氮运移转化过程，模拟所得土壤含水率剖面分布与实测值基本一致，各观测点的土壤硝态氮、铵态氮质量浓度的模拟值与实测值相对误差均在20%以内，不同观测点、不同时间的相对误差均不同。滴灌结束时，土壤含水率随离滴头距离的增加而减小。滴灌施肥后土壤铵态氮质量浓度主要分布于离滴头20 cm的范围，且随着时间的推移逐渐减小；24～120 h土壤硝态氮质量浓度各观测点随着时间的推移逐渐增大，且硝态氮质量浓度在滴头20 cm的范围内由0.442 mg·cm-3增加到1.2 mg·cm-3，48 h后土体内土壤硝态氮质量浓度空间分布上随着离滴头距离的增加呈减小的趋势。由于实际土壤中尿素态氮的转化受多种因素的影响，且文中数值模型进行了简化假设，忽略了铵态氮的挥发，土壤氮素的模拟值与实测值存在一定误差，该模型还有待进一步研究完善。