地膜覆盖具有增温保墒、增产控盐、防治杂草等优点, 被广泛应用于农业生产^[1]。我国现有农膜生产企业近千家, 总产能超过500万t, 年农膜总产量近280万t, 覆盖面积达17 657.0万hm², 并以12%的速度逐年增长^[2-3]。我国农田覆膜后平均残膜量为50~ 260 kg·hm^-2, 且主要分布在土壤耕层(0~20 cm)^[4-5]。目前残膜机械回收率低, 在土壤中难以降解^[6], 随使用年限增加土壤中残膜量越来越高^{[4, 7]}。农用残膜已成为农业环境污染的一个重点问题。2015年《农业部关于打好农业面源污染防治攻坚战实施意见》中提出要着力解决农田残膜污染^[8]。

残膜会破坏土壤团聚体结构、阻断孔隙连续性, 导致土壤容重增大、孔隙度减小及入渗阻力增加, 严重影响水和溶质在土壤中的运移过程, 阻碍农作物根系发育, 最终造成农作物品质和产量下降^[9-10]。因此近些年国内外专家学者关于残膜对土壤水分入渗影响方面做了大量研究。李仙岳等^[11]研究认为滴灌条件下土壤中残膜含量增加, 使湿润锋进程明显变缓且湿润体呈缩小趋势。解红娥等^[12]发现土壤入渗率与残膜含量呈对数关系, 当残膜量达到360 kg·hm^-2时, 入渗率显著降低。Jiang等^[13]发现残膜显著改变了土壤物理性质, 使得玉米根区水分运动受阻。王志超等^[14]研究残膜量对沙壤土和砂土水分入渗过程的影响发现残膜量越高、湿润锋运移距离越短, 其入渗速率也逐渐减缓, 残膜量对沙壤土入渗的影响更为明显, Philip入渗模型能准确描述土壤水分入渗过程。李元桥等^[15]发现残膜的存在对湿润锋横向距离影响显著, 残膜区土壤湿润体变小, Kostiakov入渗模型能较好地拟合残膜土壤入渗过程。通过已有研究可以看出, 残膜对不同土壤类型入渗过程有明显影响, 且对入渗模型的选择和拟合结果有一定差异^[16]。

目前关于银川平原不同残膜含量对土壤水分入渗影响的研究较少, 特别是入渗模型适用性评价方面的研究尤为缺乏。银川平原地处干旱、半干旱地区, 气候高寒冷凉, 年均蒸发量是降水量10余倍, 低温、干旱及盐渍化是影响该区域农业生产的主要问题^[17]。地膜覆盖具有增温、保墒和控盐等作用, 是该区域主要的农业耕作模式, 现已应用于40多种农作物的种植栽培, 并呈现持续增长的趋势, 但目前缺乏有效回收措施, 土壤中残膜量逐年增加, 造成耕地质量下降和农业面源污染等问题^[18-19]。因此, 本文选取银川平原淡灰钙土, 通过室内土柱模拟试验, 研究不同残膜含量条件下土壤水分入渗过程, 分析评价模型适用性, 揭示残膜含量对土壤水分入渗特征和湿润锋运移的影响, 为残膜污染区水分入渗模型的选择提供数据支持。

1 材料与方法 1.1 试验材料与装置

供试土壤取自银川市滩南庄农田(106°08′ E, 38°31′N), 土壤类型为淡灰钙土, 将供试土壤过2 mm筛。土壤黏粒含量4.3%, 粉粒含量7.9%, 沙粒含量87.8%, 全氮含量0.43 g·kg^-1, 全磷含量0.51 g·kg^-1, 全钾含量20.11 g·kg^-1, 全盐量1.92 g·kg^-1, pH 8.6, 容重1.45 g·cm^-3。试验残膜选用当前市场上较为常见的农膜(厚度0.008 mm)。采用室内模拟土柱法, 探讨不同残膜含量对于土壤水分入渗过程的影响。试验装置由土柱、马氏瓶、支架组成。其中土柱高70 cm, 内径15 cm, 由5 mm厚的有机玻璃制成。柱体底部布设微小排气孔, 为入渗提供气流通畅的环境。马氏瓶由型号统一的玻璃瓶改装, 其内径19 cm, 高50 cm, 截面积283.53 cm²。

1.2 试验设计与方法

由于残膜70%分布在土层0~10 cm处, 30%在10~20 cm土层中, 同时随着残膜在土壤中年限的增加, 残留在土壤中的农膜会逐渐破碎, 最终碎片面积主要为3~5 cm^{2[4, 9]}。基于此, 将农膜裁剪为2 cm×2 cm的正方形, 埋设在土壤表层0~20 cm范围内, 其中残膜总量的70%埋设在土壤表层0~10 cm处, 10~20 cm埋设残膜总量的30%。根据已有研究, 当不进行机械清理残膜时, 残膜量Q_f(kg·hm^-2)与覆膜年限T(a)之间存在线性关系Q_f=5.546T+47.840(R²=0.871), 根据该方程, 为预测长期覆膜对土壤水分入渗影响, 设置5组残膜量处理:0(T0)、100(T1)、200(T2)、400(T3)、800 kg·hm^-2(T4), 分别对应不覆膜、覆膜9年、覆膜27年、覆膜63年和覆膜136年^[20]。每个处理重复3次。

装土前, 在土柱内壁均匀涂抹凡士林以消除壁面优势流的影响。按照土壤容重1.35 g·cm^-3分层(每层5 cm)填装土柱, 填装高度60 cm。为避免土壤非均质性分布及优先流产生, 将层间打毛且土柱装好后自然沉降24 h。采用一维垂直定水头法测定土壤入渗过程, 入渗时马氏瓶控制供水水头高度为5 cm。用秒表记录数据, 需要记录的指标为各个处理的湿润锋进程、入渗所需时间及马氏瓶水位变化。前5 min每10 s记录一次数据, 5~12 min每30 s记录一次数据, 12~ 30 min每1 min记录一次数据, 30~90 min每5 min记录一次数据, 90~300 min每15 min记录一次数据, 300 min以后每30 min记录一次数据, 直到湿润锋在土柱中完成运移为止, 同时记录对应时间的湿润锋位置。为了保证试验精度, 每次记录湿润锋运移位置时, 沿土柱四周测量5个深度, 取其均值。

1.3 入渗模型

为了进一步说明残膜含量对土壤水分入渗过程的影响, 选取Kostiakov模型、Philip模型、Horton模型和通用经验模型^[21]对不同残膜含量土壤入渗过程进行拟合。

(1)

式中:i(t)为土壤入渗速率, mm·min^-1; t为入渗时间, min; a、b为根据试验求得的模型参数。

(2)

式中:s为土壤吸渗率, mm·min^-0.5; t为入渗时间, min; A为稳定入渗速率, mm·min^-1。

(3)

式中:i_c为稳定入渗速率, mm·min^-1; i₀为初始入渗速率, mm·min^-1; t为入渗时间, min; k为试验求得的模型参数。

(4)

式中:a、b、-n为根据试验求得的模型参数; t为入渗时间, min。

利用均方根误差(RMSE)、平均偏差(MBE)及决定系数(R²)作为模型的评价指标。

1.4 数据处理

利用Excel 2003、SPSS 18.0软件进行数据处理和分析, 采用单因素方差分析及LSD多重比较来确定各处理入渗速率均值的差异, 采用独立样本T检验比较入渗模型模拟值和实测值的差异, 土壤入渗模型拟合使用Origin 8软件。

2 结果与分析 2.1 不同残膜含量对累积入渗量的影响

根据不同残膜含量下土壤累积入渗量随入渗历时的变化进程(图 1)可以看出, 土壤累积入渗量均随入渗的进程而增大。在同等时间内残膜含量越多, 其累积入渗量越小。累积入渗量按大小依次是:T0>T1> T2>T3>T4。进一步对累积入渗量进行差异性分析发现, 当入渗时间为0.5 min时, T0和T1之间、T3和T4之间差异不显著(P>0.05), 其他处理之间差异显著(P < 0.05)。当入渗时间为5、20、65 min时, 各残膜处理的累积入渗量之间均存在显著性差异(P < 0.05)。而当入渗时间超过270 min后, T3和T4之间差异不显著(P>0.05)。结果说明入渗初期少量残膜对累积入渗量影响不显著, 残膜量大于400 kg·hm^-2时过量的残膜对累积入渗量影响也不显著。随着入渗的进行, 不同残膜量间的累积入渗量差异愈发明显, 但到入渗后期, 过量的残膜对累积入渗量影响不显著。

2.2 不同残膜含量对湿润锋运移的影响

图 2为不同残膜含量对湿润锋运移的影响。可以看出T0~T4处理的湿润峰运移距离随着入渗的进行不断增加, 但行进至土柱同一位置所用时间不同, 当湿润锋运移至土柱底部时, 不同处理所需入渗时间为T0(390 min) < T1(420 min) < T2(510 min) < T3(600 min) < T4(660 min), 此时入渗过程完成, 其中T0和T1差异不显著(P>0.05), 其他处理之间存在显著性差异(P < 0.05)。可以看出土壤含有残膜越多, 湿润锋到达某一位置就需要更长的时间, 且少量残膜对入渗时间影响较小。湿润锋整体变化趋势与累积入渗量类似。同一入渗时间, 残膜含量越多, 湿润锋运移距离越短。在整个入渗过程的不同入渗时间, T0和T1的湿润峰运移距离基本呈现不显著差异(P>0.05), 而其他处理之间均呈现显著性差异(P < 0.05)。即土壤中少量残膜对湿润锋运移进程影响较小。

为了研究不同残膜含量下湿润锋进程的变化趋势, 根据曲线形状, 选择幂函数拟合湿润锋深度随时间的变化, 即:Z=A×t^B, 其中:Z为湿润锋深度(mm), t为入渗用时(min), A、B为拟合参数, A表征第一个计时单位后的湿润锋进程; B表征湿润锋进程的衰减程度, 拟合结果见表 1。由表 1可知, 随着残膜含量的增加, 初始湿润锋进程逐渐减小, 且差别较大, 其衰减程度随着残膜含量的增加而增大, 但变化幅度较小。其决定系数(R²)均大于0.98, 具有较好的相关性。因此, 含有残膜土壤的湿润锋深度随时间的关系遵从幂函数变化规律。

2.3 土壤入渗特征及其与残膜含量的关系

入渗速率是反映土壤水分入渗快慢的直接指标。土壤入渗过程通常可分为3个阶段:渗透初期分子力作用下入渗率较高的瞬变阶段、毛管力和重力作用下入渗率迅速下降的渗漏阶段和重力作用下的稳定入渗阶段。图 3显示不同残膜量对土壤入渗速率的影响。可以看出入渗速率随时间推移而减小并趋于平缓, 渗透初期入渗速率下降幅度大, 随着入渗进程推进下降幅度越来越小, 经过一定时间后达到稳定入渗。

利用土柱入渗前20 cm深度内入渗速率随时间的变化(图 4), 进一步分析入渗初期残膜量的影响。可以看出在入渗初期, 入渗速率随时间迅速减小, 不同残膜量处理的入渗速率大小依次为T0>T1>T2>T3> T4, 说明随着残膜量的增加, 入渗速率逐渐下降。在入渗时间达到100 min后各处理入渗速率均逐渐平缓并趋于一个稳定的值, 即稳定入渗速率, 说明此时进入稳定入渗阶段。

表 2为利用初始入渗速率、平均入渗速率、稳定入渗速率来表征土壤入渗过程。可以看出不同残膜量处理下, 土壤水分入渗各指标变化规律均为初始入渗速率>平均入渗速率>稳定入渗速率, 且各指标均随着残膜含量的增加逐渐减小:T0>T1>T2>T3>T4。与T0处理相比, 初始入渗速率分别减少了5.42%、16.75%、28.07%、35.61%, 稳定入渗速率分别减少了11.11%、22.22%、29.63%、40.74%, 平均入渗速率分别减少了14.18%、35.82%、50.00%、58.96%。其中各处理的初始入渗速率差异不显著(P>0.05), 稳定入渗速率在T0与T2、T3、T4之间存在显著性差异(P < 0.05), T1与T3、T4之间存在显著性差异(P < 0.05), T2与T4之间存在显著性差异(P < 0.05)。平均入渗速率在T0与T4之间存在显著性差异(P < 0.05)。

2.4 土壤水分入渗过程模拟

为进一步研究不同残膜量对土壤水分入渗过程的影响, 本文选用4种回归模型对土壤入渗过程进行拟合分析。由表 3可以看出, 在Philip入渗模型中参数s值为5.025~8.059, s值越大入渗能力越强, 即土壤水分入渗能力随着残膜含量的增加递减, 这与图 2中的入渗曲线变化相一致。Kostiakov模型中a值表征初始入渗速率, b值表征入渗速率的衰减程度, 由表 3可知a值随残膜含量的增加迅速降低, b值与之相反, 即随着残膜含量的增加, 土壤的初始入渗速率迅速降低, 入渗速率的衰减程度逐渐增大。Horton模型中i_c表征稳定入渗速率, 随残膜含量的增加而递减, 与表 2结果一致, 同时i_c与实际观测值仅相差0.009~0.024 mm·min^-1。但初始入渗速率i₀随残膜含量的增加而增加(除T2处理外), 与试验结果存在较大差异, 说明Horton模型不能很好反映残膜对土壤入渗过程的影响。通用经验模型中a和b的绝对值分别表征稳定入渗速率和初始入渗速率, 均表现为T0>T1>T2>T3> T4, 即随土壤中残膜含量增加, 初始入渗速率和稳定入渗速率均呈递减趋势, 与实测值一致。表 3还可看出通用经验模型R²值最高, 介于0.958~0.967, RMSE最小, MBE最接近于0, 拟合效果最好。Kostiakov模型次之, 其R²值在0.916~0.950, Philip模型和Horton模型R²值均小于0.922, RMSE较大, MBE与0相距较远, 拟合效果较差。综上可以看出通用经验模型能很好拟合不同残膜含量下淡灰钙土水分入渗过程。

为进一步验证模型拟合结果, 对入渗速率的模型模拟值和实测值进行独立样本T检验。如果两组样本数据t值的绝对值小于T_0.05且P值大于0.05, 则两组数据没有显著差异。由表 3可以看出:Philip入渗模型、Kostiakov模型和通用经验模型的t值绝对值小于T_0.05且P值大于0.05, 说明这3种模型模拟不同残膜含量下淡灰钙土水分入渗过程是可行的。而Horton模型的检验结果表明, T3和T4处理下t值的绝对值大于T_0.05且P值小于0.05, 表明两组数据具有显著性差异, Hor- ton模型不适合模拟含残膜淡灰钙土水分入渗过程。

由于Horton模型R²值远低于其他3种模型, 独立样本T检验也说明Horton模型不适合模拟含残膜淡灰钙土水分入渗过程, 而Philip模型、Kostiakov模型和通用经验模型的R²值均高于0.88, 因此将入渗试验实测的T4处理累积入渗量和入渗模型计算值进行对比(图 5), 分析这3种模型对含残膜土壤水分入渗过程的适应性。可以看出实测累积入渗量与3种模型拟合的数据在入渗初期差异较小, 随入渗进行差异逐渐增大, 在入渗后期Philip模型和Kostiakov模型与实测值差距明显, 而通用经验模型与实测值的差异逐渐减小。进一步对数据进行独立样本T检验, 发现实测值与Philip模型和Kostiakov模型获得的累积入渗量之间存在显著差异(P < 0.05), 而与通用经验模型没有显著差异(P=0.248), 说明通用经验模型描述含残膜的淡灰钙土水分入渗过程是可行的。

3 讨论

土壤中农用残膜会显著影响淡灰钙土水分入渗过程, 不同残膜含量下土壤入渗曲线随时间的延长先急剧降低后渐趋平缓^{[15, 22]}。不同残膜处理对淡灰钙土水分初始入渗速率差异不显著(P>0.05), 这是因为在入渗初期, 土壤表面干燥, 湿润锋锋面的水势梯度和非饱和度大, 入渗速率快, 不同残膜量处理对湿润锋和入渗量的影响较小。随着入渗进程推进, 当土壤含水率达到最大分子持水量时, 残膜的阻水效应显现, 且随残膜量增多, 其阻水效应增强, 因此不同残膜量处理对淡灰钙土水分湿润锋和入渗量的影响显著(P < 0.05)。随着水分的继续下渗, 土壤趋于饱和, 达到稳定入渗阶段, 残膜的阻水效应减弱, 不同处理入渗速率差异减小, 最后逐渐趋于相近。本研究得出结论与其他不同土壤质地含残膜土壤水分入渗趋势基本一致^[15-16]。但同时发现不同质地土壤水分具体入渗过程存在明显区别。本研究发现淡灰钙土入渗约5 min后不同残膜量处理间入渗率出现显著性差异(P < 0.05), 且随入渗过程的进行差异逐渐变大, 在入渗时间达到100 min后入渗速率逐渐趋于平缓, 稳定入渗速率为0.16~0.27 mm·min^-1; 砂壤土在入渗约100 min后, 不同残膜量处理出现显著差异(P < 0.05), 360 min后达到稳定入渗阶段, 稳定入渗速率约为0.07 mm·min^-1[14]; 砂土入渗约3 min后就出现显著差异, 15 min后达到稳定入渗阶段, 稳定入渗速率约为1.8 mm·min^-1[14]; 粉砂质黏壤土在入渗约360 min后, 不同残膜量处理之间才出现显著差异^[23]。相关研究表明土壤质地是影响含残膜土壤水分入渗过程的重要因素^{[20, 23]}。银川平原淡灰钙土砂粒含量高, 有机质含量低, 大孔隙多, 因此其入渗过程和入渗率介于砂土和粉砂质黏壤土之间。

少量残膜对淡灰钙土水分入渗初期湿润锋运移距离和累积入渗量影响不显著(P>0.05), 而随残膜含量增加, 其对土壤入渗的阻滞作用逐渐增强。残膜存在破坏了土壤孔隙结构、阻断毛管孔隙及堵塞大孔隙, 减小了过水断面积, 增加了水流弯曲度, 抑制土壤水分下渗。王志超等^[24]通过CT扫描技术分析土壤中残膜面积、存在位置和分布形态后, 发现残膜存在破坏了土壤孔隙连续性, 造成了土壤中的大孔隙堵塞, 阻碍了土壤水分运移。本研究发现T3(400 kg·hm^-2)和T4(800 kg·hm^-2)处理间累积入渗量差异不显著(P>0.05), 即在高残膜含量条件下, 随着残膜含量增加其累积入渗量的增幅在统计学意义上不显著(P>0.05)。分析原因可能是当土壤残膜量过多时, 土壤大孔隙的比例增大, 甚至可能产生优势流, 导致土壤水分的优势迁移, 使水分运移距离增加。李元桥等^[15]也发现当土壤残膜量达到720 kg·hm^-2时, 产生优势流, 导致湿润峰的运移加快。综上可以看出土壤中残膜含量的增加, 一方面降低土壤孔隙度和通透性, 阻断土壤孔隙连续性, 抑制土壤水分下渗; 另一方面土壤残膜量过多时可能会产生优势流, 导致残膜量达到400 kg·hm^-2以上时, 入渗量差异不显著(P>0.05)。因此有必要利用土壤切片或CT扫描等技术深入分析残膜对土壤孔隙结构(孔隙度、孔隙数量、成圆率及分形维数等)的影响, 以便能通过土壤孔隙特征变化, 揭示不同残膜量对土壤入渗的影响机理。

Horton模型拟合得出T0~T4处理i₀依次为3.69、7.701、5.349、11.392、13.713 mm·min^-1, 即除T2处理外, i₀随残膜含量增加而增加, 而试验结果得出i₀随残膜含量增加而减少, 拟合结果和实测结果不一致, 且拟合效果不理想(R²值为0.515~0.922), 说明Horton模型不能很好反映含残膜土壤入渗过程, 相关学者在进行含残膜土壤入渗模型拟合时也都没有选择Hor- ton模型^{[11, 16, 20]}。Philip模型拟合R²值为0.880~0.919, 低于Kostiakov模型和通用经验模型, 邹小阳等^[23]也发现Philip模型对参数的精度和灵敏性要求更高, 且在模拟长历时入渗过程时, Philip入渗模型偏差较大。Kostiakov模型和通用经验模型对不同残膜含量土壤水分入渗拟合均较好, 但是Kostiakov模型中当t趋于0时, i(t)趋于无限大; 而t趋于无限大时, i(t)趋于0, 这与实际状况不符。通用经验模型添加了参数a, 克服了Kostiakov模型的不足, 即t趋于无限大时, i(t)达到稳定入渗率^[25]。

4 结论

(1) 随残膜含量递增, 在同一入渗历时内累积入渗量逐渐减小。入渗初期100 kg·hm^-2残膜量对累积入渗量影响不明显, 随入渗过程的进行不同残膜量处理的累积入渗量差异变大, 入渗后期过量的残膜对累积入渗量影响不明显。

(2) 残膜含量越少, 湿润锋进程越快。100 kg·hm^-2残膜量对湿润锋运移距离影响不显著。含残膜土壤的湿润锋深度随时间遵从幂函数变化规律。

(3) 不同残膜含量土壤入渗曲线随时间先急剧降低后渐趋平缓, 土壤水分入渗各指标变化规律均为初始入渗速率>平均入渗速率>稳定入渗速率。不同残膜量处理对初始入渗量影响不明显, 对稳定入渗速率和平均入渗速率影响明显。各入渗速率均随着残膜含量的增加逐渐减小, 说明随残膜量增多, 其阻水效应增强。

(4) 通用经验模型更适合描述含残膜土壤的水分入渗过程。