不同沉陷应力区土壤水分和溶质运移的模拟试验

2023-10-08尚海丽黄显武李建伟李依临杨宏宇

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

温欣，尚海丽，黄显武，李建伟，李依临，杨宏宇

（1.中国矿业大学（北京），煤炭资源与安全开采国家重点实验室，北京 100083；2.内蒙古科技大学矿业与煤炭学院，内蒙古包头 014010）

煤炭资源长期占据我国能源主体地位。中国山西-陕西-内蒙古交界地长期高强度煤炭开采，造成地表沉陷，发育不同规模的裂缝，土壤内部发育微裂隙和陷落结构。由于土壤结构破坏，导致地表蒸发量增大，地表水也会沿裂缝向深部渗漏，造成地表水资源流失。同时由于“盐随水走”，地表沉陷不仅造成土壤水分空间变异，而且影响土壤盐分运移，加大了土壤盐碱化风险。开采沉陷对土壤水盐运移的改变会对土壤微生物群落代谢活动、地表植被和土壤养分循环等地表生态环境产生巨大影响，是山西-陕西-内蒙古交界地矿区地表生态退化的重要因素［1-3］。因此，研究开采沉陷地土壤水盐运移规律是矿区生态修复的重要组成部分，对科学评价开采沉陷的环境效应具有重要意义。

近年来，国内外学者［4～7］利用田间试验、物理模型、数值模拟等方法研究采煤沉陷地土壤水盐运移规律，主要集中于沉陷裂缝对土壤水盐和植被根系生长的影响、沉陷盆地土壤水盐时空变异特征等方面。其中，Tripathi 等［4］以印度热带干旱地区煤田为研究对象，通过对沉陷斜坡、沉陷洼地、相邻未受损地区的对比研究发现，沉陷导致土壤容重、土壤持水能力在沉陷洼地增加，而在沉陷斜坡降低。Xiao等［5］以山东省济宁市东滩煤矿为研究对象，对土壤含水量的空间变异性和相关性进行分析评价，表明土壤含水量和容重具有强烈的空间依赖性。张健［6］通过室内物理试验与数值模拟结合的方式模拟不同尺度地裂缝周围土壤水分空间分布变化，研究结果表明地裂缝影响周围土壤水分空间分布，加速裂缝周围土壤水分散失，且地裂缝尺度越大，其影响范围越大；越靠近裂缝壁，土壤含水量越低；随土壤深度增加，裂缝的影响作用逐渐降低。李雯雯［7］选择典型黄土矿区大佛寺煤矿为研究区，采用遥感影像反演、实地监测、光谱数据建模等方法，从不同尺度分析采煤沉陷区土壤湿度的时空演变特征，研究结果表明采煤沉陷引起地表拉伸变形和裂缝发育，造成土壤密实度降低、土壤水分入渗和蒸发作用增强。

随着土壤水盐运移理论模型的深入研究以及大量田间试验数据的支撑，前人开展了大量水盐运动过程定量化研究。HYDRUS 水盐运移模拟模型是目前应用最广泛的定量水盐运移方法，因其边界条件灵活，主要适用于室内土柱模拟［8-9］。其中，Liu等［10］为了提高暗管排水中HYDRUS-2D模拟的准确性，提出阻力层、实际入渗面和等效半径处理3种方法，来代替HYDRUS-2D地下管道的标准渗流边界，大大提高了土壤含盐量、含水率的模拟精度。王世明等［11］通过HYDRUS 模拟防护林建设初期不同矿化度地下水灌溉条件下土壤水盐分布特征和时间变化，结果表明表层（0～30 cm）土壤含水量和含盐量受灌溉影响大，数值波动剧烈，且土壤盐分表聚强烈；深层（50～150 cm）则受灌溉的影响小，数值波动小。于洋等［12］通过GMS和HYDRUS的结合，构建包气带-饱和带的水盐耦合模拟模型，以半干旱地区宁夏中卫某速生林高阶地为例，采用历史数据进行了验证，结果表明该模型可有效模拟包气带-饱和带的水流运动和氯离子的运移情况。邹慧［13］结合室内试验HYDRUS-2D水盐模型对神东矿区煤炭开采中土壤水分动态、土壤水分入渗、蒸发及其对土壤物理性质的影响进行研究，结果表明受地形、植被和开采的影响，表层土壤水分呈现差异显著的斑块状分布格局，沉陷和裂缝引起了土壤含水量等土壤物理性质的变化。HYDRUS-2D可以模拟土壤水盐在水平方向和任意复杂地形条件下的迁移问题。因此，对于有效预测点源、面源程度上对陕西-山西-内蒙古交界地的采煤沉陷地土壤水盐运移研究具有重要意义。

在前人研究的基础上，为了揭示开采沉陷微裂隙对土壤水分和溶质运移的影响，本文采用不同沉陷应力区土壤容重的变化来描述沉陷微裂隙的变化，利用自主研发的土壤沉陷模型进行模拟试验，研究蒸发过程中土壤水分、总盐分、不同溶质离子在不同沉陷应力区随土壤深度和时间变化的规律；并通过采集沉陷剖面不同应力区土壤容重数据，建立基于HYDRUS-2D的不同沉陷应力区土壤水盐运移数值模型，为系统揭示开采沉陷对土壤水盐运移的影响机制、科学评价开采沉陷的环境效应提供理论依据。

1 材料与方法

1.1 沉陷物理模型建立

本研究利用课题组自行设计的沉陷模拟装置进行沉陷模拟试验（图1a）。该装置为一敞口长方体有机玻璃，尺寸：150 cm×30 cm×120 cm。A 底面由15个活动面板组成。每个活动面板尺寸：10 cm×30 cm×1 cm。每个活动面板由B底面4个活动螺钉支撑，可通过扭动螺钉来调节活动面板升降高度，最大可调高度15 cm。在A 底面活动面板上方分层填充土壤，通过调节螺丝高度模拟沉陷过程，A底面两端各留3 个活动面板，高度不变（高度为15 cm）。从左边第4个活动面板向右依次调节9个活动面板高度，模拟沿剖面走向沉陷盆地依次向右推进，上覆土层从左至右依次沉陷，形成沉陷剖面。土壤剖面中不同沉陷应力区的划分是基于开采沉陷盆地的理论模型（图1b）而来。沉陷盆地中心部位分布挤压应力，定义为挤压区。该挤压区两侧分布拉张应力，开采裂缝集中发育于该部位，定义为拉张区。靠近模型边界部位受到物理模型边界效应影响较大，因此该拉张区域数据未采用，定义为边缘区。

图1 不同沉陷应力区特征模拟装置及开采沉陷盆地理论示意图Fig.1 Simulation device of different subsidence stress regions characteristics and theoretical diagram of mining subsidence basin

1.2 试验过程和数据采集

试验在内蒙古科技大学微生物复垦试验室进行。试验设2种处理，试验组为沉陷组，对照组为无沉陷处理，每种处理设3 个重复。试验所用土壤是采自内蒙古包头地区的砂土，过2 mm 筛去除杂物后，风干、备用。根据粒度分析数据，该砂土与山西-陕西-内蒙古交界地采煤沉陷地的风积砂粒度组成相近（表1）。

表1 供试土壤基本理化性质Tab.1 Basic physical and chemical properties of the tested soil

模型搭建：利用图1所述沉陷模型，为避免无支撑条件下降低活动面板造成土层垮塌的问题，首先在A底面活动面板上铺设0.95 cm厚度的石膏板，为上覆土体提供近似于覆岩的支撑力。然后将供试土壤与矿化度3 g·L-1（离子浓度比：Na+:Cl-:CO32-:Ca2+:SO42-:Mg2+=8:4:1:1:1:1）的水溶液充分混合，维持最大饱和含水率为21%。分层装土，每层厚10 cm，控制试验土壤容重为1.6 g·cm-3，层面打毛，防止分层效应。同时安装25 个监测点的土壤水分采集器。静置24 h，待水分平衡后，对试验组进行沉陷模拟。

沉陷模拟过程：从左边第4 个活动面板开始向右依次下降9个活动面板，为达到最佳沉陷效果，每隔1 h 从左至右水平推进10 cm，即下降1 个活动面板。24 h后达到稳沉阶段（图2）。

图2 稳沉阶段模拟土柱实物图Fig.2 Physical diagram of simulated soil column in steady settlement stage

蒸发过程实时监测：利用预先安装的EM50 数据采集器+5TE 传感器进行土壤水分和电导率数据实时监测。EM50 监测仪（壤博士RS485）测量间隔时间为30 min，仪器精度：电导率0～10000 μS·m-1±3%FS·m-1、含水率0～50%±2%。受5TE传感器监测范围限制，土壤含水率达到12%，试验结束。蒸发过程为自然蒸发，室内温度20～25 ℃，湿度13%～18%。

考虑物理模型的边界影响，本试验2 个边缘区数据未采用。

由于传感器所测电导率受土壤含水率影响，因此采用土水比1:5的土壤浸提液分别测定25个观测点的土壤浸提液电导率（EC1:5，dS·m-1）。土壤自动传感器所测电导率（EC5TE，dS·m-1）与含水率（θ，%）的转换关系为［14］：

土壤EC5TE转换为EC1:5的换算关系为：

土壤EC1:5转换为土壤含盐量（SSC，g·kg-1）的计算公式为［15］：

1.3 HYDRUS-2D水盐模型建立

1.3.1 土壤水分运动基本方程在采煤沉陷条件下，土壤水分运动简化成径向和垂直2个方向，属于二维运动水分入渗问题。假设试验土壤均质，各向同性。在考虑非饱和土壤水滞后效应情况下，土壤水分运动规律的数学模型［16-17］为：

式中：θ为土壤体积含水率（cm3·cm-3）；φ为土壤水负压（cm）；K(θ)为非饱和土壤导水率（cm·d-1）；z为垂向坐标，向上为正；x为横向坐标；t为时间。

1.3.2 土壤溶质运移模型在均匀介质中，用可控的对流弥散方程模拟非反应离子运移，其数学模型［18］为：

式中：i、j为x、z轴坐标；c为液体中的溶质浓度（g·cm-3）；q为土体中水的流速；Dij为水动力弥散系数；t为时间。

1.3.3 土壤水分特征方程土壤的水力特征方程一般采用VG模型，其表达式［16］为：

式中：Ks为土壤饱和导水率（cm·d-1）；θe为土壤相对饱和度（cm3·cm-3）；θr为土壤相对饱和导水率（cm3·cm-3）；θs为土壤饱和体积含水率（cm3·cm-3）；h为负压水头（cm）：α、n为试验测定经验函数；m=1-1/n；l为形状系数，通常取平均值0.5。

1.4 模型定解条件

土壤水分运移的边界条件如图3所示。

图3 不同沉陷应力区土壤水分和溶质运移数学模型Fig.3 Mathematical models of soil moisture and solute transport in different subsidence stress regions

土壤溶质运移的初始条件：

土壤盐分运动上边界为大气边界：

下边界为浓度边界：

式中：C0为土壤初始含盐量（g·kg-1）；z为垂向坐标，向上为正；qs为地表水分通量；Cs为上边界电导率浓度（g·cm-3）；Cb为下边界电导率浓度（g·cm-3）；t为时间。

1.5 模型参数

1.5.1 土壤水力特征参数采用土壤容重指标量化地表变形程度。利用供试土壤粒径、沉陷剖面不同应力区25 处观测点实测容重数据，由HYDRUS-2D的Rosetta程序自动生成土壤水力特征参数（表2）。

表2 土壤水力特征参数Tab.2 Soil hydraulic characteristic parameters

1.5.2 单元划分及离散根据试验过程中土柱沉降变化与实测土壤容重数据建立二维模型，模拟沉陷土壤剖面。试验组与对照组模拟区域与尺寸如图3所示。采用不等间距网格对模拟区域进行离散化。试验组生成节点1147 个，离散网格2500 个。对照组生成节点196个，离散网格4802个。时间离散单位为d，模拟时段为2020 年9 月18 日—2021 年11月4日，共计78 d。设定初始时间步长0.1 d，最小时间步长0.001 d。

2 结果与分析

2.1 HYDRUS-2D水盐模型精度检验

为了确保HYDRUS-2D 水盐模型的可靠性，将土壤含盐量、含水率实测值与模型的模拟值进行平均相对误差（ME）、均方根误差（RMSE）和决定系数（R2）检验（表3）。当ME≤0.5、RMSE≤0.5、R2＞0.95时，模拟精度满足试验要求。就整体而言，土壤水盐模拟值与实测值吻合度较好，模拟结果能较好地反映土壤水盐的动态变化。这证明HYDRUS-2D水盐模型可以较好地模拟不同沉陷应力区土壤剖面中的水盐运移。

表3 模型精度检验Tab.3 Tests of model accuracy

2.2 不同沉陷应力区土壤含水率时空变化特征

沉陷拉张作用增强了土壤蒸发作用，拉张区土壤含水率显著低于挤压区，土壤含水率随深度变化显著（图4）。＜40 cm深度土壤含水率随时间变化曲线呈现典型蒸发过程“三阶段”：稳定蒸发阶段（0～5 d），蒸发速率下降阶段（5～40 d）和蒸发水分扩散阶段（40 d以后）。在蒸发全过程中，0～20 cm沉陷挤压C区土壤含水率显著高于两侧沉陷拉张区，与对照区差异较小。这说明沉陷剖面2个拉张区土壤微裂隙发育，增大了土壤与大气的接触面，蒸发作用强烈。但是，20～40 cm、40～60 cm左侧拉张B区和挤压C区的土壤含水率无显著差异。类比随开采推进方向沉陷盆地的发育过程，这是由于随着沉陷向右推进，左侧20～60 cm 沉陷拉张裂缝闭合的原因。＞40 cm深度土壤含水率在蒸发作用初期显著增大，可能是上层土壤水分向下运移造成的下层土壤含水率波动。＞60 cm 深度挤压C 区土壤含水率仍然显著大于两侧拉张区。这说明沉陷剖面深部拉张区对土壤含水率的扰动时效性比浅层更久。

图4 不同深度不同沉陷应力区土壤含水率随时间变化Fig.4 Variation of soil moisture with time in different subsidence stress regions at different depths

蒸发试验结束时，不同沉陷应力区土壤含水率随深度的分布剖面图（图5）表明，0～40 cm 深度，沉陷拉张作用显著增强土壤蒸发作用，造成挤压C 区土壤含水率显著大于拉张区。40～60 cm 深度不同沉陷应力区土壤含水率无显著差异。＞60 cm 深度土壤含水率又呈现随沉陷应力区不同而分异的特征。

图5 第78 d时土壤水分剖面模拟Fig.5 Simulation of soil moisture profile on the 78th day

2.3 不同沉陷应力区土壤总含盐量时空变化特征

不同深度土壤总含盐量随时间逐渐减小，以0～20 cm 深度土壤随时间脱盐量最大（图6）。这是表层土壤强烈蒸发作用导致盐分在0～5 cm 深度形成盐壳的结果。0～60 cm 深度左侧沉陷拉张B 区土壤脱盐量随时间逐渐增大，与挤压C 区差异显著。＞60 cm 深度土壤总含盐量随时间变化曲线无显著差异。

图6 不同深度不同沉陷应力区土壤总含盐量变化Fig.6 Changes of soil total salt content in different subsidence stress regions at different depths

沉陷不同应力区土壤总含盐量随深度呈强烈变异（图7）。蒸发试验10 d，0～20 cm土壤盐分受蒸发作用影响，向土壤最表层迁移，沉陷土壤总含盐量高于对照组，右侧拉张D 区盐分积累量显著高于其他应力区，盐分最终在土壤表层0～5 cm 形成盐壳。在50 d、78 d，对比对照组，沉陷组各应力区土壤总含盐量不但随深度呈现强烈变异，具体表现为20～40 cm、60～80 cm 土壤总含盐量相对积聚；而且右侧拉张区总含盐量积聚深度有向下迁移的趋势。这表明沉陷区深部土壤水分不均匀分布会导致土壤盐分随深度的变异，拉张区土壤盐分向下迁移，减缓了沉陷地盐碱化风险。

图7 不同沉陷应力区土壤总含盐量模拟剖面Fig.7 Simulation profile of soil total salt content in different subsidence stress regions

2.4 不同沉陷应力区6种土壤溶质离子分布特征

沉陷拉张区和挤压区显著改变了土壤溶质离子随深度的变化，不同的溶质离子随深度呈现不同的积聚效应（图8）。在对照组中，Ca2+、SO42-、CO32-土壤离子浓度均随深度呈现典型双峰式积聚效应，积聚深度为20～40 cm、60～80 cm。对比对照组，沉陷拉张区Ca2+、SO42-、CO32-土壤离子浓度随深度呈单峰积聚，且Ca2+、SO42-离子在沉陷拉张区积聚深度下移至40～60 cm，挤压C区无显著积聚效应；CO32-离子在沉陷挤压C 区和拉张区具有相似分布规律，仅在60～80 cm积聚，呈现显著向下迁移的趋势。在对照组中，Mg2+、Na+、Cl-土壤离子浓度以单峰式积聚。对比对照组，在沉陷拉张区3种离子积聚深度均上移；在沉陷挤压区Mg2+、Cl-离子随深度无显著积聚效应；Na+离子在沉陷挤压区呈双峰式积聚。沉陷拉张区Ca2+、SO42-、Mg2+、Cl-离子积聚浓度大于挤压区。

图8 不同沉陷应力影响下6种溶质离子分布剖面Fig.8 Distribution profiles of six solute ions under different subsidence stress

3 讨论

3.1 沉陷应力区分布对土壤水分运移的影响机制

在中国山西-陕西-内蒙古三省交界地，高强度的井工开采造成风积沙地、黄土覆盖区大面积沉陷［19］。在开采沉陷扰动下和稳定沉陷盆地形成过程中，土体内部形成2类裂缝，一类为边缘裂缝，在开采沉陷过程中不断加深，且不会在稳沉后消失。另一类为动态裂缝，随开采沉陷推进，呈现拉张-闭合-拉张的动态变化过程［20］。而本文室内模拟试验中，通过依次调节活动面板模拟土壤受开采沉陷应力作用形成的拉张区和挤压区，最终形成稳定沉陷盆地。土壤容重指标量化沉陷裂缝分布和发育情况，不同沉陷应力区导致土壤容重发生显著变化，试验研究发现开采沉陷对土壤水盐运移的主要影响因素是土体裂缝的发育情况，包括宏观裂缝和土体内部的微裂隙。这些微裂隙为土壤水分运移创造了优先路径，加速土壤水盐运移过程。这也侧面印证了Hu等［5］的结论，土壤含水量和容重具有强烈的空间依赖性。这成为干旱少雨、光照强度大的西部地区地表生态损伤的重要因素。

毕银丽等［21］发现在沙地沉陷区，开采沉陷增强了60～80 cm土壤水分的变异性，沉陷区土壤蒸发量均明显大于未沉陷区。本文室内模拟试验结果表明，0～40 cm 深度，沉陷拉张作用显著增强土壤蒸发，造成挤压区土壤含水率显著大于拉张区。这说明沉陷拉张应力导致土壤内部微裂隙发育，增加土壤与大气的接触面积，蒸发作用增强，加速土壤水分运移。本试验中，＞60 cm深度土壤含水率又呈现随沉陷应力区不同而分异的特征。综合分析，浅煤层开采沉陷会在60～80 cm深度形成不同程度垮塌，垮塌在沉陷土壤剖面的水平位置具有随机性，这造成土壤水分在垮塌区变异性增强。

本试验研究结果还表明，开采沉陷对土壤水分影响具有时效性，土壤水分在一定时间后趋于稳定状态。这一规律表明开采沉陷对地表土壤、植物根系乃至地表生态系统的影响均随时间变化。随着动态裂缝的闭合，开采沉陷区地表生态系统表现出显著的自修复特点。Druzbicka等［22］通过对干旱、半干旱地区2 个废弃矿山的研究发现，采矿引起岩层暴露与沉陷促进了土地盐碱化，但通过长时间自然演替过程，受干扰土地自然恢复，地区生物多样性增强。因此，在开采沉陷地生态修复过程中，我们必须充分认识沉陷地土壤水分空间变异性与时效性，有机结合人工修复技术，加快沉陷区生态自修复进程。

3.2 开采沉陷促进干旱地区土壤盐分向下运移

“盐随水走”，土壤盐分的运移均随水分的运动而迁移［23］。因此，拉张应力区土壤微裂隙的发育不仅增强土壤蒸发作用，而且势必影响土壤盐分运移，这是导致山西-陕西-内蒙古三省交界地土壤盐碱化的重要因素。

Ashish等［24］发现底部即沉陷挤压应力区中心位置0～60 cm 深度土壤电导率显著增加，也证明沉陷应力增加了土壤盐分的空间变异性。本文试验研究结果表明，沉陷应力作用对土壤含盐量影响深度远大于土壤含水率。0～60 cm左侧沉陷拉张区土壤脱盐量随时间逐渐增大，与挤压区差异显著。这说明开采沉陷对土壤盐分空间异质性的影响比土壤水分更显著。史文娟等［25］通过室内土柱模拟试验，研究蒸发条件下不同层位夹砂层土壤剖面盐离子的动态变化情况，结果发现盐分离子向上迁移的速度随时间的增加而减小。此外，本试验结果还表明，沉陷应力使0～5 cm土壤积盐明显。在前人［25］研究的基础上，本文研究结果证明了开采沉陷作用下土壤表层盐壳的形成进一步阻止水分和盐分向上迁移。

沉陷应力作用影响土壤质量，导致沉陷区域内土壤物理和化学性质受影响，尤其在中国山西-陕西-内蒙古三省交界地易造成土地盐碱化［26］。土壤性质受到自然因素和人为因素影响，具有极强的空间变异性。本试验研究发现，沉陷区深部土壤水分不均匀分布会导致土壤盐分随深度的变异，拉张区土壤盐分向下迁移。这与毕银丽等［27］利用箱状试验模拟开采沉陷裂缝，进行HYDRUS-2D 数值模拟研究结果相似，研究结果表明裂缝区含盐量向下迁移，证明开采沉陷裂缝抑制了次生盐碱化发生。本研究发现沉陷拉张区Ca2+、SO42-、Mg2+、Cl-离子积聚浓度大于挤压区，Ca2+、SO42-、CO32-土壤离子浓度随深度呈单峰积聚，沉陷拉张区积聚深度均显著下移。这进一步证明不同盐分离子具有向下运移的趋势，且离子类型影响盐分的下移规律。

土壤养分的变异是土壤水盐变化的后效作用［28］。开采沉陷应力作用影响并改变土壤水盐运移，最终影响土壤养分，抑制植物生长。Zhao等［29］对黄土地区采煤沉陷区采样分析土壤有机质与全氮，结果表明沉陷降低了0～20 cm深度上土壤有机质与全氮含量。地面沉降增大了有机质和全氮的空间异质性。裂缝宽度、裂缝深度和取样点与裂缝边缘的距离是导致土壤有机质和全氮水平空间异质性的主要影响因素。

因此，开采沉陷对地表生态系统的影响程度很大。针对开采沉陷不同应力区的生态恢复，要开展“分区”生态修复模式，以土壤自修复为基础，建立差别化生态修复模式。