盐结皮厚度对土壤水盐分布特征的影响

2023-09-28李新虎王弘超李佳琳

干旱区地理(汉文版) 2023年8期

郭敏，李新虎，王弘超，李佳琳

（1.中国科学院新疆生态与地理研究所荒漠与绿洲生态国家重点实验室，干旱区生态安全与可持续发展重点实验室，新疆乌鲁木齐 830011；2.新疆阿克苏绿洲农田生态系统国家野外科学观测研究站，新疆阿克苏 843017；3.中国科学院大学，北京 100049）

盐结皮是土壤盐渍化的极端表现，对土壤性质以及地表水文过程等具有重大影响［1-2］。干旱区土壤水盐运动是一个复杂的动态过程［3］，并且土壤水盐分布特征会影响土壤物理化学及生物过程［4-5］，进而影响土壤重金属再分布［6］、碳排放［7］等，而在干旱区盐结皮普遍存在，因此研究盐结皮土壤水盐分布特征具有重要的理论意义。

目前国内外学者对于土壤水盐分布的研究主要集中在非盐结皮土壤。Zhang 等［8］对蒸发过程中浅层裸土水盐分布进行室内试验，发现由于水盐迁移不同步，盐分剖面呈波浪形分布。当土壤表面被不同材料覆盖时会显著降低土壤蒸发，减轻土壤盐分的表聚［9-11］。盐结皮对土壤的影响类似于秸秆覆盖［12］，因此盐结皮必然会影响土壤水盐分布特征。然而，由于目前的技术手段对高盐区域土壤含盐量测量的局限性以及盐结皮影响下土壤水盐热耦合作用的复杂性，导致目前仅有少数学者对盐结皮土壤水盐分布特征进行研究，如Gran 等［3］通过室内试验，发现盐结皮影响下不同质地土壤含水量变化曲线有很大差异，盐分在0～4 cm 土层累积，并且越接近土壤表面盐分浓度越高。莫治新等［13-14］研究表明盐结皮土壤含水量均大于无盐结皮土壤，不同盐结皮厚度会影响水分的迁移导致不同土层含水量变化规律不同，并且盐结皮越厚的土壤区域其各层土壤中的盐分含量越高。王晓静等［15］对塔里木沙漠公路防护林带土壤进行分析，发现盐结皮厚度会逐渐增长变厚并趋于稳定，土壤电导率呈现先增加后逐渐降低的趋势。Guglielmini等［16］、Dashtian等［17］通过建立模型预测盐分在土壤中的迁移过程，结果表明，盐分浓度最大值出现在液-汽界面附近。这一液-汽转化面被称为蒸发锋［18-19］，能直接反映土壤水盐分布状况，其变化会影响盐分累积的深度和水分运动状态，是描述蒸发过程中水盐分布的重要指标［20］。然而，当前研究针对盐结皮土壤蒸发锋提出的不同观点，使盐结皮土壤水盐分布特征的研究仍然存在争议。一些学者［3,21］认为，蒸发过程中盐结皮土壤蒸发锋区域较窄（2～4 cm），而另一些学者［19-20］在研究中发现有较宽的蒸发锋区域（8～12 cm）。上述研究结果的差异，可能原因是不同的试验条件，如：盐分浓度［22-23］、盐分类型［3］、不同地下水深度［19］等，导致形成不同厚度的盐结皮，从而对土壤水盐分布特征产生不同的影响，但上述研究均未考虑盐结皮厚度这一重要的物理特征。

在干旱半干旱地区土壤盐渍化程度不同，其土壤含盐量有很大的空间变异性，各种厚度盐结皮广泛分布，而盐分浓度是影响盐结皮形成发育的重要因素［19,23］。研究表明，蒸发是土壤水盐运动的主要驱动力［8］，而不同盐结皮厚度会使土壤蒸发产生较大差异［24-25］，因此盐结皮的厚度必然会影响土壤水盐分布特征，但目前盐结皮厚度对土壤水盐分布特征的作用机理还不明确。因此，本研究设置了4 种不同初始浓度以得到不同厚度的盐结皮，分析其土壤水盐的动态变化，明确不同厚度盐结皮土壤水盐分布特征的差异，阐明盐结皮厚度对土壤水盐分布特征的影响机理，为定量描述干旱区盐渍土壤的地表水文过程提供理论基础。

1 材料与方法

1.1 试验材料

供试土壤质地为砂土，土壤容重为1.65 g·cm-3，饱和含水率为39.67%。采用直径10 cm 高33 cm 的聚氯乙烯（PVC）管作为容器，并包裹隔热材料。将试验采用的砂土洗净，去除杂质并将电导率控制在100 μS·cm-1以下，烘干后过0.5 mm 土壤筛均匀装入PVC 管，在PVC 管底部填充3 cm 高的滤石（直径5 mm）。试验土柱自下而上充分饱和24 h。在土柱正上方38 cm 处放置一盏卤素灯（50 W）用于驱动蒸发，整个试验过程平均环境温度为25±3.5 ℃，湿度为（47.2±13.9）%。

1.2 试验设计与方法

试验设置了4 个处理，分别为0 g·L-1（CK）、10 g·L-1（LC）、150 g·L-1（MC）和250 g·L-1（HC），其中CK 为无盐处理，采用蒸馏水饱和，LC、MC 和HC 采用不同浓度NaCl溶液进行初始饱和，以获得不同厚度的盐结皮。试验蒸发采用自动称重平台（型号WP20，北京时域通科技有限公司）测定，并连接数据采集器（型号i-logger，北京时域通科技有限公司）记录数据。采取部分重复逐步退出的方法［3］，即在选定的时间点将部分土柱进行切割，将盐结皮从土壤表面剥离，10 cm 深度以上每隔2 cm 取土，10～30 cm每隔5 cm 取土，对其含水量及含盐量进行测定。根据之前研究［25］所得不同初始浓度各蒸发阶段的时间，本文选取5个关键节点作为退出时间，其中第一阶段选取特征时间a，第二阶段选取特征时间b和c，第三阶段选取特征时间d和e，不同处理各特征时间如图1 所示。采用烘干法测定土壤含水量，将烘干土样研磨、过1 mm土壤筛，按土水比为1:5配置成溶液，用玻璃棒搅拌、静置，之后使用电导率仪测定其电导率，计算土壤含盐量与电导率之间的关系式如公式（1）所示。采用游标卡尺测量盐结皮厚度，之后将盐结皮加蒸馏水溶解、过滤、烘干后测量盐结皮质量。

图1 各处理逐步退出时间节点的蒸发速率Fig.1 Evaporation rate of each treatment at gradually exit time nodes

式中：St为土壤含盐量（g·kg-1）；EC 为土水比1:5 浸提液电导率（mS·cm-1）。

土壤进气值指饱和土壤脱水过程中开始进入空气时的吸力值，与Van Genuchten模型中进气吸力相关参数和形状系数有关［26］。本研究通过引用土壤进气值来分析不同处理的水力学参数。土壤进气值的计算公式如下：

式中：H为土壤进气值；α、n分别为VG 模型中与进气值有关的经验参数。

受试验条件限制，所有处理蒸发锋的位置变化及下移时间并不能被完全观测，因此选择HYDRUS-1D 模型延长蒸发时间，继续观测土壤水分动态过程，以明确各处理蒸发锋的变化。

1.3 土壤水盐动态模型

HYDRUS-1D模型中土壤水分运动方程：

式中：θ(h)为土壤体积含水率（cm3·cm-3）；h为压力水头（cm）；t为水分运移时间（d）；z为垂直方向坐标变量（cm）；K(h)为非饱和导水率（cm·d-1）。

盐分运移方程：

式中：C为溶质浓度（g·cm-3）；D为弥散系数（cm2·d-1）；zij为垂直方向的变量（cm）；q为水通量（cm·h-1）。

模型高度为30 cm，分为31 个节点，节点间距1 cm，分别在土层2 cm、4 cm、6 cm、8 cm、10 cm、15 cm、20 cm、25 cm 和30 cm 处设观测点。模型模拟时间同试验时间相同，之后延长蒸发时间到120 d，初始时间步长为0.01 d，最小时间步长为0.001 d，最大时间步长为1 d。

通过HYPROP 软件实测土壤水分特征曲线得到各处理土壤水力参数，溶质运移参数基于文献确定［27-28］，之后通过HYDRUS-1D 模型采用实测土壤含水量与含盐量对参数进行修正，最终不同处理土壤水力参数与溶质运移参数如表1所示。本试验中采用同种土壤质地，不同初始盐分浓度饱和土壤，导致各处理土壤含盐量有很大不同，而土壤盐分含量会影响土壤持水性和土壤孔隙［29］，从而使各处理土壤水力参数产生较大差异，随着初始浓度的增加，α逐渐减小，n逐渐增大［30］。

表1 不同处理土壤水力参数Tab.1 Soil hydraulic parameters under different treatments

1.4 模型初始条件和边界条件

对于土壤水分运动，试验开始时土壤初始饱和，初始条件为饱和含水量，上边界设置为变通量边界，将试验实测的蒸发通量作为上边界，下边界为零通量边界。对于溶质运移，不同处理初始条件为初始饱和浓度，上下边界条件均设置为浓度通量边界。

1.5 数据处理

采用Excel 2019和SPSS 25软件进行统计分析，使用Origin 2018 软件绘制图形。模型验证使用决定系数（R2）和均方根误差（RMSE），RMSE 越接近于0，表明模型模拟精度越高，R2越接近1，说明模型可以较好地捕捉到实测值的变动趋势。

式中：N为样本容量；Si为模拟值；Mi为实测值。

2 结果与分析

2.1 不同厚度盐结皮影响下土壤水分动态变化

图2 所示为试验结束时含盐处理土壤盐结皮。含盐处理在不同特征时间盐结皮的质量与厚度存在差异（图3），含盐处理盐结皮随着时间的推移质量逐渐增加，厚度逐渐变厚，土壤含盐量对盐结皮厚度的影响较大［31］。李胜辉等［23］研究发现盐浓度从0.5 mol·L-1增加到2.5 mol·L-1时，盐结皮厚度逐渐增加，本试验研究结果与其相似，在试验结束时LC、MC 和HC 盐结皮厚度分别为4.5 mm、6.6 mm 和7.3 mm。通过分析不同处理土壤含水量随时间的动态变化（图4），MC 和HC 土壤含水量相比CK有显著差异（P＜0.05）。对于本研究定义的表层土壤（0～2 cm），盐结皮越厚土壤含水量越大。在特征时间a，CK、LC、MC 和HC 土壤含水量分别为0.123 cm3·cm-3、0.202 cm3·cm-3、0.329 cm3·cm-3和0.293 cm3·cm-3，这是因为土壤中盐溶液使土壤渗透势发生变化，饱和水汽压减小［32］，抑制水分的向上传输，导致含盐处理含水量大于CK。在特征时间b，HC土壤含水量显著性高于CK（P＜0.05），相比特征时间a 仅减少了0.019 cm3·cm-3，可能原因是此时HC 土壤表面形成一层较薄的盐结皮，抑制了水分向上传输，导致土壤含水量变化较小。在特征时间c，CK 表层土壤含水量降低到0.009 cm3·cm-3，表层土壤基本干燥（土壤含水量小于风干土含水量0.050 cm3·cm-3定义为干燥土［18］），而含盐处理由于盐结皮逐渐变厚导致水分传输阻力的增加，此时LC、MC 和HC 土壤含水量分别为0.096 cm3·cm-3、0.194 cm3·cm-3和0.261 cm3·cm-3。在特征时间d～e阶段，各处理表层土壤含水量仅有小幅变化，试验结束时LC、MC 和HC 土壤含水量均大于CK（0.005 cm3·cm-3），分别为0.008 cm3·cm-3、0.131 cm3·cm-3和0.176 cm3·cm-3，盐结皮越厚，表层土壤含水量越大。

图2 试验结束时含盐处理土壤盐结皮Fig.2 Soil salt crust on saline treated at the end of the test

图3 不同初始浓度盐结皮质量与厚度Fig.3 Mass and thickness of salt crust with different initial concentrations

图4 不同处理各深度土壤含水量随时间的变化Fig.4 Variation of soil water content at different depths with time in different treatments

对于本研究定义的上层土壤（2～10 cm），含盐处理土壤含水量在该层变化幅度较大。在特征时间a～c 阶段，由于表层水分的减少，上层土壤水分开始向表层迁移，HC 由于较低的蒸发驱动力［25］，导致其上层土壤含水量相比其他处理变化较小。在特征时间d～e 阶段，LC 其2～4 cm 土层含水量相比0～2 cm 突然增加，土壤含水量的突变表明蒸发锋移动到了2～4 cm 深度，当蒸发锋发生明显下移，其上方土壤水力联系被切断，此时LC土壤表面盐结皮厚度与质量趋于稳定，这与Nachshon 等［33］研究结果相似。然而MC 和HC 上层土壤含水量均大于0.14 cm3·cm-3，仍与表层土壤保持水力联系，未明确观测到蒸发锋的下移，主要原因是试验前期土壤溶液中较多的盐离子影响了土壤水分的运动状态［11］，并且盐结皮的快速形成阻碍了水分的传输，导致土壤含水量相比CK 变化较小。试验结束时，上层土壤含水量表现为HC＞MC＞LC＞CK，盐结皮越厚，上层土壤含水量越大。

对于下层土壤（10～30 cm），不同处理均表现出随着深度的增加土壤含水量逐渐增加。在试验周期内，含盐处理下层土壤与CK 相比变化较为平缓，盐结皮越厚，下层土壤含水量越大，主要原因是随着时间的推移，盐结皮逐渐生长变厚，导致蒸发驱动力降低，下层土壤水分向上迁移速率减缓。在试验结束时，LC下层土壤含水量相比CK较小，而MC和HC 显著大于CK（P＜0.05），可能原因是在试验后期LC盐结皮较薄，对土壤水分的抑制作用远小于MC和HC，并且盐结皮形成增加了蒸发的交换面积［24］，在一定程度上增加了蒸发速率［20］，这一特殊效应导致LC下层土壤含水量小于CK，而对于HC 和MC，这一特殊效应未对盐结皮的抑制作用产生明显影响，从而导致HC和MC下层土壤含水量远大于CK。

2.2 不同厚度盐结皮影响下土壤盐分动态变化

图5所示为含盐处理各深度土壤含盐量随时间的变化。对于表层土壤，不同处理盐分均在土壤表面累积，这是因为土壤中盐溶液在蒸发驱动下逐渐向土壤表面迁移，随着溶液中水分的蒸发，盐分被滞留在土壤表面，当盐分浓度达到饱和状态时，盐分开始在土壤表面结晶析出。在特征时间a，LC表层土壤含盐量大于初始值（3.28 g·kg-1），而MC和HC表层土壤含盐量小于初始值（29.54 g·kg-1和45.79 g·kg-1），可能原因是在特征时间a，MC 和HC 处理2～8 cm 土壤含水量相比初始值减少，并且2～8 cm 土层盐分均小于表层土壤，表明2～8 cm 土层盐分随着水分向上迁移至土壤表层，由于MC 和HC 较高的盐分浓度其表层土壤很快达到饱和状态，盐分结晶析出形成盐结皮，导致表层盐分低于初始值，而LC 盐结皮形成较慢，较低的盐分浓度使得下层盐分持续向上迁移，导致表层盐分高于初始值。随着试验的进行，LC 表层土壤含盐量的变化表现为先减小之后趋于稳定，而MC 和HC 表层土壤含盐量随着时间的推移逐渐减小。造成这种差异的原因可能是LC 在特征时间a～c 过程中蒸发速率较高［25］，下层盐分向上运动并逐渐形成盐结皮，在特征时间d之后，蒸发速率明显减小，此时水分运动主要以水汽扩散的形式向上运动，表层含盐量趋于稳定，而MC 和HC 表层土壤较为湿润，水分主要以液态水的形式传输，携带盐分向上缓慢迁移并结晶析出，表层含盐量持续稳定减小，试验结束时，MC 和HC 表层土壤含盐量相比初始值分别减小了30.7%和15.9%。

图5 不同处理各深度土壤含盐量随时间的变化Fig.5 Variation of soil salt content at different depths with time in different treatments

对于上层土壤，含盐处理盐分在该层变化幅度较大，且土壤剖面最小含盐量位于该层。在特征时间a～c 阶段，LC 最小含盐量所处深度在2～6 cm，MC和HC最小含盐量在6～10 cm。在特征时间d～e阶段不同处理最小含盐量所处深度趋于稳定。在试验结束时，LC 最小含盐量为0.31 g·kg-1，位于2～4 cm深度，相比初始含盐量减小了90.5%，该处土壤含盐量的突变与含水量的突变相吻合，进一步验证了蒸发锋移动至2～4 cm 深度，此时盐分主要在蒸发锋上方累积。而MC和HC最小含盐量在6～8 cm深度，分别为15.87 g·kg-1和31.08 g·kg-1，相比初始含盐量分别减小了46.3%和32.1%。结果表明，盐结皮越厚，上层土壤盐分随时间的减小幅度越小。

对于下层土壤，土壤含盐量随着深度的增加逐渐增加，盐结皮越厚，土壤含盐量减小幅度越小。MC 和HC 下层土壤含盐量相比LC 具有显著性差异（P＜0.05），LC 盐分峰值位于表层土壤，而MC 和HC盐分峰值位于下层土壤。盐分分布产生较大差异主要原因是LC初始浓度较低，土壤表层盐分达到饱和需要下层盐分的补给，盐结皮形成较为缓慢，对土壤水盐抑制作用较小，下层盐分可持续向上迁移，在试验结束时，下层土壤含盐量相比初始值减小了67.3%以上。而MC 和HC 由于较高的初始浓度，下层少量的盐分向上迁移使上层盐分很快达到饱和状态并结晶析出，盐结皮很快形成并逐渐覆盖整个土壤表面，进而抑制了水盐的向上传输，使下层土壤含盐量变化幅度较小，试验结束时MC 下层25～30 cm和HC下层20～30 cm土壤含盐量相比初始值仅有很小的变化。

盐分均衡可以反映一定区域盐分输入量和输出量之间的关系。本研究中含盐处理输入盐分为初始饱和不同浓度盐溶液，盐分结晶析出形成的盐结皮为输出盐分。在试验结束时，LC 表层、上层和下层土壤含盐量分别为1.04 g·kg-1、0.44 g·kg-1和2.13 g·kg-1，输出盐分为7.26 g·kg-1；MC 表层、上层和下层土壤含盐量分别为5.30 g·kg-1、17.86 g·kg-1和33.20 g·kg-1，输出盐分为24.36 g·kg-1；HC 表层、上层和下层土壤含盐量分别为2.13 g·kg-1、64.00 g·kg-1和100.93 g·kg-1，输出盐分为33.04 g·kg-1。结果表明，初始盐分浓度越高，输出盐分质量越大，由于不同含盐处理蒸发速率的差异，初始浓度越高蒸发速率越小［25］，下层盐分向表层迁移速率越小，导致输出的盐分占总盐分质量的比例越小。

2.3 不同处理土壤剖面水盐动态变化模拟

通过HYDRUS-1D模型对不同处理含水量随时间的变化进行了模拟，模拟结果如图6，模拟精度见表2。由表2 可知，CK 在特征时间b 土壤含水量的模拟精度最高，R2和RMSE 分别为0.99和0.006 cm3·cm-3。LC 在特征时间c 土壤含水量模拟精度最高，R2和RMSE分别为0.99和0.006 cm3·cm-3。MC和HC土壤含水量模拟值随着时间的推移模拟精度逐渐变高，在特征时间e，R2分别为0.98 和0.99，RMSE 均为0.008 cm3·cm-3。从图6 可以看出随着初始浓度的增加，HYDRUS-1D 模型对试验初期土壤含水量模拟精度较小。土壤含水量的动态变化采用VG 模型进行模拟，通过分析VG 模型参数发现，不同初始浓度对参数α和n的影响较大。随着初始浓度的增加，α逐渐减小，n值逐渐增大，而α和n是与土壤进气值有关的参数［26］，根据公式（2）计算出CK、LC、MC 和HC 土壤进气值分别为16.66、42.82、157.97 和170.45，可以看出随着初始浓度的增加，土壤进气值增大，土壤中的水分较难排出，土壤含水量变化越小，这可能是导致试验初期含水量模拟精度较差的原因。

表2 不同处理水分模拟精度评价Tab.2 Evaluation of moisture simulation accuracy of different treatments

图6 不同处理土壤含水量实测值和模拟值对比Fig.6 Comparison of measured and simulated values of soil water content in different treatments

HYDRUS-1D 模型盐分模块没能很好地反映盐分在表面形成结皮之后的盐分迁移状况，因此采用分段模拟的方法，以更好模拟盐结皮土壤盐分动态变化过程，根据盐结皮在不同蒸发阶段的形成过程及作用机制［32］将试验过程分为3个阶段，0～a阶段盐结皮开始形成并覆盖土壤表面，a～c 阶段盐结皮逐渐生长变厚阶段，c～e阶段盐结皮稳定阶段。图7所示为含盐处理土壤含盐量的模拟值，模拟精度和参数见表3，其中LC 各阶段以及MC 和HC a～c 阶段、c～e 阶段模拟精度较好，R2均大于0.85，而MC 和HC 0～a 阶段模拟精度相对较差。对比含盐处理弥散系数（D），随着初始盐分浓度的增加，D逐渐增大，可能原因是不同盐分处理其溶液流速不同，从而导致D产生较大差异。对比含盐处理不同阶段，在试验a～c、c～e 阶段D小于0～a 阶段，可能是试验后期盐结皮改变了土壤孔隙度，从而对溶液D产生影响。

表3 含盐处理不同时间溶质运移参数及模拟精度Tab.3 Solute transport parameters and simulation accuracy at different times of salt treatments

图7 含盐处理土壤含盐量模拟值与实测值对比Fig.7 Comparison of simulated and measured soil salt content in saline treatments

由于试验条件限制，试验仅观测到MC 和HC 处理0～2 cm 土层的含水量，无法明确蒸发锋是否下移。因此，通过HYDRUS-1D 模型对表层0.1 cm 深度土壤含水量进行模拟，在试验结束时MC 和HC 土壤含水量分别为0.143 cm3·cm-3和0.185 cm3·cm-3，2个处理在0.1 cm 处土壤含水量远大于风干土含水量。结果表明，在试验周期内MC 和HC 土壤在0.1 cm 及以下深度仍保持着水力连接，未发现蒸发锋的明显下移。

通过模型延长蒸发时间，继续对MC 和HC 土壤含水量进行监测。结果发现，MC 在71 d 时2 cm 深度土壤含水量发生突变，此时蒸发锋位于该深度，并且随着时间的推移蒸发锋继续下移，在80 d 时蒸发锋下移至4 cm 深度，90 d 下降到6 cm 深度，在100 d 蒸发锋基本稳定在15 cm 深度。对于HC，模拟时间到100 d 时，发现蒸发锋位于土壤2 cm 深度，120 d 蒸发锋下移并稳定在8 cm 深度。上述结果表明随着盐结皮厚度的增加，表层土壤变干的速度明显减缓，蒸发锋下移的时间被极大的延长。

为明确含盐处理蒸发锋区域土壤含水量的具体变化，通过模型对含盐处理土壤含水量发生突变的土层按0.1 cm划分剖面网格对其土壤含水量进行模拟，发现LC、MC 和HC 土壤含水量突变的土层区域仅为1 cm 左右。结果表明，即使盐结皮厚度存在差异，但含盐处理蒸发锋区域均较窄（1～2 cm），该结果与Konukcu等［19］提出的较宽蒸发过渡区的观点不一致，试验与模型的结果均验证了Gran 等［3］提出的较窄蒸发锋的观点。