考虑成层土阻水减渗效应的地表水入渗过程分析模型

张兴毅,曹文贵,崔鹏陆,李慧鑫

(湖南大学 土木工程学院,长沙 410082)

1 研究背景

地表土由于沉积作用而普遍具有成层结构,由于不同土层的土壤特性不同,因此具有不同的基质吸力[1-2]。基质吸力是非饱和土壤水分运动的驱动力[3],致使水从高吸力土层进入低吸力土层时,由于吸力不足,下层土会对下渗水起到阻水的作用[4],从而对成层土入渗过程产生重要的影响。如何考虑成层结构对地表水入渗的水分运移过程的影响正是本文研究的出发点。

在一维常水头入渗中,湿润土与未被湿润的原始土的分界面称为湿润锋,由于湿润锋处土壤含水量梯度大,因此湿润锋处的土壤吸力作用明显,对入渗影响很大[5]。研究成层土水分入渗过程即研究湿润锋与入渗时间的关系。国内外学者建立了许多描述水分入渗过程的模型[6-9]。其中由Green和Ampt[9]于1911年基于毛细管理论提出的Green-Ampt土壤水分入渗模型(GA模型),也称传统GA模型,如图1所示,虽然有物理意义明确、形式简单、模型参数少等优点,但由于GA模型是假设湿润锋饱和并基于均质土一维常水头入渗提出的,因此国内外学者对成层土的入渗规律进行了更进一步的研究:王文焰等[1]通过室内试验发现当湿润锋到达界面时,湿润锋为非饱和且吸力大于下层砂土的基质吸力,湿润锋下移受阻发生停滞现象;Fok[10]假定损失水头为上下土层之和,并根据各土层厚度推求平均导水系数,以此计算土壤入渗情况;王文焰等[11]根据土壤水的能量原理,研究了夹砂层黄土的入渗特性,发现砂层使得入渗受阻;王全九等[12]根据GA模型的特性,分析了GA模型应用于层状土可行的计算方法及困难;范严伟等[13]利用Hydrus-1D模拟了不同情况下入渗过程,通过导水系数改进了GA模型,提出了H-GA模型;王春颖[14]通过层状夹砂土柱一维室内试验,进一步改进H-GA模型,提出了有针对性的S-GA模型。

图1 Green-Ampt模型简图Fig.1 Brief description of the Green-Ampt model

有学者通过均质土一维入渗试验发现湿润锋为非饱和状态,并且针对含水量随深度变化曲线进行了一定的改进[15-17]。因此基于假设湿润锋饱和提出的H-GA模型、S-GA模型等成层土入渗模型中,入渗进行到分界面时只有减渗作用,没有停滞受阻现象,模型与实际入渗情况不符。综上,研究成层土入渗时湿润锋饱和情况及其在入渗过程中变化的影响具有重要的理论意义。

针对上述已有研究的不足,本文以地表水入渗成层土为研究对象,考虑入渗时湿润锋非饱和及其与阻水减渗效应的关系,提出成层土入渗全过程分析模型并进行深入研究,以期完善地表水入渗过程分析的理论方法。

2 考虑成层土阻水减渗的地表水入渗模型

在研究地表水入渗时,湿润锋是否饱和是研究的关键。结合Richards模型[8]的模拟结果和文献[16]中利用VG-M模型[18]计算土体含水量的结果,发现当采用椭圆曲线拟合含水量分布时,与Richards方程计算结果接近,并取得了不错的效果。因此本文沿用此假设:入渗过程中,在含水量饱和的饱和区与含水量为初始状态的原始区之间存在一个含水量随深度呈1/4椭圆弧线变化的非饱和过渡区,得到了一个分层含水量分布模型,见图2。该假设与Coleman和Bodman[19]提出的分层思想是一致的。

图2 湿润锋分层简图Fig.2 Schematic diagram of the wetting front stratification

为了建立新的地表水入渗模型,便于本文研究,对实际入渗进行简化,如图3所示。并作出如下假设:

图3 地表水一维入渗模型简图Fig.3 Schematic diagram of one-dimensional infiltration model of surface water

(1)地表有一层深度不变且均匀的积水,深度为H;黏性土初始含水量均为θd,l;砂土初始含水量均为θd,s,砂土埋深为Z1。

(2)上层黏性土吸力大于下层砂土吸力。

(3)湿润锋存在非饱和过渡区,且过渡区含水量随深度变化呈1/4椭圆弧线变化。

当湿润锋到达分界面前,土体含水量情况如图4所示,因此土体总入渗量为饱和区与过渡区含水量之和,即

图4 湿润锋到达分界面前土体含水量剖面Fig.4 Profile of soil water content before the wetting front reaches the interface

式中:I为累计入渗量(cm);θs,l为饱和体积含水量;Zf为湿润锋深度(cm);Zd为饱和区深度(cm)。

饱和区深度与湿润锋深度呈一定比例关系[16],因此则有

Zd=γZf。

(2)

式中γ为饱和区占总湿润锋深度的比例。

根据GA模型[9],入渗速率为

(3)

式中:i(t)为入渗速率(cm/min);t为入渗时间(min);ks,l为黏性土的饱和渗透系数(cm/min);hf,l为黏性土湿润锋基质吸力(cm)。

联立式(1)和式(3)并积分可得湿润锋与时间的关系为

其中,

当湿润锋到达分界面时,如图5实线部分,Zf=Z1,入渗历经时间t1为

图5 湿润锋到达分界面时土体含水量剖面Fig.5 Profile of soil water content when the wetting front reaches the interface

其中,

(7)

此时的入渗速率为

(8)

当湿润锋在分界面处入渗时,由于黏性土吸力更高,因此湿润锋行进发生停滞。但入渗持续进行,湿润锋吸力逐渐减小,湿润锋进入砂层时,上层黏性土接近饱和[1,14]。该变化即为图5中实现部分向虚线变化的过程,假设在该阶段入渗速率线性降低。累计入渗变化量ΔI是线性变化的入渗速率对时间的积分,则滞水时间内可以表示为

(9)

式中:ΔI为湿润锋由非饱和到饱和过程中增加的水量(cm);ist为最终稳渗率(cm/min)。

化简积分后可得此过程所需时间Δt为

(10)

此过程中累计入渗量及入渗速率变化分别为:

I=I1+(i1+i)(t-t1)/2 ;

(11)

(12)

其中,

(13)

当t=t2时,湿润锋接近饱和,吸力接近最小,认为此时湿润锋吸力与下层土进水吸力接近,湿润锋进入砂层,t2表达式为

t2=t1+Δt。

(14)

此时,Zf≈Z1,下层砂土仍处在初始状态。Colman和Bodman[19]认为成层土壤的入渗特性应由细质土来控制,根据水量平衡原理,上下界面入渗速率应相同,此时的入渗速率为

(15)

式中:kst,s为砂土的渗透系数(cm/min);hf,s为砂土湿润锋基质吸力(cm)。

又由于湿润锋经过界面时,上下土层此时的吸力值很小,则将式(15)化简得

(16)

通过计算可知,此时砂土导水系数并没有达到饱和导水系数,砂土仍然处于非饱和状态。若表面积水深度可忽略的话,所得入渗速率接近上层土的饱和导水率,符合Hill和Parlange[20]的试验结果。

因此考虑砂土基质吸力的稳定入渗速率(稳渗率)表达式应为

(17)

式(17)所得的稳渗率即此时砂土的非饱和渗透系数,利用VG-M模型[18]可反求得到砂土此时的含水量为

(18)

其中,

(19)

式中:ks为饱和渗透系数;Se为砂土的相对饱和度;θst为砂土稳渗时的含水量;θr为残余含水量;m、l均为VG-M模型参数。

湿润锋在砂土中的行进仍然可以视为非饱和区与过渡区的组合,其含水量分布如图6所示。此时的总入渗量为:

图6 湿润锋穿越分界面后土体含水量剖面Fig.6 Profile of soil water content after the wetting front crosses the interface

I=Il+Is;

(20)

Il=(θs,l-θi,l)Z1;

(21)

(22)

由于此时稳渗率为常数,联立式(14)、式(17)、式(22),可知地表水入渗模型湿润锋深度与时间的变化关系为

(23)

因此地表水入渗模型如表1所示。

表1 地表水入渗模型Table 1 Surface water infiltration model

3 模型结果与验证

王春颖等[14]通过试验,在传统GA模型上建立了夹砂层一维入渗模型,简称为S-GA模型。为验证模型可靠性,取文献[14]中的试验数据和S-GA模型的模拟结果,对本文的地表水入渗模型进行验证。建立土壤深度为42.5 cm的模型,积水深度为2 cm,其中模型上层为黏性土,下层为砂土,砂土在地下22.5 cm处,土壤特性参数见表2。

表2 模型参数Table 2 Model parameters

根据王文焰等[1]研究中甘肃秦王川地区实验资料,饱和区与过渡区深度相同,γ=0.5。本文将从累计入渗量、入渗深度及入渗速率变化3个方面对模型进行对比验证。

利用最小均方偏差(RMSD)评价模拟的结果,RMSD可表示为

(24)

式中:Yi计算表示通过模型计算得到的参数;Yi试验表示试验数据测量的相关结果。

RMSD计算结果见表3,模型模拟结果与试验数据的对比如下文所示。

表3 RMSD计算结果Table 3 RMSD calculation results

3.1 累计入渗量变化

试验与模拟的成层土柱的累计入渗量随时间变化规律如图7所示。入渗试验进行到55 min左右,累计入渗量随时间变化斜率减小,并且最终接近常数。S-GA模型未考虑下层砂土的基质吸力影响,因此地表水入渗模型斜率大于S-GA模型,导致地表水入渗模型累计入渗量更接近试验结果。从斜率变化上来看,本文模型斜率是55 min开始减小,62 min时减小到固定值,而S-GA模型则是55 min时突变为固定值,这是由于本文模型假设湿润锋非饱和,因此模拟了湿润锋在分界面停滞时的入渗速率减小现象。本文模型RMSD为0.141,小于S-GA模型的0.217。因此本文模型更准确地描述了成层土入渗过程中累计入渗量随时间的变化。

图7 实测与模拟累计入渗量随时间的变化Fig.7 Measured and simulated cumulative infiltration with time

3.2 入渗深度变化

试验与模拟的成层土柱的入渗深度随时间变化规律如图8所示。湿润锋在55 min左右时,存在明显的行进停滞现象,直至62 min左右行进继续。S-GA模型在模拟入渗过程中没有停滞受阻现象,仅仅在62 min时斜率发生了突变。而本文模型很好地模拟了行进停滞的现象。这是由于本文模型假设湿润锋非饱和,且下层土吸力低于上层土,因此当湿润锋到达分界面时,由于下层土吸力不足,湿润锋行进发生停滞。湿润锋到达分界面前,地表水入渗模型入渗速度更快,这是由于本文模型假设湿润锋为非饱和,湿润锋在相同入渗速率时入渗深度更深。本文模型RMSD为0.916,小于S-GA模型的1.403。地表水入渗模型在模拟出阻水作用的情况下,与试验数据吻合良好,能够很好地描述成层土入渗时深度与时间的关系。

图8 实测与模拟湿润锋位置随时间的变化Fig.8 Measured and simulated advance of the wetting front with time

3.3 入渗速率变化

试验与模拟的成层土柱的入渗速率随时间的变化如图9所示。在55～62 min时,入渗速率呈现线性下降趋势,62 min后入渗速率趋于稳定。S-GA模型在55 min时,入渗速率发生突变,但本文模型入渗速率发生线性变化,符合实际变化情况。这是因为S-GA模型没有考虑非饱和区的存在,因此当湿润锋到达分界面时,湿润锋直接穿越界面使得入渗速率发生了突变。试验数据中稳渗率为0.037 7,地表水入渗模型的最终稳渗率为0.037 8,而S-GA模型为0.037 0,这是由于地表水入渗模型考虑了砂土基质吸力,因此稳渗率偏大。本文模型RMSD比S-GA模型的更小,为0.037。综上可以认为地表水入渗模型更加符合实际情况,有较为理想的可靠性和准确度。

图9 入渗速率随时间变化Fig.9 Change of infiltration rate with time

从以上结果可知,两个模型均能描述成层土中的入渗规律和砂层减渗的特性,但是地表水入渗模型能够更好地描述湿润锋到达分界面时的停滞现象,并且地表水入渗模型模拟了停滞过程中入渗速率线性降低的变化趋势,符合实际情况。因此本文提出的模型在模拟成层土入渗中更具优势。

4 相关因素对地表水入渗模型的影响

Hydrus-1D是国际地下水模型中心公布的,利用一维饱和-非饱和带水分运移基本方程(Richards方程)计算水分、盐分运移规律的软件,其内核为一维的有限元模型。为了使本文模型具有一定的普适性,本文基于范严伟等[13]利用Hydrus-1D研究夹砂层土壤入渗的模拟结果,对比Hydrus-1D模拟结果与本文模型计算结果,进一步研究本文模型在各个影响因素条件下的入渗规律,分析不同因素对模型的影响。模拟中的土壤Van Genuchten-Mualem(简称VG-M)模型水力特性参数分别取自参考文献[13]、文献[21]和文献[22],如表4所示。

表4 不同土质VG-M模型水力特性参数Table 4 Hydraulic parameters in Van Genuchten-Mualem model of different soils

4.1 初始含水量的影响

采用S1下层砂土和A1上层黏性土,上层土厚30 cm,下层土厚15 cm,积水深度为6 cm,计算得出不同上层土壤初始含水量时入渗特性曲线,数据点为Hydrus-1D模拟结果,如图10所示。

图10 初始含水量对模拟结果的影响Fig.10 Influence of initial moisture content on simulation results

当初始含水量增加时,累计入渗量略有减少,但入渗速度加快。这是由于初始含水量的增加减小了湿润锋基质吸力,降低了入渗速率。入渗完成时间随初始含水量的增大而减小,这是由于初始含水量的增加使得土体更易达到饱和。初始含水量对滞水时间影响很小,对最终的稳渗率没有影响。综上可知初始含水量对积水入渗的影响较小,主要影响上层土中入渗速率的变化,从而影响入渗时间。本文模型与模拟结果吻合良好,证明了模型的合理性。

4.2 积水深度的影响

取S1下层砂土和A1上层黏性土,上层土厚30 cm,下层土厚15 cm,计算不同积水深度时入渗特性曲线,数据点为Hydrus-1D模拟结果,结果如图11所示。

图11 积水深度对模拟结果的影响Fig.11 Influence of water depth on simulation results

当积水深度增大,土壤的累计入渗量不变,入渗速度加快,入渗时间减少。这是由于积水深度增加使入渗压力水头增大,根据达西定律,入渗速率增大。水头差增大使停滞时间略微减小,但进入下层砂土的水量增加,因此砂土稳渗时含水量增加,稳渗率增大。积水10 cm时,出现较为明显的误差,这是由于Hydrus-1D模拟时,黏性土与砂土初始含水量不同,在分界面处存在一定吸力差,使分界面处发生水量交换,使上下界面吸力接近,因此滞水作用减弱,湿润锋提前进入砂层,较计算结果偏大,水头差较大时,这种影响也较大。因此在积水较浅时地表水入渗模型有更好的效果。

4.3 土层埋深的影响

取S1砂土和A1黏性土,分别取上层黏性土厚10 cm与30 cm,砂层厚15 cm,积水深度为6 cm,计算不同埋深时入渗特性曲线,数据点为Hydrus-1D模拟结果,如图12所示。

图12 埋深对模拟结果的影响Fig.12 Influence of burial depth on simulation results

砂层埋深越深,滞水作用越明显,但减渗效果越小。这是因为砂层埋深越深,非饱和区域增大,因此湿润锋在分界面到达饱和的时间增加,滞水作用更加明显。埋深增加使湿润锋更迟到达分界面,到达界面时入渗速率更小,更接近黏性土的饱和导水率,因此减渗效果降低。所以较浅的砂层埋深,在短期内能快速减渗,但无法阻滞水分下渗,从长期来看,稳渗率反而会更大。因此砂层埋深对地表水入渗模型结果影响较大,埋深直接影响了阻水减渗作用的效果,合适的埋深对于工程中阻水减渗有着重大意义。

4.4 土层种类的影响

分别取S1、S2、S3下层砂土和A1上层黏性土,上层土厚30 cm,下层土厚15 cm,积水深度为6 cm,计算不同砂土种类时双层土壤入渗特性曲线,数据点为Hydrus-1D模拟结果,如图13所示。

图13 砂层种类对模拟结果的影响Fig.13 Influence of sand type on simulation results

当湿润锋进入下层土后,不同砂土的稳渗率不同,但是相差不大,并且累计入渗量也各不相同。这是由于不同砂土的基质吸力不同,因此稳渗时砂土的含水量不同,稳渗率也不相同。但相对于水头(即黏性土厚度与积水深度之和)而言,砂土的基质吸力很小,因此对稳渗率的影响很小。残余含水量和饱和含水量的差异,也使得不同砂土间最终入渗时间有较大的差别。因此不同的低吸力下层土对入渗有不同的影响,地表水入渗模型很好地体现了这一点。

4.5 湿润锋过渡区占比的影响

采用S1下层砂土和A1上层黏性土,上层土厚30 cm,下层土厚15 cm,积水深度为6 cm,计算得出不同湿润锋过渡区占比时入渗特性曲线,如图14所示。

图14 湿润锋过渡区占比对模拟结果的影响Fig.14 Influence of transition zone proportion on simulation results

在其他条件相同时,随着湿润锋过渡区占比减小,湿润锋行进需要更多入渗量,因此到达分界面的速度减慢,并且阻水时间有明显的减少,入渗速率开始改变的时间也延后。但累计入渗量和稳渗率没有变化。因此在使用地表水入渗模型进行模拟预测时,选择合理的湿润锋过渡区占比对模型的精确度有重要的意义。

5 结 论

本文以地表水入渗的成层土为研究对象,假设成层土入渗时湿润锋含水量为非饱和,提出了一种适用于从高吸力土到低吸力土的地表水入渗模型。结合前人的夹砂层土柱一维入渗试验和Hydrus-1D模拟结果分析研究,得出以下结论:

(1) 通过与试验数据及S-GA模型数值解对比,地表水入渗模型更好地模拟了阻水减渗作用,验证了本文模型在求解成层土地表水入渗规律方面的正确性与适用性。

(2) 成层土入渗过程中,当湿润锋行进发生停滞时,入渗速率呈线性下降的趋势,地表水入渗模型很好地体现了这个入渗特征。

(3)稳渗率不仅受上层土入渗特性控制,不同下层低吸力土的稳渗率不同,滞水能力也不同,下层土对土壤入渗有较大影响。

(4)地表水入渗中,下层低吸力土的种类、埋深及积水深度对入渗特性有较大的影响;土体初始含水量仅对入渗时间影响较大;湿润锋过渡区占比对滞水作用有较大影响。