基于HYDRUS-2D的滨海地区膜下滴灌土壤水盐运移模拟研究

马 波，周青云，张宝忠，王 航，尹林萍

(1.天津农学院水利工程学院，天津 300384；2.中国水利水电科学研究院流域水循环模拟与调控国家重点实验室，北京 100083；3.国家节水灌溉北京工程技术研究中心，北京 100083)

土地盐碱化是全球农业生产面临的一个严峻问题，不仅影响作物正常生长，还会造成大量土地退化，粮食减产[1]。我国约有920.9×104hm2土地受到不同程度盐碱化危害[2]，天津滨海地区土壤盐分主要来自海水，地势低洼、排水不畅、地下水埋深浅、矿化度高、水资源开发及利用不合理等因素导致该地区土壤盐碱化问题严重[3]。膜下滴灌是将覆膜技术与滴灌技术结合起来开发利用盐碱地的灌溉技术，自推广以来取得显著效果[4]。覆膜不仅可以有效减少植物棵间蒸发，实现植物根区土壤的增温保墒，而且能够减少土壤积盐[5]，改善植物生长环境，提高产量。滴灌具有小流量、高频率、能够直接将水分输送到植物根部进行灌溉的特点，较传统灌溉水分损失减少，利用滴头在植物根部形成的淡化湿润区实现精确灌溉，大大提高了农业用水效率，并且能够有效减少深层渗漏和次生盐渍化问题发生[6]。

水是土壤脱盐的原动力，也是盐分运动的载体[7]，不合理的灌溉方式以及灌水定额都会影响土壤盐分的分布与运移。国内外许多学者对土壤水盐运移规律进行了深入研究，郑晓辉等[8]对新疆干旱地区盐碱地不同灌水方式下水盐运移的研究发现，相比于沟灌、漫灌，滴灌是最节水的一种灌水方式，并且局部压盐效果最好。牟洪臣等[9]基于田间实测数据，分析新疆干旱地区棉田土壤水盐运移规律得出，水平方向上滴头位置相较于膜边盐分累积较少，竖直方向上0～20 cm土层淋洗效果明显，而20～60 cm土层处于缓慢积盐状态，60 cm以下土层盐分累积严重。汪志荣等[10]对不同质地土壤滴灌研究发现，土壤湿润体体积与滴头流量、灌水量及土壤性质有关。利用室内试验，分析不同灌水定额1.4 L、2.8 L、4.2 L对土壤盐分淋洗效果，李远等[11]发现灌水量越大，淋洗效果越好，湿润峰位置越低。HYDRUS模型已被广泛运用于模拟土壤水盐运移过程。李显溦等[12]对新疆地区暗管排水、排盐条件的土壤水盐运移进行模拟研究发现，模拟值与实测值具较高精度。Chen等[13]利用HYDRUS-2D模型对咸水滴灌下土壤水盐运移进行模拟研究，验证了该模型用于田间试验的可靠性。陆阳等[14]对地下水浅埋区土壤水盐运移规律研究发现，滴灌可以有效减少宁夏银北灌区耕作层土壤积盐，控制地下水位。

上述研究主要针对西北干旱地区不同灌水量、灌水方式土壤水盐运移进行研究，且多偏重于微咸水灌溉及暗管排水研究，而针对北方滨海地区不同灌水定额的土壤水盐运移规律研究较少。本文针对滨海地区淤泥质盐渍土[15]粘度大、地下水埋深较浅、矿化度较高等特点，通过田间试验研究不同灌水定额下的土壤水盐运移变化，并利用HYDRUS-2D模型对土壤水盐运移进行模拟，揭示滨海地下水浅埋区土壤水盐运移规律，以期为滨海地区制定合理的灌溉制度，以及盐碱地预防和治理提供一定理论依据。

1 材料与方法

1.1 试验区概况

试验田位于天津市东南部的葛沽镇(117°38′E，38°98′N)，该地区属于海河流域下游海积与河流冲积形成的平原区，平均海拔为3 m，地下水埋深约为1 m，矿化度为5 g·L-1，土壤为淤泥质盐渍土，pH值为8.4。该区属暖温带季风型大陆性气候，四季分明，夏季炎热多雨，七月份温度达到最高，平均气温26～28℃，多年平均年降雨量为556.4 mm，6、7、8月降雨较多，占全年降雨的75%，多年平均蒸发量为1 809.6 mm[16]。

1.2 试验设计

试验以春玉米(郑丹958)为研究对象，试验地面积50.0 m9.3 m，设置2个不同灌水定额处理(FI20和FI10)，每个处理设置3个重复，每个重复小区面积为8.3 m×9.3 m,其中FI20处理单次灌水定额为20 mm，FI10处理单次灌水定额为10 mm。种植方式为一膜两管两行等间距平作，行距60 cm，株距30 cm。灌溉方式为膜下滴灌，滴灌带间距为60 cm，滴头间距为30 cm，滴头布置在玉米茎秆处，滴头流量1.38 L·h-1。覆膜方式为半覆膜，膜宽80 cm。在行距方向上布置2个取样点，分别位于距离玉米茎秆15 cm和30 cm处，分别记做1号取样点和2号取样点，共12个取样点，试验地布置如图1所示。2017年和2018年试验期间分别进行了3次灌水，灌水日期分别是2017年6月10日、6月15日和7月10日，2018年6月1日、6月22日和7月6日，试验土样的采集在灌水后的1～2 d。单点滴灌属于三维水分运动问题，垂直平面上，土壤水分沿滴头所在平面两侧呈对称分布，行距方向的两滴头中间位置形成零通量面，故可以将滴灌问题简化成中心对称的二维水分运动问题[17]。试验地设置暗管排水设施，可将地下水位控制在100 cm处，故取玉米左右行距一半，选择60 cm100 cm典型区域进行模拟计算，如图2所示。

1.3 观测项目与方法

1.3.1 土壤水分、盐分及土壤水分特征曲线 土壤含水率：将60 cm深度土层划分为0～20、20～40、40～60 cm 3层，利用试验前埋在土壤中的PR2管原位测定玉米根区各层土壤含水率，并用烘干法对土壤含水率进行率定。

土壤含盐量：土壤样品自然风干，过2 mm筛，将处理后土壤样品按土∶水比例1∶5混合，振荡、过滤，利用DDS-307型电导率仪测定浸提液电导率，换算得土壤含盐量。

图1 试验地布置图Fig.1 Diagram of test layout

图2 模拟区域Fig.2 Simulation area

土壤水分特征曲线：挖土壤剖面，用环刀取0～20 cm、20～40 cm、40～60 cm、60～80 cm、80～100 cm土壤，通过日系R11D2高速恒温冷冻离心机测定土壤吸力与含水率的关系，绘制土壤水分特征曲线。

1.3.2 玉米根系分布情况 玉米根系分布密度函数通过LA-S植物根系分析系统测定。

1.3.3 气象因子与参考作物蒸发蒸腾量 试验地设有自动气象站，降雨量直接来自气象站数据，参考作物蒸发蒸腾量利用FAO56-PM公式计算得到[18]。2018年降雨和参考作物蒸发蒸腾量如图3所示。

(1)

式中,ET0为参考作物蒸发蒸腾量(mm·d-1)；Rn为输入冠层净辐射量(MJ·m2d-1)；G为土壤热流强度；es、ea为饱和与实际水汽压(kPa)；Δ为饱和水汽压与气温关系斜率(kPa·h-1)；γ为干湿计常数(kPa·℃-1)；u2为2 m高处风速(m·s-1)；T为平均气温(℃)。

1.4 数据处理

运用Excel进行数据分析，Excel以及Origin 9.0画图，SPSS 24.0对模拟值和实测值吻合程度进行分析。

2 HYDRUS-2D模型

图3 2018年降雨和参考作物蒸发蒸腾量变化Fig.3 Variation in rainfall and reference cropevapotranspiration in 2018

2.1 数学模型构建

2.1.1 土壤水分运动模型 假设试验地土壤均质、各向同性，忽略气相和温度对水分运动的影响，考虑质量守恒定律和达西定律，采用Richards水分运动控制方程[20]，其公式为

(2)

式中,θ为土壤体积含水率(cm3·cm-3)；h为土壤基质势(cm)；z为垂向坐标，向上为正；x为横向坐标；t为时间(d)；K(h)为非饱和土壤导水率(cm·d-1)；S(h)为根系吸水项。

2.1.2 土壤溶质运移模型 溶质在土壤中的运移受对流和水动力弥散作用影响，模型采用对流—弥散方程来描述溶质的运移[20]，其公式为：

(3)

式中,c为土壤溶液浓度(g·cm-3)；qi为水流通量(cm·d-1)；DT为横向弥散系数(cm2·d-1)；DL为纵向弥散系数(cm2·d-1)。

土壤含盐量与土壤溶液中盐分含量转换关系式为[21]：

(4)

式中,SSC为土壤含盐量(g·kg-1);γ为土壤堆积密度(g·cm-3)。

2.1.3 根系吸水模型 根系吸水采用Feddes模型[22]，计算公式为:

S(h)=α(h)b(x,z)TpL

(5)

(6)

(7)

式中,α(h)为水分胁迫系数；b(x,z)为根系吸水分配密度函数，按照实际根系分布确定；Tp为潜在蒸腾速率(cm·d-1)；L为根区分布最大土壤表面宽度(cm)；h1、h2、h3和h4为根系从土壤中吸水的不同压力水头(cm)，具体参数参考Wesseling等[23]提出的参考值。

2.2 初始条件与边界条件

模拟区域上边界为半覆膜，滴灌过程中滴头处为饱和水头边界；覆膜部分为零通量边界，未覆膜部分为大气边界。模拟区域地下水埋深为100 cm，故下边界为饱和水头边界，左右边界两侧为对称面，为零通量边界。

(1)初始条件。即模拟区各位置初始的土壤水分盐分分布情况。

h(x,z,t)=h0(x,z),
(0≤x≤60,0≤z≤100,t=0)

(8)

式中,h0为初始土壤基质势(cm)。

c(x,z,t)=c0(x,z),
(0≤x≤60,0≤z≤100,t=0)

(9)

式中,c0为初始含盐量(g·kg-1)。

(2)上边界条件。滴灌过程中滴头处为饱和水头边界，故

h(x,z,t)=0, (z=0,t≥0)

(10)

(11)

式中,q为地表水分通量(L·h-1)，cs为上边界盐分浓度(g·kg-1)。

覆膜部分上边界条件为零通量条件，故

(12)

(13)

未覆膜部分上边界条件为大气边界条件，故

(14)

(15)

(3)下边界条件。下边界为饱和水头边界，故

h(x,z,t)=0, (z=100,t≥0)

(16)

c(x,z,t)=cb(t), (z=100,t≥0)

(17)

式中,cb为地下水含盐量(g·kg-1)。

2.3 模型参数

将试验测得土壤理化性质与HYDRUS-2D中的Rosetta模型相结合，得出土壤水力参数；将2017年土壤含水率和含盐量实测值输入HYDRUS-2D模型反求溶质运移参数。利用2017年土壤水分盐分实测数据对模型参数进行率定，土壤含水率R2为0.664，均方根误差(RMSE)为0.025 cm3·cm-3，平均绝对误差(MAE)为0.02 cm3·cm-3，土壤含盐量R2为0.718，均方根误差(RMSE)为0.34 g·kg-1，平均绝对误差(MAE)为0.25 g·kg-1，率定精度较高，故最终确定土壤水力参数如表1所示，土壤溶质运移参数如表2所示。

2.4 模型模拟

对2018年6月1日—2018年7月21日共计50 d春玉米土壤水盐运移进行模拟，模拟区域如图2所示。沿竖直方向将100 cm土层根据土壤质地及容重分成3层，共划分101个节点，水平方向划分61个节点，采用三角形网格将模拟区域离散化，共生成12 000个网格，模拟的时间单位为d。

2.5 模型统计量

利用SPSS 24.0对土壤含水率、含盐量模拟值与实测值进行T检验，并计算其均方根误差(RMSE)、平均绝对误差(MAE)来验证模型用于大田模拟的可靠性[24]。

(18)

(19)

式中,Mi为第i个实测值；Si为第i个模拟值；n为数据个数。RMSE和MAE越接近0，表明模拟值与实测值吻合程度越好。

3 结果与分析

3.1 距离滴头不同位置土壤水盐运移过程

以2018年6月1日灌水前土壤含水率和含盐量为初始值对模拟区域进行土壤水盐模拟计算，滴灌前后FI20处理土壤水盐运移情况如图4、图5所示，可以看出整个模拟期内，灌水前后距离滴头不同位置取样点土壤水盐运移趋势基本一致。6月1日灌水前0～20 cm土层土壤含水率较低，最小值出现在滴头正下方，为0.15 cm3·cm-3，含盐量较高，为4.85 g·kg-1。

盐随水动，土壤水分运动决定着盐分的运移。在灌溉过程中，土壤盐分受到水分淋洗作用向外扩散，土壤含盐量在滴头处出现最小值。水平方向上，随着距滴头距离的增加土壤含盐量呈现增大趋势，土壤水分将盐分带到膜边，并在膜边形成积盐区，距离滴头更近的1号取样点洗盐效果明显好于2号取样点；竖直方向上，随着土层深度增加盐分逐渐增大，0～20 cm土层盐分明显降低，40～60 cm土层盐分降低不明显，土壤盐分随着水分向下运移，下层土壤出现盐分累积现象，洗盐区域为0～60 cm。

在水分重分布过程中，随着玉米生长、作物蒸腾耗水，土壤下层的盐分随水分向上运动，水去盐留，土壤表层出现盐分积聚现象[25]，使得6月22日灌水前表层土壤含盐量较高。6月22日与7月6日进行灌水处理，土壤受灌溉水淋洗，0～20 cm土层盐分明显减小，40～60 cm土层含盐量增大，灌溉水将盐分淋洗到60 cm土层以下。80～100 cm土层土壤饱和导水率(Ks)较小，土壤干容重(γd)较大，且试验地设置有地下排水暗管能够将地下水位控制在100 cm以下，使得地下水对0～60 cm土层土壤含水率和含盐量影响很小。

表1 土壤水力参数

表2 土壤溶质运移参数

注：图中线条表示土壤体积含水率等值线。Note: The line in the figure represents the contour of soil volume water content.图4 FI20处理灌水前后土壤体积含水率变化Fig.4 Variation of soil volumetric water content before and after irrigation of FI20 treatment

注：图中线条表示土壤含盐量等值线。Note: The lines in the figure represent the contour lines of soil salt content.图5 FI20处理灌水前后土壤含盐量变化Fig.5 Variation of soil salt content before and after irrigation of FI20 treatment

3.2 不同灌水量土壤盐分淋洗效果

图6、图7为FI20处理和FI10处理各取样点土壤含盐量变化，对比图6、图7可知不同灌水定额下各层土壤含盐量变化趋势基本一致，0～20 cm土层淋洗效果明显好于20～60 cm土层。2018年6月1日、6月22日和7月6日灌水后0～20 cm土层土壤含盐量下降幅度最大，其土壤含盐量比灌水前降低16.1%，随着土层深度的增加含盐量降低幅度有所减小，20～60 cm土层土壤含盐量比灌水前降低约9.6%，并且下层土壤含水率和含盐量变化幅度维持在较小范围内。结果表明滴灌对0～20 cm土层盐分淋洗效果明显，对20～60 cm土层盐分有抑制作用。

图8为1号取样点和2号取样点各深度土壤含盐量变化，对比1、2号取样点各深度土壤含盐量可知不同深度高灌水定额FI20处理土壤含盐量均小于低灌水定额FI10处理，在0～20 cm、20～40 cm、40～60 cm土层中，FI20处理对土壤盐分的抑制作用分别比FI10处理高26%、11%、19%。从图6、图7可以看出，整个模拟期内FI20处理40～60 cm土层土壤含盐量相对较低，而FI10处理40～60 cm土层土壤含盐量相对较高，说明高灌水定额FI20处理将土壤盐分淋洗到60 cm土层以下，而低灌水定额的FI10处理未能将盐分淋洗到60 cm以下，在40～60 cm土层出现盐分累积现象，结果表明高灌水定额FI20处理对土壤盐分淋洗效果优于低灌水定额FI10处理。

图6 FI20处理各取样点土壤含水率和含盐量变化Fig.6 Variation of soil water content and salt content at each sampling point in FI20 treatment

图7 FI10处理各取样点土壤含水率和含盐量变化Fig.7 Variation of soil water content and salt content at each sampling point in FI10 treatment

图8 1号取样点和2号取样点各深度土壤含盐量变化Fig.8 Variation of soil salt content of sampling points 1 and 2 at different soil depths

3.3 模型可靠性验证

利用2018年土壤水分、盐分实测数据对模型可靠性进行验证，由图6、图7可以看出土壤含水率和含盐量模拟值与实测值个别点存在一定的偏差，但总体上拟合效果较好。土壤含水率和含盐量T检验的p值分别为0.592和0.752，均大于0.05，RMSE与MAE值分别为0.022、0.014 cm3·cm-3和0.357、0.288 g·kg-1，相对较小，与Phogat等[26]的研究统计结果范围一致，如表3所示，表明模型可靠，能够准确反映滴灌前后土壤水盐运动规律，可用于田间土壤水盐运移模拟。

表3 HYDRUS-2D模型统计量

4 讨 论

对滨海地下水浅埋区不同灌水定额土壤水盐运移进行模拟研究，结果表明，滴灌使土壤水分含量发生变化，从而导致盐分再分布，土壤水分将表层盐分带到下层，使表层土壤含盐量明显下降，而下层土壤出现盐分累积现象；植物蒸腾耗水使得土壤盐分随水分向上运动并在地表聚集。

距离滴头的深度及灌水量的大小都会影响土壤盐分的分布，灌溉水将土壤表层盐分淋洗到下层，0～20 cm土层土壤淋洗效果明显，20～60 cm土层淋洗效果一般；距滴头距离越近淋洗效果越好，灌溉水形成的湿润区将玉米根区盐分带到膜边，与牟洪臣等[9]、张治[27]的研究结论一致。水分重分布过程中，由于覆膜原因，降雨少量进入土壤，地表蒸发较小，土壤含水率的变化主要取决于植物蒸腾耗水量，植物蒸腾作用与根系分布密度呈正相关[28]，玉米根系主要吸水部分分布在0～20 cm土层[29]，所以在土壤水分重分布过程中0～20 cm土层相比其他土层含水率和含盐量变化幅度更大，这与周和平等[30]研究结论一致。高灌水定额的FI20处理将盐分淋洗到60 cm土层以下，而FI10处理未能将土壤盐分淋洗到60 cm土层以下，代表了滴灌淋洗深度，这与张治等[27]研究结论一致。

HYDRUS-2D模型能够准确模拟滨海地区土壤水盐运移情况，灌水量的多少决定淋洗效果，谭军利等[31]认为盐碱地滴灌灌水量较小会影响压盐效果。对于地下水埋深较浅的干旱地区，朱文东等[32]认为大量灌水提高地下水位，棵间蒸发和根系吸水作用使地下水向上运动，加剧土地盐碱化，这使得灌水量的多少成为滨海地下水浅埋区研究的关键。本试验采用的灌溉方式为小灌水定额滴灌，对土壤中盐分有较好的抑制作用且没有抬高地下水位，试验地80～100 cm土层是粘土，土壤致密，并且设置有排水暗管能够将地下水位控制在100 cm处，使得高矿化度的地下水对0～60 cm土层土壤含盐量影响较小。

5 结 论

1)灌溉过程中，不同取样点不同深度土壤含水率和含盐量沿水平、竖直方向变化趋势基本一致。滴灌使土壤盐分发生重分布，土壤水分将表层盐分带到下层，使土壤表层含盐量明显下降，而土壤下层出现盐分累积现象；作物生长、蒸腾耗水使得土壤盐分向上运动并在地表聚集。0～20 cm土层较其他土层淋洗效果明显，灌水对下层土壤影响较小；距滴头距离越近淋洗效果越好。水分重分布过程中，0～20 cm土层相比其他土层土壤含水率和含盐量变化幅度最大。这对于滨海地区盐碱地预防与治理具有一定的指导意义。

2)灌水定额对土壤盐分分布影响显著，高灌水定额FI20处理对土壤盐分的淋洗效果明显好于低灌水定额FI10处理，FI20处理将土壤盐分淋洗到60 cm土层以下，而FI10处理由于灌水量小，未能将土壤盐分淋洗到60 cm土层以下。

3)利用HYDRUS-2D对滨海地区不同灌水定额处理下土壤水盐运移进行模拟，将模拟值与田间实测数据对比分析，结果表明，HYDRUS-2D模型可靠，能够很好模拟滨海地下水浅埋区土壤水盐运移。