包气带井回灌引起的非饱和-饱和流分析.pdf

第4 2卷 第4期2 0 2 3年 7月 地 质 科 技 通 报B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y V o l.4 2 N o.4J u l.2 0 2 3祁翠婷,詹红兵,郝永红.包气带井回灌引起的非饱和-饱和流分析J.地质科技通报,2 0 2 3,4 2(4):1 1 8-1 2 9.Q i C u i t i n g,Z h a n H o n g b i n,H a o Y o n g h o n g.A n a l y s i s o f u n s a t u r a t e d-s a t u r a t e d f l o w i n d u c e d b y a v a d o s e z o n e w e l l i n j e c t i o nJ.B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y,2 0 2 3,4 2(4):1 1 8-1 2 9.基金项目:国家自然科学基金地质联合基金项目(U 2 2 4 4 2 1 4)作者简介:祁翠婷(1 9 9 2),女,讲师,主要从事含水层人工补给及含水层热能储存模拟研究。E-m a i l:q i c u i t i n g 1 2 6.c o m E d i t o r i a l O f f i c e o f B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y.T h i s i s a n o p e n a c c e s s a r t i c l e u n d e r t h e C C B Y-N C-N D l i c e n s e.包气带井回灌引起的非饱和-饱和流分析祁翠婷1 a,詹红兵2,郝永红1 b(天津师范大学 1 a.地理与环境科学学院;1 b.天津市水资源与水环境重点实验室,天津 3 0 0 3 8 7;2.美国德克萨斯农工大学地质与地球物理系,德克萨斯 7 7 8 4 3-3 1 1 5)摘 要:包气带井回灌是人工补给含水层的重要方法。在回灌过程中,非饱和带水力参数的准确描述对于提高评价精度、增强回灌管理具有重要意义。目前已有的许多包气带井回灌的解析模型是基于G a r d n e r土壤水分特征曲线模型(含二参数)建立的,随着三参数模型(MB模型)和四参数模型(MN模型)的提出,一个值得关注的问题是使用更加灵活的土壤水分特征曲线模型是否能提高回灌系统中非饱和-饱和流的模拟精度。本研究采用MN模型(G a r d n e r模型和MB模型为其子集)建立了包气带井回灌模型,使用C OM S O L M u l t i p h y s i c s对模型进行求解,所得解用于研究在回灌过程中非饱和带水力参数和地表通量对非饱和-饱和流的影响,并比较在回灌过程中基于不同土壤水分特征曲线模型所得到的水力响应。试验结果分析发现包气带井回灌引起的水力响应以及地表通量对回灌的影响均会受到非饱和带渗透性以及储水能力的影响,相对渗透系数指数k影响非饱和带渗透系数的变化,持水指数c影响非饱和带的储水能力,将指数参数简化为k=c=会在计算和预测包气带井回灌引起的水力响应时带来一定的误差。压力头阈值差b1=a-k的绝对值较小时,对水头增量的影响也较小,此时将其简化为b1=a-k=0所带来的误差也较小。本研究成果可以帮助学者们提高对包气带井回灌过程的认识,对于回灌方案设计、实施和管理具有重要的实际意义。关键词:包气带井;人工回灌;非饱和-饱和流;土壤水分特征曲线模型2 0 2 2-1 2-2 2收稿;2 0 2 3-0 3-1 6修回;2 0 2 3-0 4-1 2接受中图分类号:P 6 4 1 文章编号:2 0 9 6-8 5 2 3(2 0 2 3)0 4-0 1 1 8-1 2d o i:1 0.1 9 5 0 9/j.c n k i.d z k q.t b 2 0 2 2 0 7 0 3 开放科学(资源服务)标识码(O S I D):A n a l y s i s o f u n s a t u r a t e d-s a t u r a t e d f l o w i n d u c e d b y a v a d o s e z o n e w e l l i n j e c t i o nQ i C u i t i n g1 a,Z h a n H o n g b i n2,H a o Y o n g h o n g1 b(1 a.S c h o o l o f G e o g r a p h i c a n d E n v i r o n m e n t a l S c i e n c e s;1 b.T i a n j i n K e y L a b o f W a t e r R e s o u r c e s a n d W a t e r E n v i r o n m e n t,T i a n j i n N o r m a l U n i v e r s i t y,T i a n j i n 3 0 0 3 8 7,C h i n a;2.D e p a r t m e n t o f G e o l o g y a n d G e o p h y s i c s,T e x a s A&M U n i v e r s i t y,T X 7 7 8 4 3-3 1 1 5,U S A)A b s t r a c t:O b j e c t i v eV a d o s e z o n e w e l l(V ZW)i n j e c t i o n i s a n e f f e c t i v e m e t h o d o f m a n a g e d a q u i f e r r e c h a r g e(MA R).T o i m p r o v e V ZW i n j e c t i o n m a n a g e m e n t,i t i s o f g r e a t i m p o r t a n c e t o a c c u r a t e l y d e s c r i b e t h e u n-s a t u r a t e d z o n e p r o p e r t i e s.S e v e r a l a n a l y t i c a l m o d e l s h a v e b e e n d e v e l o p e d f o r V ZW i n j e c t i o n b a s e d o n t h e t w o-p a r a m e t e r c o n s t i t u t i v e m o d e l(G a r d n e r m o d e l).A s t h e t h r e e-p a r a m e t e r m o d e l(MB m o d e l)a n d f o u r-p a r a m e t e r m o d e l(MN m o d e l)h a v e b e e n p r o p o s e d,i t i s o f i n t e r e s t t o k n o w w h e t h e r t h e a p p l i c a t i o n o f m o r e f l e x i b l e c o n s t i t u t i v e m o d e l s i s a b l e t o i m p r o v e t h e a n a l y s i s o f u n s a t u r a t e d-s a t u r a t e d f l o w i n d u c e d b y 第4期祁翠婷等:包气带井回灌引起的非饱和-饱和流分析V ZW i n j e c t i o n.M e t h o d sI n t h i s s t u d y,t h e MN m o d e l(w h i c h i n c l u d e s t h e G a r d n e r m o d e l a n d MB m o d e l a s s u b s e t s)w a s e m p l o y e d t o e s t a b l i s h a V ZW i n j e c t i o n m o d e l.T h e m o d e l w a s s o l v e d u s i n g C OM S O L M u l t i p h y s i c s.T h e r e s u l t s a r e u t i l i z e d t o i n v e s t i g a t e t h e u n s a t u r a t e d-s a t u r a t e d f l o w i n d u c e d b y V ZW i n j e c-t i o n f o r d i f f e r e n t v a l u e s o f t h e u n s a t u r a t e d z o n e p r o p e r t i e s t o a n a l y s e t h e i n f l u e n c e s o f t h e g r o u n d s u r f a c e f l u x(G S F),a n d t o c o m p a r e t h e h y d r a u l i c r e s p o n s e s b a s e d o n d i f f e r e n t e x p o n e n t i a l c o n s t i t u t i v e m o d e l s.R e s u l t sT h e a n a l y s i s d e m o n s t r a t e s t h a t t h e h y d r a u l i c r e s p o n s e i n d u c e d b y V ZW i n j e c t i o n a n d t h e i n f l u-e n c e o f g r o u n d s u r f a c e f l u x a r e a f f e c t e d b y t h e h y d r a u l i c c o n d u c t i v i t y a n d t h e w a t e r s t o r a g e c a p a c i t y o f t h e u n s a t u r a t e d z o n e.T h e r e l a t i v e h y d r a u l i c c o n d u c t i v i t y e x p o n e n t k a f f e c t s t h e c h a n g e i n h y d r a u l i c c o n d u c-t i v i t y o f t h e u n s a t u r a t e d z o n e.T h e m o i s t u r e r e t e n t i o n e x p o n e n t c a f f e c t s t h e w a t e r s t o r a g e c a p a c i t y o f t h e u n s a t u r a t e d z o n e.T h e a p p r o x i m a t i o n o f k=c=w i l l r e s u l t i n s o m e e r r o r s i n t h e c a l c u l a t i o n a n d p r e d i c-t i o n o f t h e h y d r a u l i c r e s p o n s e c a u s e d b y V ZW i n j e c t i o n.Wh e n t h e a b s o l u t e v a l u e d i f f e r e n c e o f t h e p r e s s u r e h e a d t h r e s h o l d b1=a-k i s s m a l l,i t s e f f e c t o n t h e w a t e r h e a d i n c r e m e n t i s s m a l l.I n t h i s c a s e,i t i s r e a-s o n a b l e t o a s s u m e t h a t b1=a-k=0.C o n c l u s i o nT h i s s t u d y c a n h e l p s c h o l a r s i m p r o v e t h e u n d e r s t a n d-i n g o f t h e V ZW i n j e c t i o n p r o c e s s a n d h a s i m p o r t a n t p r a c t i c a l s i g n i f i c a n c e f o r t h e d e s i g n,i m p l e m e n t a t i o n a n d m a n a g e m e n t o f i n j e c t i o n s c h e m e s.K e y w o r d s:v a d o s e z o n e w e l l;m a n a g e d a q u i f e r r e c h a r g e;c o u p l e d u n s a t u r a t e d-s a t u r a t e d f l o w;s o i l c o n s t i t u-t i v e m o d e lR e c e i v e d:2 0 2 2-1 2-2 2;R e v i s e d:2 0 2 3-0 3-1 6;A c c e p t e d:2 0 2 3-0 4-1 2 含水层人工回灌(m a n a g e d a q u i f e r r e c h a r g e,简称MA R)指人们主动将水回灌入地下含水层中储存起来以备后续使用,或用来补给含水层并解决一些环境问题1。MA R是一项跨多个领域的技术,它可以用于保障水资源供给2-4、管理和控制地下水水质与水量5-7、保护和维持依赖地下水的生态系统8-1 0。随着世界人口的日益增长,人们对于水资源的需求量也在不断增大,同时人们也需要开发新空间并增强蓄水能力,在雨季等水源丰富的时期将水储存起来用于满足干旱时期的用水需求,这将使MA R的 应 用 与 推 广 变 得 愈 加 重 要1 1-1 6。B o u w-e r1 1总结了常见的MA R系统,包括地表回灌系统、包气带回灌系统以及组合回灌系统。包气带回灌技术是替代MA R地表回灌技术的一种很有前景的方法,它具有减少回灌过程中水的蒸发损失量和可利用地下高渗层提高回灌效率等优点1 7。包气带井(v a d o s e z o n e w e l l)也被称为干井,是在包气带中钻孔用于补给含水层的井,它使得回灌能够开始于地下较深层从而避开渗透性较差的地表层。包气带井回灌是在地下水位较深、缺乏渗透性较好的土壤层的半干旱和干旱地区,以及没有足够大的土地面积进行地表回灌的地区实施人工补给含水层的一项重要技术1 1。目前包气带井回灌已在全球各国得到了广泛使用,包括美国1 8-2 1、中国2 2-2 3、奥地利2 4、芬兰1,2 5、葡萄牙2 6等。为了提高对包气带井回灌的评估与管理,人们需要加强对回灌过程中非饱和带流动过程的认识,并能够较为准确地描述非饱和带水力特性。为了评估回灌效果,一方面学者们基于有限元法、有限差分法及有限体积法等提出了各种数值模型,另一方面,解析模型也是回灌评价的重要方法之一。尽管解析模型使用时做出的一些简化假设使其不适用于非均质含水层回灌等复杂情况的分析,但是解析模型具有应用简单、易于实现、计算省时等优点,使得解析模型仍是回灌评价的有力工具。在早期研究中,用于解释非承压含水层试验的解析模型都假定非饱和带的影响可以忽略不计2 7。后来,Nw a n k w o r等2 8尝试通过类型曲线分析和体积平衡方法估算给水度,他们发现将非饱和带的影响引入模型中有助于得到更为准确的给水度估计值。几年后Nw a n k w o r等2 9使用张力计和土壤湿度测量直接证明了在非承压含水层试验中,非饱和带中存在着排水过程。此后人们为了提高对非饱和带流动过程的认识开展了许多研究,并对利用参数估计方法确定土壤水力特性的可行性产生了兴趣3 0-3 2。由于非饱和带渗透系数很难通过直接测量获得,人们通常假设土壤水力特性可以用一个包含少量未知参数的相对简单的确定性模型来描述3 3。G a r d n e r3 4提出了一种简单的指数模型来描述渗透系数、含水量和毛细压力水头之间的关系,K r o s z y n s k i等3 5应用该模型得到了线性化的R i c h a r d s方程,并建立了刚性含水层的二维(2 D)轴对称解析模型,这种线性化方法已广泛用于求解非饱和带-饱和带流体流动耦合问题3 6-3 9。然而,M o e n c h指出,使用一个本构指数来表征相对渗透系数和含水量可能不适用于分析无承压含水层试验4 0。M a t h i a s等4 1提出了采用2个独立的指911h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 数来表征相对渗透系数和含水量的模型,M i s h r a等4 2提出了通过不同的指数和阈值分别表征相对渗透系数和含水量的模型。这些研究表明,采用更灵活的土壤水分特征曲线模型可以更好地表征非饱和带的水力特性,提高了非承压含水层试验的分析精度。对于包气带井回灌,L i a n g等4 3基于G a r d-n e r模型和K r o s z y n s k i等3 5的研究,将R i c h a r d s方程线性化,建立了半解析模型来研究回灌中的水力响应和包气带井的补给能力。Q i等4 4使用同样的方法将R i c h a r d s方程线性化,在模型中存在由地表入渗或蒸散发产生的地表通量的情况下,为包气带井回灌建立了半解析模型,研究了随时间变化的地表通量对包气带井回灌的影响。目前已有的研究没有详细讨论过三参数或四参数土壤水分特征曲线模型在包气带井分析中的应用与影响,尤其是存在地表通量的情况下使用不同土壤水分特征曲线模型的影响。使用不同的土壤水分特征曲线模型对所得水力响应具有怎样的影响,以及如何根据实际情况选择最合适的模型,仍有待进一步研究。笔者拟采用不同土壤水分特征曲线模型对包气带井回灌所引起的非饱和-饱和流分析的影响。首先建立包气带井回灌的概念模型,其中非饱和带流体流动用R i c h a r d s方程来描述,饱和带流体流动用地下水控制方程来描述。由于四参数土壤水分特征曲线模型可以通过设置特定参数值转换为二参数或三参数模型,使用更加灵活,因此我们建立的包气带井模型使用的是四参数土壤水分特征曲线模型。耦合后的非饱和带-饱和带流体流动方程用C OM-S O L M u l t i p h y s i c s偏微分方程求解器进行求解,所得结果用于研究不同非饱和带水力特性值下包气带井回灌引起的非饱和带-饱和带流体流动,并分析地表通量对回灌的影响,比较在使用不同土壤水分特征曲线模型时包气带井回灌的水力响应。1 包气带井回灌概念模型和数学模型 图1为无承压含水层中包气带井的横切面,非饱和带和饱和带都是均匀、各向同性的,并在横向上无限延伸。在含水层中初始地下水位是水平的,饱和带厚度为b,其底层为隔水层,非饱和带厚度为a,其顶层为地表,在地表上有由地表入渗或蒸散发所引起的地表通量I(t),I(t)是随时间变化的。包气带井部分穿透非饱和带,以恒定流量Q进行回灌,回灌井段顶部到初始地下水位的垂直距离为j,底部到初始地下水位的垂直距离为d。本研究中使用柱坐标系,r轴沿水平径向水流方向,z轴为垂直向上,坐标原点在回灌井的中轴线与初始潜水面相a为非饱和带厚度;b为饱和带厚度;d为回灌井段底部到初始地下水的垂直距离;I(t)为地表通量;Q为流量;r为沿水平径向水流方向的坐轴标;j为回灌井段顶部到初始地下水位的垂直距离图1 包气带井回灌以及地表通量所引起的非饱和-饱和流示意图F i g.1 S c h e m a t i c d i a g r a m s o f u n s a t u r a t e d-s a t u r a t e d f l o w i n d u c e d b y v a d o s e z o n e w e l l i n j e c t i o n a n d g r o u n d s u r f a c e f l u x交处,初始地下水位高程为z=0。非饱和带中流体流动是由R i c h a r d s方程描述的,由于R i c h a r d s方程具有非线性,并且2种不同流动状态(饱和和非饱和)之间存在移动界面(地下水位),求解图1所示包气带井回灌系统是很困难的。为了解决这一问题,本研究采用了K r o s z y n s k i等3 5提出的扰动展开方法来简化R i c h a r d s方程,该方法假设在回灌过程中流体流动仅稍微偏离平衡状态,然后将R i c h a r d s方程中的因变量展开为扰动级数,忽略其中的高阶项,将所余项代入方程,最终得到一个一阶线性化方程。描述包气带井回灌和地表通量引起的非饱和带流体流动的线性化控制方程为:Krk0(z)1rrrur +Kzzk0(z)uz =C0(z)ut,0za(1 a)u(r,z,0)=0(1 b)Kzk0(z)uz(r,z,t)z=a=I(t)(1 c)l i mru(r,z,t)=0(1 d)l i mr0ru(r,z,t)r=-Q(z)2 Kr(l-d)(0za)(1 e)k0(z)=k(0)C0(z)=C(0)(1 f)式中:u=hu-h0为非饱和带水头增量;h0为非饱和带初始水头;hu为非饱和带水头;Kr和Kz分别为径向(r)和垂向(z)饱和渗透系数;Q为回灌速率;(z)是B o x函数,定义为(z)=H(z-d)-H(z-j);H()是H e a v i s i d e阶跃函数;k0(z)和C0(z)分别是在初始含水量0时相对渗透系数021 第4期祁翠婷等:包气带井回灌引起的非饱和-饱和流分析k()(0k1)和土壤比持水量C()=d/d(0)的零阶近似;是非饱和带压力头;是体积含水量;I(t)是由地表入渗或蒸散发所引起的地表通量;I(t)为正值时表示流动方向是垂直向下的。式(1 c)定义了地表的边界条件,将u=hu-h0代入 式(1 c)得 到Kzk0(z)uzz=a=Kzk0(z)huzz=a-Kzk0(z)h0zz=a=I(t),其中Kzk0(z)huzz=a为地表上的总通量;Kzk0(z)h0zz=a为地表上的初始通量;I(t)可以看作是相对于初始状态下地 表 通 量 的 变 化。这 里 我 们 假 设Kzk0(z)h0zz=a=0,回灌前含水层水力条件处于稳态。由包气带井回灌引起的饱和带地下水流动控制方程为:Kr1rrrsr +Kz2sz2=Ssst,-bz0(2 a)s(r,z,0)=0(2 b)sz(r,z,t)z=-b=0(2 c)l i mrs(r,z,t)=0(2 d)l i mr0rs(r,z,t)r=0,-bz0(2 e)式中:s=hs-h0为饱和带水头增量;hs为饱和带水头;Ss为贮水率。在潜水面上,非饱和带和饱和带流体流动通过压力连续性和流量连续性耦合3 5,3 9,线性化后的潜水面上的界面条件可表示为4 3-4 4s-u=0,z=0(3 a)sz-uz=0,z=0(3 b)2 土壤水分特征曲线模型含水层中流体流动的模拟与预测通常依赖于数学模型,而这些模型需要对水力参数准确描述。对于非饱和带,关键参数包括体积含水量,相对渗透系数k,和非饱和带压力头。学者们为建立通过或来估计k的函数模型做了许多研究3 4,4 5-4 7。G a r d n e r3 4提出了k和随毛管压力头的增加而呈指数变化的土壤水分特征曲线模型:k()=e(-a)(4 a)()=r+Sye(-a)(4 b)式中:r为残留含水量,是包气带充分重力释水而又未受到蒸发、蒸腾消耗时的含水量(在此我们不考虑蒸发、蒸腾消耗的含水量),表示地下水面上方毛管作用所能支持的最大含水量4 8;s为饱和含水量;a是含水层由饱和变为非饱和时的压力头(a0);Sy=s-r为给水度;0为本构指数,反映了非饱和带的持水能力,其取值范围是0.25 m-14 9。土 壤 容 水 度 表 达 式 为C=d/d=Sye(-a)。基于G a r d n e r模型,K r o s z y n s k i等3 5得到零阶近似函数k0(z)和C0(z)分别为:k0(z)=e-z(4 c)C0(z)=Sye-z(4 d)M a t h i a s等4 1提出了一个含有3个拟合参数的土壤水分特征曲线模型(MB模型),该模型使用2个独立的本构指数来分别表征相对渗透系数和容水度:k()=ek(-a)(5 a)C()=Sycec(-a)(5 b)式中:k为相对渗透系数指数,k越大,相对渗透系数越小;c是持水指数,反映了非饱和带的持水能力,c越大,非饱和带持水能力越弱。持水能力可以理解为在毛细压力和重力共同作用下非饱和带吸持水分的能力,更强的持水能力(与较大的1/c值或较小的c值相关)意味着大部分水会被非饱和带吸附并保持“不动”状态,只有小部分水是“可动”的,可以通过重力排水。M i s h r a等4 2提出了一个含有4个拟合参数的土壤水分特征曲线模型(MN模型),该模型通过不同的阈值和本构指数来表征相对渗透系数和容水度:Se=()-rSy=ec(-a),c0(6 a)k()=ek(-k),k1,k k0(6 b)式中:Se是有效饱和度;k为相对渗透系数等于1时的压力头(k0),k并不总是等于a。在使用MN模型时零阶近似函数k0(z)和C0(z)为:k0(z)=ek(b1-z)(6 c)C0(z)=Syce-cz(6 d)式中:b1为a和k的差值,b1=a-k。以上3种土壤水分特征曲线模型具有内在联系,MN模型可以通过设置特定参数值转换为其他两种模型,具有更高的灵活性。当a=k时,MN模型简化为MB模型;当a=k以及k=c=时,MN模型简化为G a r d n e r模型。因此,G a r d n e r模型和MB模型可以看作MN模型的2个子集。因此本研究采用MN模型来线性化R i c h a r d s方程并求解包气带井回灌引起的水力响应。121h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 3 方程求解及水力响应分析3.1概念模型求解无量纲化是用数学工具研究物理问题时常用的方法,运用无量纲化能够化简微分方程,使方程易于求解,也 可 通 过 无 量 纲 化 求 得 微 分 方 程 的 近 似解5 0。在工程学和物理学中,使用无量纲变量能够将参数关联起来,一方面能够减少需要讨论的参数数量,避免进行大量实验,从而减少冗余数据,提高分析效率,另一方面能够降低参数间相关性的影响,更好地反映参数对结果的影响5 1。在工程问题的分析中,有许多无量纲参数都具有重要意义,如P e-c l e t数、R e y l e i g h数、B o u s s i n e s q数等。定义无量纲变量和参数如下:sD=2 Kra sQ,uD=2 Kra uQ,rD=ra,zD=za,bD=ba,tD=KrSsa2t,KD=KzKr,jD=ja,dD=da,D k=ka,D c=ca,=D cbDSy De-D kb1 D,Sy D=SyS,S=Ssb,=KDD k,b1 D=b1a=a-ka,ID=2 a2QI(t),=KDeD k(b1 D-zD)(7)式中:下标“D”代表无量纲项。使用无量纲项能够减少式中的参数总量,并使得余下的参数为相互独立的自由参数,避免参数相关性的影响,加深对参数影响的认识。将式(7)和式(6 c,6 d)代入式(13),得到非饱和带控制方程的无量纲形式:1rDrD(rDuDrD)+KD2uDzD2-uDzD=e(D k-D c)zDuDtD,0zD1(8 a)uD(rD,zD,0)=0(8 b)uDzD(rD,zD,tD)zD=1=ID(tD)(8 c)l i mrDuD(rD,zD,tD)=0(8 d)l i mrD0rDuD(rD,zD,tD)rD=-(zD)jD-dD,0zD1(8 e)饱和带控制方程的无量纲形式:1rDrD(rDsDrD)+KD2sDzD2=sDtD,-bDzD0(9 a)sD(rD,zD,0)=0(9 b)sDzD(rD,zD,tD)zD=-bD=0(9 c)l i mrDsD(rD,zD,tD)=0(9 d)l i mrD0rDsD(rD,zD,tD)rD=0,-bDzD0时,非饱和带的储水能力减小,使得回灌过程中uD和sD增大;当b1 D减小到0以下时,非饱和区渗透系数降低,导致uD和sD减小。以b1 D=0时的水头增量曲线为参考,当b1 D的绝对值较小时,其对uD和sD的影响也较小。综上所述,无量纲相对渗透系数指数D k对回灌过程中水头增量的影响较大,D k的变化会影响非饱和带渗透系数,随着D k的增大,非饱和带渗透系数减小的速度逐渐加快,导致uD和sD的增量逐渐减小并发生延迟。无量纲持水指数D c对包气带井回灌引起的水力响应的影响也较大,D c的变化a.取样点rD=0.1,zD=0.1,参数D k=3,D c=1;b.取样点rD=0.1,zD=-0.1,参数D k=3,D c=1图5 不同b1 D=D a-D k对应的无量纲水头增量uD(a)和sD(b)随tD的变化F i g.5 D i m e n s i o n l e s s h y d r a u l i c h e a d i n c r e m e n t s uD(a)a n d sD(b)v s.t h e d i m e n s i o n l e s s t i m e tD f o r d i f-f e r e n t v a l u e s o f b1 D=D a-D k主要影响非饱和带持水能力,随着D c的增大,非饱和带持水能力减小,导致在回灌过程中uD和sD的增量逐渐变大。将指数参数简化为D k=D c=D会给包气带井回灌引起的水力响应的计算和预测带来一定的误差。2个无量纲压力头阈值的差值b1 D=D a-D k的绝对值较小时,对uD和sD的影响也421 第4期祁翠婷等:包气带井回灌引起的非饱和-饱和流分析较小,此时将其简化为b1 D=D a-D k=0对水力响应的计算和预测带来的误差较小。3.3由地表入渗或蒸散发引起的地表通量对包气带井回灌水力响应的影响在包气带井回灌过程中,非饱和带和饱和带水头响应也受到由地表入渗或蒸散发引起的地表通量(g r o u n d s u r f a c e f l u x,简称G S F)的影响。Q i等4 4基于G a r d n e r模型,在考虑了G S F的情况下,求得了包气带井回灌的半解析解,研究了G S F和非饱和带水力参数对回灌的影响,并发现G S F对回灌的影响也受非饱和带水力参数控制。为了进一步研究G S F和非饱和带参数D k、D c及b1 D=D a-D k的影响,我们在模型中引入G S F,其值为ID,假设ID为常数,正值表示由地表入渗产生的G S F,负值表示由蒸散发产生的G S F。a.取样点rD=0.1,zD=0.1,参数D c=3,b1 D=D a-D k=0;b.取样点rD=0.1,zD=-0.1,参数D c=3,b1 D=D a-D k=0。图中虚线代表施加相同ID情况下基于G a r d n e r模型所求得的不同D对应的解图6 在tD2 0 0,2 5 0)期间施加ID=1时,不同无量纲相对渗透系数指数D k对应的时间无量纲水头增量uD(a)和sD(b)随tD的变化F i g.6 F o r d i f f e r e n t v a l u e s o f t h e d i m e n s i o n l e s s r e l a t i v e h y d r a u l i c c o n d u c t i v i t y e x p o n e n t D k:T h e d i m e n-s i o n l e s s h y d r a u l i c h e a d i n c r e m e n t s uD(a)a n d sD(b)v s.t h e d i m e n s i o n l e s s t i m e tD w i t h ID=1 d u r i n g tD2 0 0,2 5 0)假设在tD2 0 0,2 5 0)施加G S F为ID=1,在D k=0.1,1和2时(D c=3,b1 D=D a-D k=0)所求得的uD(rD=0.1,zD=0.1)和sD(rD=0.1,zD=-0.1)随时间变化的曲线如图6中实线所示。图中虚线代表在相同的ID下,基于G a r d n e r模型求得D=0.1,1和2时的水头增量随时间的变化。图6显示由地表入渗引起的G S F会使uD和sD增大,这是因为更多的水通过地表入渗进入了非饱和带,然后继续向下进入饱和带,使得回灌过程中水头增量变大。在图6中,D k增大时,非饱和带渗透系数降低,导致水通过含水层向下入渗变得更加困难,ID对非饱和带和饱和带水头变化的影响也随之减弱。图6中实线与虚线的比较表明,D k与D在ID存在的情况下对回灌中水头增量的影响是相似的,并且当D k的值接近D c时,相同颜色的实线与虚线间的差异减小。当ID=-1,由蒸散发引起的G S F对uD和sD的影响与图5相反,会使uD和sD减小,这是由于部分水分通过蒸散发作用进入空气,使得回灌过程中水头增量减小,D k增大时,非饱和带渗透系数减小,地表水分蒸发难度增大,ID对uD和sD的影响减小,此时D k与D对回灌的影响相似。a.取样点rD=0.1,zD=0.1,参数D k=0.5,b1 D=D a-D k=0;b.取样点rD=0.1,zD=-0.1,参数D k=0.5,b1 D=D a-D k=0。图中虚线代表施加相同ID情况下基于G a r d n e r模型所求得的不同D对应的解图7 在tD2 0 0,2 5 0)期间施加ID=1时,不同无量纲持水指数D c对应的无量纲水头增量uD(a)和sD(b)随tD的变化F i g.7 F o r d i f f e r e n t v a l u e s o f t h e d i m e n s i o n l e s s r e l a t i v e h y d r a u l i c c o n d u c t i v i t y e x p o