有入渗补给的层状非均质含水层圆岛潜水井流模型.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 5-2 6.C h e n C h o n g x i,T a n g Z h o n g h u a,X i e Y o n g h u a,e t a l.A w e l l f l o w m o d e l f o r a s t r a t i f i e d h e t e r o g e n o u s u n c o n f i n e d a q u i f e r i n a r o u n d i s l a n d w i t h i n f i l t r a t i o n r e c h a r g eJ.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 5-2 6.基金项目:国家自然科学基金项目(4 1 9 7 2 2 6 3)作者简介:陈崇希(1 9 3 3),男,教授,博士生导师,主要从事地下水运动解析模型与数值模拟方法的研究与教学工作。E-m a i l:c c x 3 31 6 3.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,唐仲华1,谢永桦2,王旭升2(1.中国地质大学(武汉)环境学院,武汉 4 3 0 0 7 8;2.中国地质大学(北京)水资源与环境学院,北京 1 0 0 0 8 3)摘 要:作为圆岛状均质含水层稳定井流模型,经典的D u p u i t井流模型既没有考虑普遍存在的降水入渗补给,也不适用于层状非均质含水层系统,有必要加以完善。在考虑降水入渗补给改进 D u p u i t井流模型的基础上,进一步将其拓展到层状非均质潜水含水层。引入 (吉林斯基)势函数,根据水均衡原理建立极坐标下的地下水流微分方程,再依边界条件解析得到了相应的流量方程、水位方程和分水岭公式。以双层结构为例,观察3 0组不同参数条件的典型水位曲线组,发现不同渗透系数的水位曲线交于一点的特殊现象并从理论上给出了证明。解析模型仍引入 D u p u i t 假定,且没有考虑抽水井的井壁“水跃”现象,为判断这些条件对解析公式适用能力的影响,建立轴对称剖面二维流数值模型并做了对比研究。除抽水井附近外,水位解析方程产生的相对误差一般低于4%,在最偏离 D u p u i t 假定的分水岭处,距离和水位的解析误差均小于0.1%。D u p u i t假定并没有严重影响解析模型的适用性。关键词:D u p u i t公式;降水入渗补给;层状非均质;分水岭;水跃;解析解;数值模拟2 0 2 2-1 2-3 1收稿;2 0 2 3-0 2-2 2修回;2 0 2 3-0 2-2 4接受中图分类号:P 6 4 1 文章编号:2 0 9 6-8 5 2 3(2 0 2 3)0 4-0 0 1 5-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 2 3 开放科学(资源服务)标识码(O S I D):A w e l l f l o w m o d e l f o r a s t r a t i f i e d h e t e r o g e n o u s u n c o n f i n e d a q u i f e r i n a r o u n d i s l a n d w i t h i n f i l t r a t i o n r e c h a r g eC h e n C h o n g x i1,T a n g Z h o n g h u a1,X i e Y o n g h u a2,W a n g X u s h e n g2(1.S c h o o l o f E n v i r o n m e n t S t u d i e s,C h i n a U n i v e r s i t y o f G e o s c i e n c e s(Wu h a n),Wu h a n 4 3 0 0 7 8,C h i n a;2.S c h o o l o f W a t e r R e s o u r c e s&E n v i r o n m e n t,C h i n a U n i v e r s i t y o f G e o s c i e n c e s(B e i j i n g),B e i j i n g 1 0 0 0 8 3,C h i n a)A b s t r a c t:O b j e c t i v eT h e D u p u i t m o d e l o f w e l l f l o w i s a c l a s s i c a l s t e a d y-s t a t e w e l l f l o w m o d e l f o r a h o m o-g e n e o u s u n c o n f i n e d a q u i f e r i n a r o u n d i s l a n d.H o w e v e r,i t d o e s n o t c o n s i d e r t h e w i d e l y e x i s t i n g i n f i l t r a t i o n r e c h a r g e f r o m p r e c i p i t a t i o n,a n d i s a l s o i n a p p l i c a b l e f o r s t r a t i f i e d h e t e r o g e n e o u s a q u i f e r s y s t e m s.T h e r e-f o r e,i t m u s t b e m o d i f i e d t o a d d r e s s t h e s e i s s u e s.M e t h o d sO n t h e b a s i s o f t h e r e v i s e d D u p u i t w e l l f l o w m o d e l w h i c h i n c o r p o r a t e s i n f i l t r a t i o n r e c h a r g e,t h i s s t u d y f u r t h e r e x t e n d e d i t s a p p l i c a t i o n t o a s t r a t i f i e d h e t e r o g e n e o u s u n c o n f i n e d a q u i f e r.T h e G i r i n s k i is p o t e n t i a l f u n c t i o n w a s u s e d t o c o n s t r u c t t h e d i f f e r e n t i a l e q u a t i o n s f o r t h e r a d i a l g r o u n d w a t e r f l o w a c c o r d i n g t o t h e w a t e r b a l a n c e p r i n c i p l e,a n d t h e a n a l y t i c a l s o l u-t i o n s s a t i s f y i n g t h e b o u n d a r y c o n d i t i o n s a r e t h e n o b t a i n e d a s f o r m u l a s o f t h e f l o w r a t e,w a t e r t a b l e a n d g r o u n d w a t e r d i v i d e.T a k i n g t h e b i l a y e r s t r u c t u r e a s a n e x a m p l e,t y p i c a l g r o u n d w a t e r l e v e l c u r v e s w i t h r e-s p e c t t o 3 0 s c e n a r i o s o f d i f f e r e n t p a r a m e t e r v a l u e s w e r e i n v e s t i g a t e d.A s p e c i a l p h e n o m e n o n w a s f o u n d i n w h i c h t h e c u r v e s o f d i f f e r e n t h y d r a u l i c c o n d u c t i v i t i e s i n t e r s e c t a t a s i n g l e p o i n t,w h i c h c o u l d a l s o b e p r o v e n h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 i n t h e o r y.T h i s a n a l y t i c a l m o d e l s t i l l a d o p t e d t h e D u p u i t a s s u m p t i o n a n d d i d n o t c o n s i d e r t h e h y d r a u l i c j u m p p h e n o m e n o n o n t h e w a l l o f t h e p u m p i n g w e l l.T o c h e c k t h e i m p a c t o f t h e s e c o n s t r a i n t s o n t h e a p p l i-c a b i l i t y o f t h e a n a l y t i c a l f o r m u l a s,a t w o-d i m e n s i o n a l n u m e r i c a l m o d e l f o r a x i a l l y s y mm e t r i c s e e p a g e w a s b u i l t f o r c o m p a r i s o n.R e s u l t sA s i n d i c a t e d b y t h e r e s u l t s,t h e r e l a t i v e e r r o r o f t h e g r o u n d w a t e r l e v e l e s t i-m a t e d f r o m t h e a n a l y t i c a l s o l u t i o n i s g e n e r a l l y l e s s t h a n 4%,e x c e p t f o r t h e z o n e n e a r t h e p u m p i n g w e l l.O n t h e g r o u n d w a t e r d i v i d e,w h e r e t h e D u p u i t a s s u m p t i o n i s m o s t l y i n v a l i d,t h e r e l a t i v e e r r o r s o f t h e a n a-l y t i c a l s o l u t i o n t o b o t h t h e d i s t a n c e a n d h e i g h t o f t h e d i v i d e a r e s m a l l e r t h a n 0.1%.C o n c l u s i o nT h e D u-p u i t a s s u m p t i o n d o e s n o t s i g n i f i c a n t l y i n f l u e n c e a p p l i c a b i l i t y o f t h e a n a l y t i c a l m o d e l.K e y w o r d s:D u p u i t e q u a t i o n;p r e c i p i t a t i o n i n f i l t r a t i o n r e c h a r g e;s t r a t i f i e d h e t e r o g e n e o u s;g r o u n d w a t e r d i-v i d e;h y d r a u l i c j u m p;a n a l y t i c a l s o l u t i o n;n u m e r i c a l m o d e l l i n gR e c e i v e d:2 0 2 2-1 2-3 1;R e v i s e d:2 0 2 3-0 2-2 2;A c c e p t e d:2 0 2 3-0 2-2 4 法国科学家D u p u i t1在1 8 6 3年最初提出的潜水井流模型,即D u p u i t 模型,是“圆岛状含水层潜水稳定井流模型(外边界为圆柱形定水位条件)”,而不是某些文献(如卡明斯基2,贝尔3等)所述的“无界含水层中影响半径稳定井流模型”。也就是说,D u-p u i t 公式中的R是圆岛半径、而非影响半径。无界含水层中的井流是不可能形成稳定状态的4。关于恢复D u p u i t模型为“圆岛模型”的历史过程,参见文献5。上述D u p u i t模型只有侧向的圆岛外边界(湖海)及抽水井的井壁内边界,而不涉及上边界降水入渗补给。如此一旦抽水,经过一定时间之后必导致湖海水入侵。然而自然条件下大气降水的入渗补给是普遍现象。这样一来,D u p u i t 模型基本上只能在旱季用以地下水井流试验求取含水系统的参数,而不能够用以预测。海岛地下水资源开发利用需要的是合理截取大气降水入渗补给的理论,然而奉为经典的 D u p u i t 模型却没有这个功能,有些遗憾!2 0 2 0年,陈崇希6提出改进 D u p u i t 稳定潜水井流模型,考虑了地面均匀稳定入渗补给(蒸发排泄视为其负值)作用,可用来预测相应条件的地下水开釆效果和地下分水岭的偏移。除均匀稳定入渗之外,改进模型的其他条件与D u p u i t模型相同。该模型提出后,笔者重新对国内外文献进行了挖掘式的搜索,试图检查此前文献调研是否遗漏了同类研究。最终发现,国外学者对D u p u i t模型也做了类似的改进:前苏联学者吉林斯基7的专著把入渗补给加入到D u p u i t模型,推导出了通用流量公式,但没有对分水岭进行研究;前苏联学者阿拉文等8的专著只给出了均匀入渗条件的水位公式,没有提供流量和分水岭公式;美国学者T o d d9写的教材推导了均匀入渗条件的水位和流量公式,但没有研究分水岭;美国学者H a i t j e m a1 0在其专著中提供了均匀入渗条件的水位和分水岭公式;挪威学者K i t t e r d1 1也对均 匀 入 渗 条 件 进 行 了 解 析 解 的 推 导。只 有K i t t e r d1 1引用了H a i t j e m a1 0的研究,其他学者都没有引用前人文献,说明考虑入渗补给改进D u p u i t公式在学术界一直没有得到重视,只是被少数学者独立进行了研究。在D u p u i t模型中增加降水入渗做出改进,提出更加具有普遍意义的潜水井流水位、流量公式和分水岭公式,这项工作对完善地下水动力学的基础理论很有价值。不过,现有的改进模型仍然属于均质含水层模型,不适应自然界含水系统多数具有层状非均质特征的情况。本研究的目标,是在此基础上再次做出改进,将其拓展到水平层状非均质含水系统。笔者拟提出一个新的解析模型,并通过数值模拟对照,分析讨论D u p u i t假定的引入对解析方程适用性的约束等问题。1 模型与控制方程1.1概念模型与假设条件隔水底板水平的层状非均质圆岛状潜水层,上边界具稳定均匀入渗补给,外边界是定水位条件,圆岛中心有一口定流量或定水位的抽水井。此条件下,地下水可形成径向(轴对称)稳定流。假定渗流服从D a r c y定律并满足D u p u i t假定。以潜水层底面为基准面,取柱坐标系统(图1),且以井心隔水底面处为原点;由于地下水井流中的流量通常定义抽水流量为正值,故记流速向井(原点)为正。1.2吉林斯基势函数方程对于水平层状非均质含水层,引入 势函数(见:阿拉文等8;陈崇希1 2;贝尔3):G=ni=1KiMi(h-zi)(1)式中:Ki和Mi分别为第i层介质的渗透系数与饱和带厚度(图2);zi为第i层饱和带中点的高程。如果第i层的顶部位于潜水面以下,其饱和带厚度就是实际的介质层厚度。61 第4期陈崇希等:有入渗补给的层状非均质含水层圆岛潜水井流模型K0,K1,K2,Kn分别为第1,第2,第n介质层渗透系数(m/d);Qw为抽水流量(m3/d);hR为圆岛外边界水位(m);hw为开采井中水位(m);R为圆岛半径(m);a为无分水岭水位曲线;b为有分水岭水位曲线图1 层状非均质含水层具入渗补给的圆岛状潜水井流模型图F i g.1 M o d e l m a p o f w e l l f l o w i n a r o u n d i s l a n d w i t h a s t r a t i f i e d h e t e r o g e n e o u s u n c o n f i n e d a q u i f e r o b t a i-n i n g i n f i l t r a t i o n r e c h a r g ez1,z2,zn分别为抽水井所揭露第1,第2,第n介质层底部标高(m);M1,M2,Mn分别为第1,第2,第n介质层的水带厚度(m)图2 层状非均质含水层吉林斯基势函数定义要素图F i g.2 E l e m e n t s o f G i r i n s k i is p o t e n t i a l f u n c t i o n f o r a s t r a t i f i e d h e t e r o g e n e o u s a q u i f e r如此,对于层状非均质含水层的潜水井流问题,用 势函数表示任意断面的流量为:Q(r)=2 rdGdr(2)2 解析解2.1建立水位方程和流量方程依据质量守恒原理,并考虑抽水井流量在地下水动力学中通常定义为正值,故此流速向井(原点)为正值,因此任意断面r的流量Q也可以写为:Q=Qw-(r2-r2w)(3)式中:Qw为抽水井流量;rw为抽水井的半径;为入渗强度。式(3)引入式(2)的 势函数G并积分,得任意断面的势函数方程:G=GR-Qw2+2r2w l nRr+4(R2-r2)(4)其中:GR=n-1i=1KiMihR-zi+12KnhR-n-1j=1Mj 2(5)式中:R为圆岛半径;GR为外边界(湖海)的吉林斯基势函数。对式(4),取r=rw,记对应的G为Gw,得抽水井中水位势函数方程:Gw=GR-Qw2+2r2w l nRrw+4(R2-r2)(6)式中:Gw为抽水井中水位的吉林斯基势函数。由式(4),可计算任一点r的G,进而可获得h,由此得潜水位曲线。这里计算涉及G与h间的转化。若已知h,则可利用G函数的定义式(1)直接计算G值。若已知G,则可按下列两步骤计算其h值:首先确定地下水位位于哪一层?记水位位于第n层。先计算水位位于各层顶面水位h(i)对应的G(i)。若 有G(i)0,则外边界补给地下水(水位线如图1的a曲线)。若QR=0,则外边界与地下水互相没有补排。若QR0,则外边界排泄地下水,这时存在地下分水岭(水位线如图1中的b曲线)。2.2地下分水岭地下水分水岭位置rd满足方程:dhdr=dhdGdGdr=0(1 5)只要dG/dr=0即有dh/dr=0。另外,由式(4),有:dGdr=Qw2 r-r 2-2rr2w(1 6)令dh/dr=0,即dG/dr=0,则r=rd,根据式(1 6)得地下水分水岭位置:r2d=Qw-r2w(1 7)其物理意义十分明确:抽水井的流量等于分水点与井壁之间圆环范围内单位时间获得的入渗补给量。当抽水井半径相对很小,式(1 7)等号右边可近似为Qw/。建立地下水分水点处的水位方程:在式(4)中取r=rd,记对应的G为Gd,得:Gd=GR-Qw2 l nRrd+4(R2-r2d)-2r2wl nRrd(1 8)由此结合 势函数的定义(1),可以得到分水岭的水位hd。2.3未抽水时的天然水位分布已知无入渗且未抽水时圆岛中心的势函数为:Gr=0(=0)=GR(1 9)有入渗、未抽水时,则圆岛中心的势函数为:Gr=0(=0)=GR+4R2(2 0)此即未抽水条件下,位于圆心处的势函数,属于最大值Gm a x=Gr=0。从式(7)和(2 0)可得有入渗未抽水条件下r=0处的地下水位峰值:hp=-n-1i=1(Ki-Kn)Mi/Kn+(n-1i=1(Ki-Kn)Mi)2-2Kn12Kn(n-1i=1Mi)2-n-1i=1KiMi(i-1j=1Mj+Mi2)-GR-4R2Kn(2 1)对于双层结构,上式可以简写为:hp=-(K1-K2)M1+(K1-K2)(K2-K1)M21+2K2(G+GR+4R2)K2(2 2)其中:GR=K1M1hR-M12 +K2(hR-M1)2/2(2 3)3 典型条件下水位分布特点3.1双层结构含水层天然状态水位线为了对层状非均质含水层潜水井流的水位线分布有个具体的概念,这里对两层结构含水层潜水井流水位线给出3 0组典型曲线组(图38)。其中采取的参数:R=2 0 0 0 m;hR=5 0 m;K1=1 0 m/d;rw=0.1 m;M1=3 0 m;M2取值使含水层具足够厚度,以保证潜水位不会溢出地面。不同情景下的K2和取值见表1,天然状态下Qw=0 m3/d,抽水状态下Qw=5 0 0 0 m3/d。得到3 0组情景的特征点的信息,包括井中水位hw,潜水位分水岭的位置rd及分水岭的水位hd,以及不同值的水位交叉点的位置rc及其水位hc。当Qw为零时,无分水岭,但在r=0处存在地下水位峰值(表13),其数值随入渗强度的增大和K2值的减小而升高。表1 Qw=0 m3/d时不同K2值r=0处水位峰值hpT a b l e 1 M a x i m u m g r o u n d w a t e r l e v e l,hp a t r=0 f o r d i f f e r-e n t K2 v a l u e s w h e n Qw=0 m3/d/(md-1)水位峰值hp/mK2=1K2=5K2=1 0K2=2 0K2=1 0 00.0 0 0 7 55 2.3 45 1.8 55 1.4 85 1.0 65 0.3 20.0 0 15 3.1 15 2.4 65 1.9 65 1.4 05 0.4 30.0 0 1 55 4.6 55 3.6 75 2.9 25 2.0 85 0.6 4 注:K2单位为m/d3.2双层结构含水层抽水状态水位线抽水状态的水位线表示在图4中,水位hw表示在表2中。这些特征值随有关参数的变化是符合基本概念的。在给定的条件下(抽水井流量Qw不变),抽水井井中水位hw随K2值的减小而降低,随81 第4期陈崇希等:有入渗补给的层状非均质含水层圆岛潜水井流模型图3 双层结构含水层天然状态水位曲线计算案例(.入渗强度;K2.水平渗透系数;下同)F i g.3 E s t i m a t e d r e s u l t s o f g r o u n d w a t e r l e v e l c u r v e s i n t h e n a t u r a l s t a t e f o r a n a q u i f e r o f t h e b i l a y e r s t r u c t u r e图4 双层结构含水层在Qw=5 0 0 0 m3/d抽水状态水位曲线计算案例F i g.4 E s t i m a t e d r e s u l t s o f g r o u n d w a t e r l e v e l c u r v e s i n t h e p u m p i n g s t a t e o f Qw=5 0 0 0 m3/d f o r a n a q u i f e r o f t h e b i-l a y e r s t r u c t u r e表2 Qw=5 0 0 0 m3/d时不同K2值抽水井的井中水位hwT a b l e 2 W a t e r l e v e l w i t h i n t h e p u m p i n g w e l l f o r d i f f e r e n t K2 v a l u e s w h e n Qw=5 0 0 0 m3/d/(md-1)井中水位hw/mK2=1K2=5K2=1 0K2=2 0K2=1 0 00.0 0 0 7 52 6.7 22 9.5 63 2.7 73 7.6 34 6.6 60.0 0 12 7.6 43 0.4 03 3.5 23 8.1 74 6.7 80.0 0 1 52 9.3 93 2.0 33 4.9 83 9.2 34 7.0 3 注:K2单位为m/d入渗强度的增大而升高。在同一入渗强度条件下,地下水位分水岭的位置rd不变,分水岭水位hd随K2值的减小而升高。显然,地下水位分水岭的位置rd随入渗强度的增大(抽水井流量不变)而缩小,分水岭的水位hd随入渗强度的增大而升高。很有意思,所有不同K2值水位曲线都交于同一个点。这一现象,以前在均质含水层中未曾出现过,值得探究。交叉点的位置rc随入渗强度的增大而缩小,但交叉点的水位hc与哪些因素有关还需要进一步证明。3.3水位线交叉点可以预计,当抽水流量小于区内单位时间的总入渗补给量的条件下,会存在地下水位分水岭。也可预计,2条不同K2值水位线存在交叉现象。然而,如图4所示5条不同K2值的水位线竟会交于一点(表3)。这具有重要意义,我们将给予重点分析。关于交叉点的位置,可以针对以下2个问题开展分析。表3 Qw=5 0 0 0 m3/d时分水岭和交叉点特征T a b l e 3 C h a r a c t e r i s t i c s o f g r o u n d w a t e r d i v i d e a n d i n t e r s e c t i o n p o i n t f o r Qw=5 0 0 0 m3/d/(md-1)分水岭位置rd/m不同K2值(m/d)对应分水岭处水位hd(m)K2=1K2=5K2=1 0K2=2 0K2=1 0 0交叉点位置rc/m交叉点水位hc/m0.0 0 0 7 51 4 5 6.7 3 25 0.3 15 0.2 55 0.2 05 0.1 45 0.0 49 7 4.9 55 0.0 0 00.0 0 11 2 6 1.5 6 75 0.7 35 0.5 95 0.4 75 0.3 35 0.1 06 5 0.0 25 0.0 0 00.0 0 1 51 0 3 0.0 6 55 1.7 95 1.4 25 1.1 45 0.8 15 0.2 53 1 8.5 55 0.0 0 0 注:K2单位为m/d91h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 (1)第1层是否可能有交点?若第1层(n=1)存在交点,并记交点的位置为(rc,hc),依式(4)和式(1),有:12K1h2=GR-Qw2 l nRr+4(R2-r2)-r2w2l nRr(2 4)当K2分别取K2-1和K2-2时,应满足:K1h2c2=GR(K2=K2-1)-Qw2 l nRrc+4(R2-r2c)-r2w2l nRrc(2 5)K1h2c2=GR(K2=K2-2)-Qw2 l nRrc+4(R2-r2c)-r2w2l nRrc(2 6)已知GR(K2=K2-1)=K1M1(hR-M12)+K2-1(hR-M1)2/2(2 7)GR(K2=K2-2)=K1M1(hR-M12)+K2-2(hR-M1)2/2(2 8)于是要求:K2-1(hR-M1)2/2=K2-2(hR-M1)2/2(2 9)由于K2-1肀K2-2故此式不成立,即在此条件下,在第一层中不可能存在交点。(2)第2层是否可能有交点?若第2层(n=2)存在交点,并记交点的位置为(rc,hc),依式(4)和式(1),有:K1M1h-M12 +12K2(h-M1)2=GR-Qw2 l nRr+4(R2-r2)-r2w2l nRr(3 0)当K2分别取K2-1和K2-2时,应满足:K1M1hc-M12 +12K2-1(hc-M1)2=GR(K2=K2-1)-Qw2 l nRrc+4(R2-r2c)-r2w2l nRrc(3 1)K1M1hc-M12 +12K2-2(hc-M1)2=GR(K2=K2-2)-Qw2 l nRrc+4R2-r2c -r2w2l nRrc(3 2)把式(2 7,2 8)分别代入式(3 1,3 2),得到:(K2-1-K2-2)(hc-M1)2/2=(K2-1-K2-2)(hR-M1)2/2(3 3)于是有:hc=hR(3 4)由此证明了2条曲线的交点的水位与外边界水位等高。由于K2-1和K2-2的取值是任意的,所以任何不同K2值的曲线相交于同一点,且交点处水位与外边界水位等高。此结论可推广到所有其他因素相同条件下,不同K2值的曲线都交于同一点,并等于外边界水位。明确了hc=hR之后,就可以直接根据势函数的解析式(4)来确定rc的值,因为交叉点的势函数与外周边界势函数相等,即:GR-Qw2+2r2w l nRrc+4(R2-r2c)=GR(3 5)等价于以下方程:4(R2-r2c)-Qw2+2r2w l nRrc=0(3 6)可以进一步改写为:(rc/R)2-12 l n(rc/R)=Qw+r2w R2=r2d+r2w R2=r2d+r2wR2(3 7)为了方便求出rc,定义以下关于r/R的函数:f(rR)=(r/R)2-12 l n(r/R)(3 8)该函数的曲线形态如图5所示。只要(Qw+r2w)(R2),即潜水面具有分水岭,必定在某个r=rc处满足f=(Qw+r2w)/(R2)。根据f(r/R)曲线的单 调 特 征,方 程f(r/R)=(Qw+r2w)/(R2)只有一个根rc满足rc1 0 0rw)时,式(3 6)可近似为:02 第4期陈崇希等:有入渗补给的层状非均质含水层圆岛潜水井流模型f(rR)=(rdR)2(3 9)4 数值模拟检验D u p u i t假定对解析解的影响4.1问题与目标在建立解析方程的过程中,假设满足D u p u i t 假定。引入D u p u i t假定,对本问题解析研究可以降维(略去z变量)而使有关方程变得简单明了。然而实际上,在某些r断面上可能出现见明显偏离D u-p u i t 假定的,例如在地下水分水岭附近,其解析解的水位是否会明显偏离实际情况,需要研究。另外由于解析解模型没有考虑水跃(出渗面)的存在 这也是由于D u p u i t 假定所引起。在抽水井附近,解析解的水位会明显偏离实际情况,它有多大的程度,多大的范围。为了具体了解D u p u i t 假定对解析结果带来什么样的误差,对轴对称剖面二维流做了数值模拟研究。4.2数值模型设计地下水潜水井流解析模型,包括D u p u i t模型和T h e i s模型1 3在内,通常都是已知抽水井流量,来计算含水层中的水头和水位分布,包括井中的水位在内,而且认为不同深处的井中水头是相同的hw。但是按当前数值模拟普遍的做法就有困难了。因为人们不知道抽水井的流量在井中如何分配,能使得井中的水头保持同一个数值hw。这次的模拟,为了与解析模型对比,我们采取的是把解析模型井流的流量边界改为定水头hw边界,在其上面加了一段出渗面,后者的水头等于位置高程z。采取这个方法可以与解析模型相对比。井中水位hw是根据上述解析方程得出来的。然后我们模拟井壁出流的总流量(水位之下的与出渗面出流的之和),与解析法抽水井的流量作对比。下面所做的模型,按这种方法进行。所研究模型,除了未引入D u p u i t 假定及抽水井井壁和外边界(在某些条件下)具出渗面之外,其他都与上述解析模型一致,所取参数与上述第3节算例一样。即,模型如图6所示:A B段为抽水井井壁出渗面,为水头已知边界,其位置高程z(i)即为其水头值;B C段为抽水井中定水头边界,其水头值等于井中水位hw;C D段为隔水底面;D E段为圆岛外边界,其水头为hR;E F为圆岛外边界水位之上壁面的出渗面,为水头已知边界,其位置高程z(i)即为其水头值;F G曲线段为具均匀稳定入渗补给的潜水面。b-b 水平线为1和2含水层的层界面。数值模型的参数与3.1节模型的参数一致(表1)。数 值 模 拟 采 用 的 是 有 限 差 分 法(MO D F-图6 双层结构含水层轴对称剖面二维流概念模型示意图F i g.6 S c h e m a t i c d i a g r a m o f t w o-d i m e n s i o n a l f l o w c o n-c e p t u a l m o d e l i n a x i s y mm e t r i c s e c t i o n o f b i l a y e r s t r u c t u r e a q u i f e rL OW1 4),建立直角坐标系的等效剖面二维流模型,用 于 模 拟 柱 坐 标 系 的 轴 对 称 剖 面 二 维 流 现象1 5。对于本研究的模型,以下模拟要素做了特殊的处理。(1)网格设计模型垂直z方向上,考虑到模型为层状非均质含水系统,具K值突变层界面,且边界条件具出渗面与井水位之下定水头条件的差别,故模型采取1 m上下。模型水平x方向上,总体均匀分布,考虑到抽水井和地下水分水岭附近流网的复杂性,这两处距适当加密。如此,模型垂直z方向上,最多共分为5 4层(由于不同K2其hw和出渗面高程不同而分层数不同),水平x方向上最多划分1 7 6个网格。(2)渗出面边界MO D F L OW1 4提供了D r a i n