脱硫石膏添加下盐碱土水分入渗特征分析

高小龙, 王幼奇, 白一茹, 阮晓晗, 钟艳霞

(1.宁夏大学地理科学与规划学院,750021,银川; 2.宁夏大学生态环境学院,750021,银川;3.西北土地退化与生态恢复国家重点实验室培育基地,750021,银川; 4.西北退化生态系统恢复与重建教育部重点实验室,750021,银川)

土壤盐碱化严重制约着干旱、半干旱地区生态恢复及农业可持续发展[1]。银北平原地处中温带干旱区，降雨稀少、蒸发强烈，加之长期大水漫灌，使该区域耕地盐碱化问题严重，盐碱面积达14.8万hm2[2]。盐碱化使土壤颗粒分散、质地黏重、结构性差[3]，破坏作物生长环境[4]，导致农作物品质和产量急剧下降，严重危害区域粮食生产与安全。脱硫石膏在银川宁东化工基地产生量大，其富含植物生长所需的矿物质元素Ca、S、Si等，被广泛应用于盐碱地改良，在改善土壤物理性质[5]、提高作物产量[6]、降低盐碱土pH[7]等方面均有显著效果。

近年来国内外学者对脱硫石膏改良盐碱地做了大量研究，表明脱硫石膏可以有效置换出土壤中的Na+[8]、调节pH、提高土体中Ca2+含量[9]等，对盐碱土化学性质改良效果显著；同时，一些学者也开始关注脱硫石膏对盐碱土物理性质及水力学性质的改良，发现脱硫石膏施入盐碱土中会显著影响土壤水盐运移过程[10]、改善土体团聚体结构[11]并抑制土壤返盐过程。如程镜润等[7]发现不同脱硫石膏质量配比下，土体中离子交换过程活跃，土壤团粒结构得到明显改善。Blum等[12]通过野外实验发现石膏能改善盐碱土理化特性。张继红等[5]通过室内模拟实验发现施加石膏能改变盐碱土孔隙状况，改善盐碱土水分运动特性。综上可知脱硫石膏不仅可以降低盐碱土pH和盐离子浓度，而且可以改善土壤物理及水力学性质，促进作物增产。现阶段关于脱硫石膏添加下盐碱土水分入渗特性及数值模拟的研究较少，针对银北地区盐碱土水分运动特性的研究亦相对缺乏；而土壤入渗特性及水分分布特征直接影响作物生长发育及水资源利用效率，是制定合理灌溉措施必须参考的物理数据。

利用脱硫石膏对银北平原盐碱土进行改良，既可以实现该区域资源的合理开发和利用，又能降低盐碱化改良成本[10]。笔者通过室内一维垂直入渗试验，研究脱硫石膏添加下银北平原盐碱土水分入渗特征，比较脱硫石膏添加下不同入渗模型的适用性，并利用一维代数模型模拟土壤剖面含水率，分析模拟值与实测值之间的差距，进而明确脱硫石膏对银北平原盐碱土水分入渗过程的影响，以期为研究脱硫石膏及入渗模型更好的应用于土壤改良及生产生活等提供参考。

1 材料与方法

1.1 试验材料

供试土样取自宁夏回族自治区平罗县盐碱荒地(E 106°31′，N 38°57′)。取0～40 cm的表层土壤，去除土样中的杂物，将土壤放在阴凉、干燥的地方自然风干，过2 mm筛备用。土样类型为灰钙土，风干后土壤板结成块，土表非常坚硬，其理化性质见表1。

表1 土壤基本理化性质

脱硫石膏取自宁东煤化工基地，其主要成分是二水硫酸钙(CaSO4·2H2O)，经风干、晾晒后呈白色粉末状，粒径为30～60 μm，含水率在1%左右，纯度为90%～93%。将其碾磨过筛后放置于室内干燥处备用。

1.2 试验方法

分别将脱硫石膏与所备土样按0(CK)、1.5%(T1)、3.0%(T2)、4.5%(T3)、6.0%(T4)和7.5%(T5)的质量比混合均匀，然后将混合好的土壤样品按密度1.5 g/cm3每5 cm分层填装到内径为10 cm、高为50 cm的垂直土柱中，为消除壁面优先流的影响，装土前在土柱内壁均匀涂抹少量凡士林。每层填装完成后对土表进行打毛，以防入渗过程中出现分层现象，土柱填装高度为40 cm。土柱填装完成后在土壤表面覆一层滤纸，用相同规格的马氏瓶对其进行供水，供水过程中保持3 cm稳定水头。供水开始后，观察不同时刻湿润锋运移进程，读取前后左右4个湿润锋点，取平均数记录湿润锋运移距离，同时记录马氏瓶凹液面数值。前5 min每10 s记录1次，5～10 min每30 s记录1次，10～20 min每1 min记录1次，20～30 min每5 min记录1次，30～60 min每10 min记录1次，60 min以后每30 min记录1次。当湿润锋到达土柱底部时立刻停止供水，记录马氏瓶最终读数，入渗试验结束。迅速吸干土柱表面积水，立即挖取不同深度处(每隔5 cm取样)土样装入铝盒，用烘干法测定土壤含水率。每个处理设置3次重复。

1.3 入渗模型

笔者通过常用的Philip模型[13]、Kostiakov模型[14]、Horton模型及通用经验模型[15]对脱硫石膏添加下土壤入渗速率进行模拟，其表达式如下：

Philip模型公式：

i=0.5St-0.5+ic。

(1)

式中：i为入渗速率，cm/min；S为土壤吸渗率，cm/min0.5；t为入渗时间，min；ic为稳定入渗速率，cm/min。

Horton模型公式：

i=ic+(i1-ic)e-kt。

(2)

式中：i1为初始入渗速率，cm/min；k为入渗模型参数。

Kostiakov模型公式：

i=ct-d。

(3)

式中：c为初始入渗速率的变化，cm/min；d为水分入渗速率衰退程度，量纲为1。

通用经验模型公式：

i=v+ft-n。

(4)

式中：v的绝对值表征稳定入渗速率，cm/min；f的绝对值表征初始入渗速率，cm/min，n为模型参数，量纲为1。

1.4 一维代数模型

垂直一维非饱和土壤累积入渗量[16]计算公式为：

(5)

式中：I为累积入渗量，cm；Zf为湿润锋，cm；θs为饱和含水率，cm3/cm3；θr为滞留含水率，cm3/cm3；θi为初始含水率，cm3/cm3；α为非饱和导水率形状系数。风干后土样初始含水率较小，假定θr=θi，所以累积入渗量可表示为：

(6)

不同深度处含水率的计算公式如下：

(7)

式中：θ为土壤剖面含水率，cm3/cm3；Z为任意深度，cm。根据式(7)得出的α可以估算入渗一定深度Zf时，不同土层含水率的变化。

1.5 数据处理与分析

采用Excel 2013软件作图，Origin 2018和IBM SPSS Statistics 26软件进行模型拟合，通过单因素方差分析确定各处理入渗速率间的差异性，利用均方根误差RMSE、平均绝对误差MAE和符合度指数D[17]对模拟结果进行评价。

2 结果分析

2.1 对湿润锋的影响

湿润锋可以表示一维入渗过程中土壤基质势和重力作用下土壤的水分运动特征[18]。不同脱硫石膏添加比例下水分入渗过程中湿润锋随时间变化如图1所示，可见脱硫石膏添加比例不同，湿润锋运移进程不同。整个入渗过程中随入渗时间推移各处理湿润锋均呈先快后慢的增加趋势。入渗初期120 min 内各处理湿润锋差别较小，随入渗时间增加，脱硫石膏配比对土壤水分入渗过程的影响逐渐显著。入渗400 min后，同一时刻，脱硫石膏添加比例越大湿润锋运移距离越小。入渗结束时各处理累积用时较CK分别增加9.09%、13.64%、22.73%、27.27%和36.36%。总体来看，入渗各阶段脱硫石膏添加对湿润锋运移过程均有抑制作用。

为更好探究脱硫石膏添加量对湿润锋运移过程的影响，根据湿润锋趋势，通过幂函数对湿润锋变化进行拟合。如图1所示，幂函数(Zf=atb)较好地反映了各处理湿润锋Zf与入渗时间t之间的关系，R2均>0.99，很好拟合不同脱硫石膏添加比例下湿润锋随时间的变化情况。可见脱硫石膏添加减缓水分入渗过程，随脱硫石膏施量增加减缓程度亦增加。

脱硫石膏添加比例：CK：0；T1：1.5%；T2：3.0%；T3：4.5%；T4：6.0%；T5：7.5%，Zf为湿润锋，a，b均为模型参数量纲为1。下同。Desulfurization gypsum addition ratios: CK: 0; T1: 1.5%; T2: 3.0%; T3: 4.5%; T4: 6.0%; and T5: 7.5%, the same below. Zf: Wetting front, cm; a, b are model parameters. 图1 脱硫石膏添加对湿润锋的影响Fig.1 Effects of adding desulfurization gypsum on wetting front

2.2 对累积入渗量的影响

由图2可见，累积入渗量随脱硫石膏添加比例变化趋势与湿润锋相似。入渗初期120 min内不同处理间土壤累积入渗量差异不大，各处理累积入渗量均迅速增加。随入渗历时推进，脱硫石膏对土壤水分入渗过程的影响逐渐增大，480 min时各处理累积入渗量存在明显差异，CK、T1、T2、T3、T4和T5的累积入渗量分别为8.22、7.66、7.43、7.19、6.90和6.55 cm，与CK相比T1、T2、T3、T4和T5的累积入渗量分别减少6.91%、9.64%、12.53%、16.08%和20.30%，可见同一时刻，累积入渗量随脱硫石膏增加逐渐减小。历时960 min时不同处理的累积入渗量分别为11.58、10.83、10.46、10.12、9.80和9.37 cm，T1、T2、T3、T4、T5较CK分别减小6.46%、9.68%、12.64%、15.37%和23.65%。累积入渗量一定程度上可以表征土壤水分入渗能力。相同入渗时间，脱硫石膏添加比例越高，累积入渗量越小，说明脱硫石膏添加减弱土壤水分入渗能力。脱硫石膏在与盐碱土相互作用后降低土壤pH值的同时[7]，土壤孔隙度发生改变，入渗能力降低，这极大减缓灌区土壤水分的快速下渗，阻碍地下水与灌溉水的连通过程，有利于防止灌区土壤次生盐碱化危害。

图2 脱硫石膏添加对累积入渗量的影响Fig.2 Effects of adding desulfurization gypsum on cumulative infiltration

2.3 对入渗速率的影响

入渗速率是表征土壤水分入渗快慢的直接指标，根据入渗不同阶段土体主要作用力的不同，一般分为3个阶段：瞬变阶段(分子力作用下入渗速率较高的初始入渗阶段)、渗漏阶段(重力及毛管力作用下入渗速率迅速下降的中期入渗阶段)和稳渗阶段(重力作用下入渗速率趋于稳定的终期入渗阶段)[19]。表2和图3表征不同脱硫石膏添加下土壤水分入渗速率的变化情况，可见随脱硫石膏施用量增加各阶段入渗速率均呈下降趋势，但不同时间段土壤水分入渗速率下降程度不同。为更清楚反映不同处理间的差异性，将图3局部放大得到20～720 min 时的入渗速率。在初始入渗阶段(即瞬变阶段)各处理入渗速率均较大，初始入渗速率差异不显著(P>0.05)。在渗漏阶段各处理入渗速率逐渐趋于稳定，呈迅速下降趋势，不同处理表现出明显的规律性，即脱硫石膏添加比例越大，土壤水分入渗速率越小，与CK相比90 min时T1、T2、T3、T4和T5入渗速率分别减小16.7%、33.3%、44.4%、47.2%和58.3%。在稳定入渗阶段，各处理入渗速率均较为平缓，呈稳中变缓的态势，稳定入渗速率在CK与T2、T3、T4、T5之间存在显著性差异(P<0.05)，T2与T4、T5之间存在显著性差异(P<0.05)。该入渗过程中除T3与T4外，其余各处理平均入渗速率均表现出显著差异性(P<0.05)。

表2 不同脱硫石膏添加比例对入渗速率的影响

图3 脱硫石膏添加对入渗速率的影响Fig.3 Effects of adding desulfurized gypsum on infiltration rate

2.4 土壤水分入渗过程模型拟合

为确定脱硫石膏添加下土壤水分入渗模型适用性及相关水力学参数，对Philip、Horton、Kostiakov和通用经验模型进行拟合分析，拟合参数如表3所示。Philip模型中决定系数(R2)均>0.840，S在0.617～0.557之间。S反映土壤水分入渗能力，其值越大土壤水分入渗能力越强，可见Philip模型较好地拟合土壤水分入渗过程，S随脱硫石膏添加比例增加逐渐减小，说明添加脱硫石膏对盐碱土水分入渗有抑制作用，且随脱硫石膏施用量增加其抑制作用增大。Horton模型中R2为0.748～0.753，随脱硫石膏添加比例增加，i1与ic均呈减小趋势，与实测数据变化规律吻合。Kostiakov模型中R2均>0.873，c值为0.312～0.281，随脱硫石膏配比增加逐渐减小，说明一维垂直入渗过程中脱硫石膏添加比例越大土壤初始入渗速率越小。d值为0.394～0.415，其大小反映水分入渗速率的衰退程度，可见随脱硫石膏添加比例增加其值越大，说明脱硫石膏可以增强土壤入渗速率的衰退程度，减缓水分入渗过程。在通用经验模型中R2均>0.911，v和f的绝对值分别表征稳定入渗速率和初始入渗速率，可见随脱硫石膏配比增加v、f的绝对值均呈减小趋势，即脱硫石膏添加量增加土壤入渗过程中稳定入渗速率和初始入渗速率降低。

表3 不同脱硫石膏添加比例下入渗速率模型拟合

2.5 土壤剖面含水率分布模拟

入渗结束时各处理不同土层含水率变化情况如图4所示，随土层深度增加不同处理土壤含水率均逐渐减小。入渗结束时(湿润锋达到40 cm深度处)，35～40 cm深度处CK、T1、T2、T3、T4和T5土壤含水率分别为26.84、26.25、25.84、25.50、24.70和24.11 cm3/cm3。与CK相比，T1～T5处理在 35～40 cm 土层深度处土壤含水率分别减少0.59%、1.00%、1.34%、2.14%、2.73%。通过式(6),利用累积入渗量曲线和湿润锋曲线获得二者间斜率，结合饱和含水率、初始含水率数据获得非饱和导水率形状系数α(表4)，将α代入式(7)求得每个处理不同土层深度的模拟含水率数值(图4)，通过土柱结束时所测不同土层含水率数值，与模拟值进行比较分析。用RMSE、MAE和D来反映预测值与实测值间的接近程度(表4)。不同处理下土壤含水率MAE在1.064%～1.828%之间、RMSE在1.299%～1.987%之间，其值均<2%，说明模拟值与实测值差距很小，拟合效果较好。D为标准化度量指标，范围在0～1之间，其值越接近于1，说明模型拟合效果越好、实测值与模拟值越接近[17]。本研究中各处理D值分别为：0.883、0.821、0.770、0.719、0.904、0.600，除T5处理外，其他各处理D值均>0.7，T4>0.9，说明一维代数模型可以较好模拟脱硫石膏添加下盐碱土不同深度处含水率的变化情况。

图4 不同处理土壤剖面含水率Fig.4 Soil profile moisture contents under different treatments

表4 一维代数模型拟合及评价

3 讨论

入渗前期脱硫石膏添加对湿润锋、入渗速率和累积入渗量的影响较小，随时间推移影响越来越显著，这可能是由于入渗前期土体干燥，土壤水吸力均较大[20]、脱硫石膏与盐碱土复合作用效果不显著，各处理间入渗过程差异不明显。随入渗过程推进，脱硫石膏与盐碱土充分反应，土壤中的Na+被石膏中的Ca2+置换出来，Ca2+与土壤胶体表面负电荷的吸引力及对土壤胶体的聚沉作用比Na+强[21]，同时脱硫石膏会填充土壤孔隙且其吸水性强[11]，所以土壤水分运移过程变缓，这与张继红等[5]研究得出的结论相一致。说明脱硫石膏添加能够减缓银北平原盐碱土水分入渗过程，减少土壤养分流失，防止农药、化肥等深层渗漏，阻碍地表水与地下水的连通作用。

Kostiakov模型属典型的指数模型，当时间轴趋于0时i值趋于无限大，与实际情况不符[19]，不适合描述入渗初期水分运动特性。且Kostiakov模型属于半经验模型，其参数没有明确的物理意义，相关学者在进行模拟入渗实验时也没有选用该模型[19,22]。Horton模型拟合的R2在0.748～0.753之间，拟合效果低于其他3种模型，上官玉铎等[23]认为Horton模型参数较多、拟合难度较大。通用经验模型拟合效果较好，但v值为负[24]。Philip模型不仅可以较好的模拟土壤水分入渗过程，而且能够求算该过程中土壤吸渗率的大小，反映脱硫石膏添加对盐碱土入渗能力的影响。综上，根据4种入渗模型参数的实际物理意义及模型拟合效果[25]，发现Philip模型更适合拟合脱硫石膏添加对银北地区盐碱土水分入渗过程的影响。

各处理含水率随土层深度增加逐渐减小，同一土层不同处理土壤含水率随脱硫石膏添加比例增大而减小，这与梁嘉平等[11]得出的结论一致。说明脱硫石膏添加改变各处理不同土层剖面含水率，这主要是因为脱硫石膏添加使土壤中钙质胶体增加，水的黏滞性变大，毛管对水的吸渗作用减弱[5]，且石膏颗粒对土壤孔隙有填充作用，影响土壤水分分布状况。准确拟合土壤含水率对农业灌溉制度、土壤水分调控技术的确定及促进作物增产等具有重要意义。本研究通过一维代数模型预测各处理任意深度处的土壤剖面含水率，模拟结果如图4所示，实测值与模拟值之间差异较小，可以用其模拟脱硫石膏添加下盐碱土的水分入渗情况。王春霞等[26]的研究也表明一维代数模型在土壤水分一维垂直入渗过程中有很好的适用性。针对重度盐碱土，赵连东等[17]也利用一维代数模型高精度模拟出不同深度处土壤含水率。综上，一维代数模型在预测盐碱土水分分布情况时有较好的适用性。

4 结论

1)对比不同处理下脱硫石膏添加对银北平原盐碱土水分入渗特性的影响，发现添加脱硫石膏对盐碱土水分运动过程有明显减渗作用。入渗结束时T1、T2、T3、T4、T5处理累积用时较CK分别增加9.09%、13.64%、22.73%、27.27%、36.36%。

2)综合4种入渗模型的实际物理意义及拟合效果，发现Philip模型更适合描述脱硫石膏添加下银北地区盐碱土水分入渗过程。

3)添加脱硫石膏改变土壤剖面含水率，同一土层深度处，脱硫石膏添加比例越大土壤剖面含水率越低。一维代数模型能很好预测脱硫石膏添加下盐碱土水分入渗过程中土壤剖面含水率的分布情况。