二维承压非稳定流水均衡区间的数值模拟.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):7 5-8 2.D o n g G u i m i n g,W a n g Y i n g,Z h a n H o n g b i n,e t a l.N u m e r i c a l s i m u l a t i o n o f t h e w a t e r b u d g e t i n t e r v a l f o r u n s t e a d y t w o-d i m e n s i o n-a l c o n f i n e d f l o wJ.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):7 5-8 2.基金项目:国家自然科学基金项目(4 1 2 0 2 1 7 9);江苏省生态环境保护地下水监测监控与污染控制重点实验室开放课题(GWK L 2 2 0 1)作者简介:董贵明(1 9 7 9),男,副教授,主要从事地下水数值模拟研究工作。E-m a i l:g u i m i n g 1 4 4 3 21 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,王 颖2,詹红兵3,田 娟4,李嘉宁1,代丽娜1(1.中国矿业大学资源与地球科学学院,江苏 徐州 2 2 1 1 1 6;2.四川省水资源调度管理中心,成都 6 1 0 0 1 7;3.美国德克萨斯农工大学地质与地球物理系,德克萨斯 7 7 8 4 3-3 1 1 5;4.江苏师范大学 地理测绘与城乡规划学院,江苏 徐州 2 2 1 1 1 6)摘 要:由于水文地质条件的复杂性和建模工作投入的有限性,地下水数值模型往往存在不确定性。近5 0年来,随机方法是其不确定性分析的主要方法之一。区间不确定性与随机不确定性不同,是将水文地质参数等看做区间(范围),而不再考虑其随机特征。从区间不确定性角度出发,以非稳定二维承压水流为例,提出了在已知水文地质参数等为区间的情况下,基于一阶摄动展开的地下水均衡项区间的数值模拟方法。基于地下水流和污染物迁移三维数值模拟程序G FM o d e l实现了这种方法。算例分析表明,本研究方法当参数变化率在0.1以内的时候,偏差相对误差整体上可以控制在1 0%以内,该方法的计算效率明显高于等间距连续采样法。该方法的结果中不包含随机统计信息,但在已知水文地质参数等为区间的情况下,计算出水均衡项的区间,将能在地下水资源的利用和保护决策中提供一定的理论依据。关键词:地下水数值模拟;水均衡;区间;G FM o d e l;不确定性;非稳定流2 0 2 3-0 1-1 7收稿;2 0 2 3-0 3-3 1修回;2 0 2 3-0 4-0 7接受中图分类号:P 6 4 1 文章编号:2 0 9 6-8 5 2 3(2 0 2 3)0 4-0 0 7 5-0 8d o i:1 0.1 9 5 0 9/j.c n k i.d z k q.t b 2 0 2 3 0 0 2 8 开放科学(资源服务)标识码(O S I D):N u m e r i c a l s i m u l a t i o n o f t h e w a t e r b u d g e t i n t e r v a l f o r u n s t e a d y t w o-d i m e n s i o n a l c o n f i n e d f l o wD o n g G u i m i n g1,W a n g Y i n g2,Z h a n H o n g b i n3,T i a n J u a n4,L i J i a n i n g1,D a i L i n a1(1.S c h o o l o f R e s o u r c e s a n d G e o s c i e n c e s,C h i n a U n i v e r s i t y o f M i n i n g a n d T e c h n o l o g y,X u z h o u J i a n g s u2 2 1 1 1 6,C h i n a;2.S i c h u a n P r o v i n c i a l W a t e r R e s o u r c e s D e p t.,C h e n g d u 6 1 0 0 1 7,C h i n a;3.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;4.D e p a r t m e n t o f E n v i r o n m e n t a l S c i e n c e,J i a n g s u N o r m a l U n i v e r s i t y,X u z h o u J i a n g s u 2 2 1 1 1 6,C h i n a)A b s t r a c t:O b j e c t i v eG r o u n d w a t e r n u m e r i c a l m o d e l s o f t e n h a v e u n c e r t a i n t i e s d u e t o t h e c o m p l e x i t y o f t h e h y d r o g e o l o g i c a l c o n d i t i o n s a n d t h e e c o n o m i c a n d t i m e c o n s t r a i n t s i n c o l l e c t i n g a s u f f i c i e n t l y l a r g e d a t a s e t a s i n p u t s f o r c o n d u c t i n g m o d e l l i n g e x e r c i s e s.I n t h e p a s t 5 0 y e a r s,s t o c h a s t i c m e t h o d s h a v e b e e n o n e o f t h e m a i n m e t h o d s o f u n c e r t a i n t y a n a l y s i s.T h e i n t e r v a l u n c e r t a i n t y i s d i f f e r e n t f r o m t h e s t o c h a s t i c u n c e r t a i n t y,a n d i t c o n s i d e r s t h e h y d r o g e o l o g i c a l p a r a m e t e r s a s t h e i n t e r v a l s(r a n g e s)w i t h o u t c o n s i d e r i n g t h e i r s t o c h a s-t i c p r o p e r t i e s.M e t h o d sF r o m t h e p e r s p e c t i v e o f i n t e r v a l u n c e r t a i n t y,a n u m e r i c a l s i m u l a t i o n m e t h o d b a s e d o n f i r s t-o r d e r p e r t u r b a t i o n e x p a n s i o n w a s p r o p o s e d f o r s i m u l a t i n g u n s t e a d y t w o-d i m e n s i o n a l c o n-f i n e d f l o w w i t h k n o w n h y d r o g e o l o g i c a l p a r a m e t e r s a s i n t e r v a l s i n t h i s p a p e r.T h e p r o p o s e d m e t h o d i s i m-p l e m e n t e d b a s e d o n G FM o d e l,a t h r e e-d i m e n s i o n a l(3 D)n u m e r i c a l s i m u l a t i o n p l a t f o r m f o r g r o u n d w a t e r h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 f l o w a n d p o l l u t a n t m i g r a t i o n.R e s u l t sT h e a n a l y s i s s h o w s t h a t t h e r e l a t i v e e r r o r c a n b e c o n t r o l l e d w i t h i n 1 0%w h e n t h e p a r a m e t e r c h a n g e r a t e i s l e s s t h a n 0.1.T h e c o m p u t a t i o n a l e f f i c i e n c y o f t h e p r o p o s e d m e t h-o d i s o b v i o u s l y h i g h e r t h a n t h a t o f t h e c o n t i n u o u s s a m p l i n g m e t h o d w i t h e q u a l s p a c i n g.C o n c l u s i o nT h i s m e t h o d a l l o w s t h e i n t e r v a l o f t h e h e a d o r w a t e r b u d g e t t o b e c a l c u l a t e d w i t h o u t t h e r e q u i r e m e n t f o r d e-t a i l e d s t a t i s t i c a l i n f o r m a t i o n(w h i c h i s u s u a l l y u n a v a i l a b l e i n a d v a n c e)i f t h e i n t e r v a l s o f h y d r o g e o l o g i c a l p a r a m e t e r s a r e k n o w n.I t p r o v i d e s a t h e o r e t i c a l b a s i s f o r d e c i s i o n s o n t h e u s e a n d p r o t e c t i o n o f g r o u n d w a t e r r e s o u r c e s.K e y w o r d s:g r o u n d w a t e r n u m e r i c a l s i m u l a t i o n;w a t e r b a l a n c e;i n t e r v a l;G FM o d e l;u n c e r t a i n t y;u n s t e a d y f l o wR e c e i v e d:2 0 2 3-0 1-1 7;R e v i s e d:2 0 2 3-0 3-3 1;A c c e p t e d:2 0 2 3-0 4-0 7 地下水数值模拟是地下水科学研究的重要方法。然而,由于水文地质条件的复杂性和建模工作投入的有限性,地下水数值模型往往存在不确定性,如何定量刻画模型的不确定性、提高模型的可靠性,是地下水科学的前沿问题之一1-4。近5 0年来,随机和统计方法是处理地下水不确定性 问 题 的 主 要 方 法。这 类 方 法 包 括 蒙 特 卡 洛(MC)方法、N e u m a n n展开法、T a y l o r级数展开法、摄动法、卡 尔 曼 滤 波 法、马 尔 可 夫 链蒙 特 卡 洛(MCMC)方法、通用似然不确定性估计(G L U E)方法、差异演化自适应框架(D R E AM)算法等及其各种改进优化的方法,并且这些方法有丰富的研究成果,在地下水数值模拟的相关理论和实践中发挥了重要作用5-1 8。1 9 6 6年M o o r e1 9提出了区间分析理论,为区间不确定性分析奠定了基础。目前区间分析已经应用到了许多学科的不确定性分析中。X i e等2 0将基于区间理论方法的区间固定混合随机规划(I F S P)模型应用于对不确定下的温室气体减排的管理。G a o等2 1结合了泰勒展开、矩阵摄动理论以及随机区间矩方法,建立了随机区间结构响应的均值和标准差的表达式。M i a o等2 2通过将区间参数与模糊集的概念融入D e n o v o编程框架,得出了一种不确定性下的区间模糊D e n o v o编程(I F D N P)方法,用于规划水资源管理系统。L 等2 3、Q i u等2 4采用区间分析方法开展了一系列结构力学数值模拟不确定性研究。在地下水科学的研究中,D o n g等2 5-2 7利用区间参数摄动法,开展了承压水地下水头的区间数值模拟研究。不确定性信息的描述方式取决于已知信息的数量和类型。在概率方法中,所有不确定参数的联合概率密度函数都应该是已知的,但是,当没有足够的数据来验证这些随机变量(或相关函数)的概率密度的正确性时,概率方法可能难以可靠地得到满足精度要求的计算结果2 8。为了反映客观实际,减少人为因素的影响,在信息不够充分的条件下,区间不确定性提供了一种处理不确定性问题的思路,这种方法是把这些不确定性参数视为在某一区间(范围)内变化的变量,而不考虑在区间内的概率分布情况,计算结果为节点水位、浓度或者均衡项的区间(范围),结果中也不包含任何概率和统计信息。这种描述更符合工程习惯也更容易为工程人员所接受。水位、浓度或者均衡项的区间(范围)可以评估不确定性的程度,为管理决策提供依据。区间不确定性和随机不确定性不是对立的,2种不确定性可以同时存在,区间不确定性也可以单独存在。区间不确定性本质上就是将参数看作是一个区间(严格意义上应该是闭区间,既可以取到端点值,这里的区间和数学中的区间概念是一样的),即参数是在这个区间内变化,但不关心参数在这个区间内是如何变化的,而关心当参数在这个区间内变化的时候,地下水位、浓度、水均衡值将在什么样的区间内变化。随机不确定性关心的则是如何变化的问题,即符合什么分布。如果能够知道参数的区间,那么就可以进行区间不确定性分析,如果同时能确定参数服从哪种分布,那么同时也可以开展随机不确定性分析;当没有办法知道某个参数服从什么分布,但如果能够知道该参数的区间,那么就只能进行区间不确定性分析了。本研究针对二维承压水非稳定流,从区间不确定性角度出发,开展在水文地质参数等为区间的情况下,地下水均衡项区间计算的数值模拟研究。1 G FM o d e lG FM o d e l(g r o u n d w a t e r f l o w m o d e l)是基于三角形剖分和隐式有限差分的三维地下水流及多组分污染物迁移模拟程序,由中国矿业大学董贵明开发,不包括任何其他程序的任何代码2 9。G FM o d e l按照模块化设计,代码3 0 0 0 0余行,由7 1个f 9 0文件组成。G FM o d e l的主要模块及其关系如图1所示。由图1可以看出,G FM o d e l主要由公共模块67 第4期董贵明等:二维承压非稳定流水均衡区间的数值模拟组、渗流模块组、水均衡模块组、流速模块组、污染物模块组和污染物均衡模块组6个模块组组成。模块组中包含若干个模块,其中公共模块组是整个程序中反复用到的一些数据和函数的集合。为了保持渗流计算和污染物迁移计算过程及代码的独立性,G F M o d e l中采用了渗流计算的结果输出到文件、污染物迁移计算读取文件的策略(图1中的虚线),不需要在每次进行污染物计算的时候都先进行渗流计算。图1 G FM o d e l主要模块及其关系图F i g.1 M a i n m o d u l e s o f G FM o d e l a n d t h e i r r e l a t i o n s h i p d i a g r a m G FM o d e l的各种水文地质参数、三角形活动性、源汇和节点类型及值都允许随时间变化。可以进行节点干湿变换,即含水层中的节点由于水位下降可以变成干节点,由于水位上升干节点又可以变成正常的节点。P h r e e q c接口模块可以实现G F-M o d e l调用P h r e e q c,以进行污染物的多组分迁移-反应模拟。为了解决污染物计算中存在的数值弥散和振荡问题,G FM o d e l中可以根据污染物迁移的大概范围定义其模拟范围,即污染物的模拟范围可以小于(有时是明显小于)渗流的模拟范围,这样可以在污染物迁移模拟范围内采用面积较小的三角形,这种处理方法不仅可以降低数值弥散和振荡,而且可以明显降低污染物迁移计算时间。另外,G F-M o d e l的系数矩阵和列向量是通过三角形单元渗透矩阵和单元列向量组装得到的。本研究提出的区间数值模拟方法是在G FM o d-e l的基础上实现的。2 地下水均衡项的区间数值模拟方法 地下水的区间数值模拟是当数值模型的各种输入,比如水文地质参数、源汇项和边界条件等为区间(范围)的条件下,计算地下水位、水均衡项、浓度的区间(范围)的一种数值模拟方法。或者说区间数值模拟中出现的不全部是确定的数,而是包含了一些区间数。其产生的背景是数值模拟存在不确定性。区间数值模拟是从区间不确定性角度来处理和分析不确定性问题的。2.1问题定义水文地质参数向量可以记为:=(1,2,n)(1)式中:i为第i个水文地质参数,i=1,2,n;n为水文地质参数的个数。在区间不确定性分析中,ii,i,i为区间的下限(最小值),i为区间的上限(最大值)。也可以有下面表达式:i=i0i(2)式中:i0为i的平均值或称为区间中点,且i0=(i+i)/2。i记为i的偏差,i=(i-i)/2。i0即是确定性数值模型使用的参数值。水均衡区间和确定性模型下的水均衡值的关系记作:Q()=Q0+Q()(3)式中:Q()为水文地质参数为区间时对应的水均衡77h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 区间;Q0为水文地质参数取平均值0时的水均衡值,也就是确定性数值模型下的水均衡值;Q()为水文地质参数发生偏差后引起的水均衡变化量。地下水均衡项的区间计算就是当水文地质参数在某一区间内变化时,求Q()区间的问题,也就是获得Q()的极大值和极小值的问题。另外,从式(2)中可以看出,这一问题可以转化为计算Q()的极大值和极小值的问题。地下水均衡项一般包括侧向补给量、降水入渗补给量、重力或弹性释水量、开采量等。下面将以弹性释水量与GH B边界水量为例,阐述基于一阶摄动法的地下水均衡项区间的数值模拟方法。2.2弹性释水量在G FM o d e l中,承压水非稳定数值模型中第k个三角形的待求节点j在第t个时段末的弹性释水量平均值QjS 0的表达式为:QjS 0=Sk0Ak(Hjt-Hjt-1)(4)式中:Sk0为第k个三角形单元的贮水率平均值(1/m);Sk0/Ak为第k个三角形单元内待求节点j的弹性释水量下的单元系数矩阵的系数(m2);Ak中包含三角形的尺寸信息,且不随时间变化(m-3);Hjt为节点j在第t个时段末的水头(m);Hjt-1为节点j在第t-1时段末的水头(m)。当贮水率为区间时,此时弹性释水量可以表示为:QjS=1Ak(Sk0+Sk)(Hjt+Hjt)-(Hjt-1+Hjt-1)(5)式中:Sk为贮水率的变化量(1/m);Hjt、Hjt-1分别为节点j在第t个时段末、第(t-1)个时段末的水头变化量(m)。式(5)称作是单元节点弹性释水量的一阶摄动展开。展开式(5),并忽略其中的二阶及以上的高阶项,得到:QjS=1Ak(SkHjt+Sk0Hjt-SkHjt-1-Sk0Hjt-1)(6)式(6)中,除了Hjt与Hjt-1之外,其余均为已知项,可直接计算得出。而Hjt与Hjt-1可根据文献2 8 中的方法得到,QjS极值的一般形式为:(QjS)m a x=|m|(m)m a x(7)(QjS)m i n=-|m|(m)m a x(8)m=1,2,p式中:p为节点j水平相连的所有三角形中渗透系数分区的个数加1,这里的加1就是指当前三角形的贮水率;m是p个参数的相应偏差;m由第k个三角形的各类单元渗透矩阵和列向量以及系数矩阵的逆矩阵计算得到。由此,第k个三角形中待求节点j在第t个时段末的弹性释水量区间为QjS-|m|(m)m a x,QjS+|m|(m)m a x,而整个研究区内的弹性释水量区间只需将所有待求节点的弹性释水量区间进行累加即可。2.3GH B边界水量GH B(g e n e r a l-h e a d b o u n d a r y)是地下水数值模型中使用非常广泛的边界条件,当边界不是定水位或者给定流量边界的时候,一般会使用该边界条件,该边界更接近第三类边界条件,本质上是通过达西公式的方式计算边界处的水量,G FM o d e l中也实现了这 一 类 边 界 条 件。承 压 水 确 定 性 数 值 模 型 中 GH B边界上某一节点j在第t个时段末的水量交换平均值QjB 0可表示为:QjB 0L-1K(Hj0-Hjt)Bt(9)式中:L-1为距离倒数(1/m);K为GH B边界渗透系数(m/d);Hj0为GH B边界在节点j对应的水头平均值(m);Hjt为第t个时段末节点j的水头(m);B为过水断面面积(m2),仅与含水层厚度和单元剖分尺寸有关,不随时间变化。当GH B边界中的距离倒数、渗透系数、水头发生变化时,相应的水量QjB则为:QjB=(L-1+L)(K+k)(Hj0+H)-(Hjt+Hjt)Bt(1 0)式中:L、k、H分别为距离倒数(1/m)、渗透系数(m/d)、GH B水头(m)的变化量;Hjt为节点j在第t个时段末的水头变化量(m)。同样地,展开式(1 0),并忽略其中的二阶及以上的高阶项,有:QjB=(L KHj0+L-1k Hj0+L-1KH-L KHjt-L-1KHjt-L-1k Hjt)Bt(1 1)在式(1 1)中,除L-1KHjt之外,其余均可直接计算得出,对于L-1KHjt,根据文献2 8 中的方法得到。同样,GH B边界水量区间偏差的上下限可以表示为:(QjB)m a x=|m|(m)m a x(1 2)(QjB)m i n=-|m|(m)m a x(1 3)待求节点j在第t个时段末的GH B边界水量区间 即 为 QjB 0-|m|(m)m a x,QjB 0+|m|(m)m a x,累加GH B边界上的所有节点的水量区间可得到整个研究区内GH B变界上水量的区间。其余水均衡项的区间计算与弹性释水量、GH B边界水量的计算方法类似,都可以得到与式(7,8)或者式(1 2,1 3)相同的表达式,区别是不同的水均衡项的m的计算是不同的,但都会用到三角形单元的87 第4期董贵明等:二维承压非稳定流水均衡区间的数值模拟渗透矩阵、列向量和整体系数矩阵的逆矩阵。2.4误差的检验在前边给出的弹性释水量和GH B边界水量的区间模拟方法是通过一阶摄动展开并忽略二阶及以上高阶项得到的,所以是一种近似方法。这里采用等间距连续采样区间分析方法对本研究方法进行误差检验。等间距连续采样区间分析方法本质上是一种穷举法,它从水文地质参数的区间中以较小的等间隔连续取样,采用确定性模型计算每一个样本值下的各个水均衡值,当采样完成的时候,就可以得到每个水均衡项的若干个值,从这些值中挑选出最大值和最小值构成该水均衡项的区间。当采样间隔越小的时候,得到的水均衡项的区间越接近理论区间。文献2 8 中采用了等间距连续采样区间分析方法进行水位区间的误差检验。因为对每一个采样值都要调用一次确定性模型,在取样次数非常多的情况下,计算地下水均衡区间可能要花费几天甚至更长的时间(使用一般的个人计算机),因此该方法计算效率非常低。等间距连续采样区间分析计算水均衡区间的步骤为:在参数i的区间I=i,i(i=1,2,n,n为参数的个数)中,以相等的距离连续采样,得到ni个节点。其中,每个节点为参数i在区间I中的取值,即参数值。而参数区间I此时被ni个节点分割成(ni-1)个长度相等的子区间;从ni个参数值中提取一个参数值,共可提取出n个参数值,接着可计算出在m个参数值下的水均衡项,共可得到ni=1ni个地下水均衡项值;从ni=1ni个水均衡项值中选择出极大值和极小值,分别作为等间距水均衡项区间的极大值与极小值,至此,等间距水均衡项的区间即可获得。对于参数区间I=i,i ,其中的子区间长度可以通过试算的方法确定。一般地,子区间长度为区间长度的5%或者更小,当子区间长度越小,方法的计算精度越高,但计算的时间也明显越长。3 算例分析建立的理想模型范围为1 0 0 m1 0 0 m,含水层水平等厚、非均质各向同性。将含水层边界设置为2种,分别为定水头边界与GH B边界,其中左、右边界及上边界为定水头边界,下边界为GH B边界(图2)。在模型中设置一抽水井W-1,坐标为(3 0 m,4 0 m),流量为5 0 0 m3/d。将理想模型分为5个区域分别赋予不同的渗透系数值,并在每个分区中设置一个观测井,共有5个观测井,编号从O-1至O-5(图2)。图的左下角作为原点(0 m,0 m)。图2 模型渗透系数分区图F i g.2 Z o n a t i o n o f h y d r a u l i c c o n d u c t i v i t i e s f o r t h e g r o u n d w a t e r f l o w m o d e l模型参数值见表1,各分区渗透系数值见表2,各观测井的坐标见表3。表1 理想模型含水层各项参数设定T a b l e 1 A q u i f e r p a r a m e t e r s e t t i n g o f t h e h y p o t h e t i c a l m o d e l参数参数值参数参数值含水层厚度/m2 0定水头边界水头/m7 0初始水头/m8 0顶板高程/m2 0底板高程/m0贮水率/m-111 0-4GH B边界水头/m8 0渗透系数/(md-1)6距离倒数/m-10.0 1表2 各分区渗透系数数值T a b l e 2 H y d r a u l i c c o n d u c t i v i t y s e t t i n g o f e a c h z o n e分区号12345渗透系数/(md-1)74953表3 各分区观测井信息T a b l e 3 O b s e r v a t i o n w e l l i n f o r m a t i o n f o r e a c h z o n e观测井编号X/mY/m观测井编号X/mY/mO-12 57 8O-47 44 5O-22 12 6O-56 67 2O-36 12 0首先根据表1,2和3中的数值进行地下水均衡项确定性模型计算,计算结果如表4所示。表中模拟了0.1 d内地下水各水均衡项的变化,分别列出了各时段的水均衡值。表4中的均衡误差是通过G FM o d e l中的水均衡模块计算出来的,均衡误差为0可以在一定程度上说明水位水均衡项的计算结果是正确的。均衡误差与单元剖分、非均质性、地下水位空间分布、方程97h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 表4 理想算例水均衡项(水量单位:m3)T a b l e 4 G r o u n d w a t e r b u d g e t i n t h e h y p o t h e t i c a l m o d e l时段抽水量GHB边界侧向补给量定水头边界侧向补给量弹性释水量均衡误差/1 0-21-5.0 05.0 6-8 8.9 1-8 8.8 50.0 02-5.0 07.2 4-4 0.3 3-3 8.0 90.0 03-5.0 08.4 8-2 3.5 9-2 0.1 10.0 04-5.0 09.2 2-1 5.4 9-1 1.2 60.0 05-5.0 09.6 7-1 1.1 0-6.4 30.0 06-5.0 09.9 3-8.6 3-3.7 00.0 07-5.0 01 0.0 9-7.2 2-2.1 30.0 08-5.0 01 0.1 8-6.4 1-1.2 30.0 09-5.0 01 0.2 3-5.9 4-0.7 10.0 01 0-5.0 01 0.2 6-5.6 7-0.4 10.0 0总计-5 0.0 09 0.3 5-2 1 3.3 0-1 7 2.9 50.0 0 注:1.水量流入模拟区为+,流出为-;2.均衡误差=边界侧向补给量(GH B边界+定水头边界)-弹性释水量+抽水量组求解精度要求等有关。下面将5个分区的渗透系数、定水头边界的水头、GH B边界的渗透系数以及贮水率共8个参数作为区间变量,计算参数变化率分别为0.1,0.2,0.3时0.1 d内地下水各个水均衡项的区间。区间变量的变化率为区间的偏差除以区间中点值,区间中点来自于表1中的值。比如,某个参数的平均值为1 0,那么当参数的区间为8,1 2 的时候,其变化率为0.2,1 0也可以称为区间中点值。计算结果如图3所示。图中分别给出了本文提出的方法和等间距连续采样方法计算出的区间上限(极大值)、下限(极小值),为了便于对比,也给出了确定性模型的结果。从图3可以看出,在各变化率下,水均衡区间的极值可以很好地包含确定型模型的水均衡值,各时段 内 水 均 衡 区 间 随 着 时 间 延 长 逐 渐 趋 于 稳 定。GH B边界水量区间先扩大后趋于稳定,定水头边界水量区间与弹性释水量区间都是逐渐变小后趋于稳定,各种水量区间模拟后期均趋于稳定的原因是流场逐渐趋于稳定流。对比本研究方法计算的水均衡区间与等间距连续采样法计算的水均衡区间,在同一变化率下,区间极大值、极小值的相对误差随模拟时间的增加而增加,当流场趋于稳定时,区间极值相对误差也基本稳定不变。本算例中,在渗透系数的变化率不超过0.1时,3种水均衡项的偏差相对误差均可以控制在1 0%以内,而区间极值的相对误差可以控制在5%以内,方法具有一定的精度。a,b,c.变化率为0.1;d,e,f.变化率为0.2;g,h,i.变化率为0.3图3 不同时段地下水均衡区间对比F i g.3 C o m p a r i s o n o f t h e g r o u n d w a t e r b u d g e t i n t e r v a l i n d i f f e r e n t t i m e08 第4期董贵明等:二维承压非稳定流水均衡区间的数值模拟 在计算过程中,分别计算了等间距连续采样法的采样次数为2,4,6,8这4种情况,结果发现求出的各水均衡项的区间上下限是相同的。这是因为本算例的承压水系统中各个区间参数与各水均衡项之间是单调关系,所以各水均衡项极值均在区间变量的两端处获得,这时确定性模型的计算次数为28。但并不是所有地下水系统的水均衡极值都能在区间变量的两端处获得。算例中只有8个参数为区间变量。随着模型中区间变量个数的增加,等间距连续采样法运行确定性模型的次数是指数增加的,但由于区间数值模拟的主要计算量是系数矩阵求逆,所以,当仅增加模型中区间变量个数的时候,本研究提出的方法的计算时间增加是十分有限的。4 结 论(1)基于一阶摄动法的承压水水均衡项区间模拟方法的数值算例分析表明,在区间变量的变化率不超过0.1时,3种水均衡项的偏差相对误差均可以控制在1 0%以内,而区间极值的相对误差可以控制在5%以内,方法具有一定的精度。因此,在具体应用时,区间变量的变化率一般应不超过0.1。(2)为了计算水均衡项的区间,模型输入变量(渗透系数等)的区间准确性也是至关重要的。另外,区间不确定性数值模拟的结果中不包含任何统计信息。在以后的研究中,应加强区间和随机等传统不确定性分析方法的结合。(所有作者声明不存在利益冲突)参考文献:1 S r e e k a n t h J,M o o r e C.N o v e l p a t c h m o d e l l i n g m e t h o d f o r e f f i-c i e n t s i m u l a t i o n a n d p r e d i c t i o n u n c e r t a i n t y a n a l y s i s o f m u l t i-s c a l e g r o u n d w a t e r f l o w a n d t r a n s p o r t p r o c e s s e sJ.J o u r n a l o f H y d r o l o g y,2 0 1 8,5 5 9:1 2 2-1 3 5.2 薛佩佩,文章,梁杏.地质统计学在含水层参数空间变异研究中的应用进展与发展趋势J.地质科技通报,2 0 2 2,4 1(1):2 0 9-2 2 2.X u e P P,W e n Z,L i a n g X.A p p l i c a t i o n a n d d e v e l o p m e n t t r e n d o f g e o s t a t i s t i c s i n t h e r e s e a r c h o f s p a t i a l v a r i a t i o n o f a q u i f e r p a-r a m e t e r sJ.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