含水层非均质性对地下水蒙特卡罗模拟结果的影响

王晶晶,樊尊荣

(1.江苏省水文水资源勘测局扬州分局,江苏 扬州 225000; 2.江苏兴水建设工程有限公司,江苏 南京 210001)



含水层非均质性对地下水蒙特卡罗模拟结果的影响

王晶晶1,樊尊荣2

(1.江苏省水文水资源勘测局扬州分局,江苏 扬州 225000; 2.江苏兴水建设工程有限公司,江苏 南京 210001)

为研究含水层非均质性对地下水蒙特卡罗模拟结果的影响,对非均质性分区渗透系数最终随机场的产生采用对应充填和排列充填两种处理方法,并对两种处理方法的蒙特卡罗模拟结果进行了比较。结果表明:对于具有强烈非均质性的含水层,蒙特卡罗模拟的实现数对于对应充填至少要取至2 500次,排列充填至少要取至10 000次,所得结果才是合理的;在绝对实现数相同的情况下,对应充填处理结果优于排列充填,但当实现数达到22 500次以上时,两种处理方法所得结果就无明显差别。

含水层;地下水;非均质性;蒙特卡罗模拟;对应充填;排列充填

地下水数值模拟是研究地下水最主要的方法和手段[1-2],确定性方法是进行地下水数值模拟的传统方法,目前已比较成熟,在解决实际问题上取得了很好的成果[3-4]。对一个实际问题建立模型最困难的往往是水文地质参数的确定,对于水文地质参数的处理,传统的确定性方法均为在研究区内根据实际的水文地质条件进行参数分区,并且认为同一参数分区内的水文地质参数是相同的。但由于天然土壤和含水层沉积过程的随机性,含水层的水文地质特征以及土壤的结构、构造和土壤各种矿物的组成等具有非常明显的空间变异性,含水介质的空间变异性是通过水文地质参数的空间变异性表现出来的。

众多研究成果表明,水文地质参数尤其是渗透系数的空间变异性是影响溶质运移的决定性因素。目前随机方法已成为研究非均质含水层中地下水流和溶质运移问题的重要手段[5-10]。地下水随机模拟方法主要有矩方程法和蒙特卡罗法。矩方程法通过求解有关均值和协方差的随机偏微分方程组获得随机问题的解,而蒙特卡罗方法则是通过平均一系列反映含水层实际性质的确定性问题来模拟随机过程的一种计算机模拟方法。采用蒙特卡罗方法进行地下水流和溶质运移的数值模拟具有易于理解、便于操作的优点,已被越来越多的研究人员接受[11-12]。陈彦等[13]用蒙特卡罗方法研究了含水层渗透系数的空间变异性对地下水数值模拟结果的影响;阎婷婷等[14]同样用蒙特卡罗方法研究了渗透系数的均值和方差对污染羽的二维空间分布的影响。上述研究涉及的含水层虽然有空间变异性,但相对还是比较均匀的(lnK的方差一般小于0.5,采用渗透系数符合对数正态分布的假设[15])。野外实际含水层的非均质性是非常复杂的,非均质性分区处理一般不可避免,此外,在砂层中夹有较大的亚砂土或亚黏土透镜体也是很常见的,当研究尺度不是很大时,不能简单做统计平均处理,必须按非均质性分区处理。本文研究含水层非均质性分区对地下水蒙特卡罗模拟结果的影响,采用对应充填和排列充填两种方法在不同实现数下产生非均质性分区渗透系数最终随机场,并对两种方法蒙特卡罗模拟的结果进行比较。本文研究不考虑研究区剖分单元的密度对模拟结果的影响,因此,所有模拟计算对研究区均采用同一种剖分密度。

1 研究区概况

研究区是在Borden试验场局部含水层的基础上进行改造得到的,为了体现含水层的非均质性,在原Borden试验场中加入了一个矩形高透水性小区域(渗透系数为K2,其他区域渗透系数为K1),以模拟分析三维的非均质各向异性承压含水层中的溶质运移问题。含水层的平面图如图1所示,整个研究区东西向距离为19 m,南北向距离为12 m,含水层厚为1.8 m。其中有一个相对高透水性的小区,其平面上东西向长为4.5 m,南北向宽为5.5 m,垂向上贯穿整个含水层厚度。

图1 研究区平面示意图(单位:m)

对于水文地质参数的确定,除了相对高透水性小区的渗透系数K2外,其他水文地质参数均采用Borden试验场局部含水层的数据[16-17]:含水层在水平方向的相关长度为2.8 m,垂向的相关长度为0.12 m;含水层孔隙度为0.1,储水率为0.25 mL/ m。本文研究基于渗透系数服从对数正态分布的假设,在Borden含水层资料的基础上突出大小区的透水性不同,在结果合理的前提下,假设lnK2的均值为4.605 170,方差为2.302 585,lnK1的均值为1.609 438,方差为0.3 147 165(即原Borden试验场的数值)。

含水层的左右两边边界设为定水头边界,水头分别为20.0 m和19.8 m，水流方向自左向右。总的模拟时间为9 d。假设图1中Q点为一口不完整注水井,井深达到地面以下0.5 m,注水速率随时间变化,在前3天为0.5 L/d,第4～6天为0.8 L/d,第7～9天为0.7 L/d。注入的水中,溶质的质量浓度为100 g/L。由于主要分析含水层的非均质性分区对地下水蒙特卡罗模拟的影响,故仅考虑渗透系数的空间变异性,不考虑孔隙度的空间变异性。

2 蒙特卡罗模拟

2.1 渗透系数场的产生

采用Robin等[18]提出的直接傅立叶变换方法来生成单一的渗透系数随机场,对于本文的高度非均质性研究区(即渗透系数变化大)必须进行非均质分区,但如何产生最终的渗透系数随机场目前还没有确定的方法,本文分别采用两种可能的方法得到最终渗透系数随机场来进行比较。

首先将研究区按非均质性进行分区,外围的大矩形区定为Ⅰ区,内部的小矩形区定为Ⅱ区;然后分别按前述数据(lnK1的均值为1.609 438,方差为0.3 147 165;lnK2的均值为4.605 170,方差为2.302 585)产生Ⅰ和Ⅱ区渗透系数随机场;最后采用两种可能的拼接方法得到最终的研究区渗透系数随机场。

将渗透系数随机场Ⅱ充填进渗透系数随机场Ⅰ是得到最终研究区渗透系数随机场的基础,而本文采用的两种拼接方法的区别也就在于充填方法的不同,下面以3次实现数为例说明这两种拼接方法。

假设已经分别得到了随机场Ⅰ和Ⅱ的3个不同实现(K1的随机场1、K1的随机场2、K1的随机场3和K2的随机场1、K2的随机场2、K2的随机场3),然后通过以下两种充填拼接方法得到最终的研究区渗透系数K的随机场:

a. 对应充填。将K2的3个随机场按照实现数一一对应充填到K1随机场中的相应位置,即

由此得到3个最终的研究区渗透系数K的随机场。

b. 排列充填。将K2的3个随机场按照实现数用排列的方法充填到K1的随机场中,即

由此得到9个最终的研究区渗透系数K的随机场。

2.2 蒙特卡罗模拟方法和步骤

蒙特卡罗方法是一种基于计算机模拟的统计试验方法,首先利用伪随机数生成足够多组服从给定分布规律的随机变量(如渗透系数),然后对每组输入变量分别进行数值模拟。本文研究含水层非均质性对蒙特卡罗模拟结果的影响大致分为以下5个步骤:①利用直接傅立叶变换方法分别在Ⅰ区和Ⅱ区空间离散的网格节点上随机产生符合给定分布特征的渗透系数场,即完成渗透系数场的1次实现;②用两种不同充填拼接方法得到最终的研究区渗透系数场的1次实现;③根据生成的最终渗透系数场,运行对应条件下的地下水流模型和溶质运移模型,以获得整个研究区的地下水流场和溶质的浓度分布场;④由确定的实现次数重复前3步;⑤最后对所有实现的结果进行统计。

将研究区离散为25行、39列、36层,其中小矩形区占据了第7～18行、第10～19列,垂向上为36层。水平面上沿x和y方向网格的间距均为0.5 m,垂向上网格间距为0.05 m。使用地下水流模拟模型MODFLOW[19]模拟地下水的流动,得到地下水流场,然后利用溶质运移模型MT3DMS[20]追踪溶质质点随水流的运动,最终得到研究区的浓度分布场。为了充分比较两种充填方法的结果,在采用蒙特卡罗方法进行随机模拟时,对应充填分别尝试了20次、50次、100次、150次和200次实现,相应的排列充填实现数分别为400次、2 500次、10 000次、22 500次和40 000次;另外为了更充分地比较模拟结果,还用对应充填的方法计算了400次、2 500次、10 000次、22 500次和40 000次实现的结果。

3 模拟结果及分析

3.1 水头场的比较

通过地下水流模型MODFLOW求解可以得到第3天、第6天和第9天的地下水流场。图2给出了150次相关实现数第9天的水头均值等值线不同方法模拟结果,可以看出两种方法以及不同实现数间差别不明显。这是由于本次模拟中只有单口井,注水量很小,且模拟时间短(只有9 d),所以模拟过程中研究区水头的变化都不大。

图2 150次相关实现数第9天的水头均值等值线对比(单位:m)

此外,从图2还可以看出,中间相对高透水性区域水头均值等值线排列疏松,即水力坡度下降较慢;两边水头均值等值线相对较密,水力坡度下降快,即为相对渗透系数小的区域,因此模拟得到的地下水流场是合理的。

3.2 水头方差的比较

图3 150次相关实现数第9天的水头方差等值线对比(单位:m2)

用蒙特卡罗方法计算出的水头和浓度都是存在方差的。由图3可以看出两种处理方法得到的水头方差的差别主要表现在中间相对高透水性小区域中,而且较大的方差也是出现在小区域中,最大方差值出现在井点附近。总体而言,排列充填的方差范围要比对应充填的方差范围小,而且随实现数的增加,方差范围逐渐缩小。从图3可以看出,当对应充填的实现数增加到400次以后,所得同样大小的方差范围比相同实现数的排列充填的范围小,也就是说所得水头均值变幅更小,即更准确。浓度方差亦有相同规律,具体见下文分析。

3.3 质量浓度均值的比较

使用特征线法追踪溶质质点随水流的运动,可以获得研究区溶质在不同模拟时间的质量浓度分布场,为简便起见,本文仅对最顶层的质量浓度均值等值线进行比较。

为了分析计算结果的可靠性,仍按传统方法将同一参数分区里的渗透系数取为同一数值,即将渗透系数按均值处理进行一次地下水流与溶质运移的模拟(以下简称为按均值处理),然后与蒙特卡罗方法进行比较。根据蒙特卡罗模拟所得的均值和方差结果可以绘制出蒙特卡罗计算结果的波动范围,这一波动范围代表了各种可能性的计算结果,也包括按均值处理的情况;因此,如果蒙特卡罗的计算结果是合理的,那么按均值处理的结果应包含在蒙特卡罗计算结果的波动范围之中。将两种方法不同实现数的蒙特卡罗计算结果波动范围与按均值处理的结果进行对比分析表明,只有当对应充填的绝对实现数取到2 500次以上,排列充填的绝对实现数取到10 000次以上时,按均值处理的结果才会完全落在蒙特卡罗计算结果的波动范围内(图4)。由此可见,对于这种具有强烈非均质性的含水层进行蒙特卡罗模拟时,以往经常采用的几百次的实现数是远远不够的,对应充填实现数要取到2 500次以上,而排列充填实现数要取到10 000次以上时所得的结果才是合理的。

图4 蒙特卡罗模拟与按均值处理运移5 d的质量浓度均值对比(单位:g/L)

另外对比蒙特卡罗模拟的结果和按均值处理的结果,可以发现按均值处理的结果质量浓度均值等值线主要为长条形,质点在纵向上随水流运移明显,而蒙特卡罗模拟得到的质量浓度均值等值线则在横向上亦有明显的分布,且在纵向上的扩散要比按均值处理慢,由于水流方向是沿纵向自左向右,可见蒙特卡罗方法的计算结果既有纵向上因对流引起的溶质扩散,又较好地体现了溶质质点在横向上的弥散现象,蒙特卡罗方法能更好地体现含水层介质的空间变异性。

比较不同实现数不同运移时间下两种处理方法的质量浓度均值等值线(图5),可以看出对应充填方法和排列充填方法得到的质量浓度场存在着明显的差别,且随着运移时间的增长,区别也愈加明显,但是随实现数的增加,二者区别则相对减小。排列充填在实现数很少(100次以内)的情况下得出的结果显然是不合理的,而当对应充填的实现数增加至与排列充填的绝对实现数相同时,总体上溶质在纵向上的扩散排列充填要比对应充填快。由图5可以看出,在实现数低于10 000次的情况下,两种方法得到的质量浓度均值等值线有差别,主要表现在纵向上的扩散快慢,而当实现数达到22 500次以上时,两种方法所得结果就无明显差别。

图5 运移9 d质量浓度均值等值线对比(单位:g/L)

当对应充填的实现数为20次、50次、100次、150次和200次时,排列充填的绝对实现数为400次、2 500次、10 000次、22 500次和40 000次,此时两种方法所产生的基础K1和K2的随机场是相同的,但是对于排列充填而言等于是增加了蒙特卡罗模拟实现的次数。根据随机理论,增加实现数最终结果的波动会变小,排列充填得到的结果会更准确。当对应充填的实现数增加到2 500次之后,所得的结果比相同绝对实现数下排列充填得到的结果更准确,因为此时排列充填的绝对实现数虽然和对应充填相同,但是其并未增加基础K1和K2的随机场数目,而对应充填在增加了实现数的同时还增加了产生的K1和K2的随机场数目,这样得到的最终K的随机场更为符合实际情况,如表1所示。随实现数的增加,所得到的K1和K2随机场的均值和方差越接近给定的均值和方差,误差越小,得到的最终渗透系数场也就更为准确,因而所得到的质量浓度结果亦更为准确(前述水头方差随实现数增加而变小原因同此)。需要注意的是当绝对实现数达到22 500次以上时,两种方法所得结果就无明显差别。

表1 不同实现数时lnK1和lnK2的均值与方差

3.4 质量浓度方差的比较

模拟结果表明,质量浓度方差的分布规律与质量浓度均值的分布规律相同。当对应充填的实现数为20次、50次和100次时,质量浓度方差对应于质量浓度值也存在一个狭长带,随实现数的增加,该狭长带也就逐渐缩短并消失。质量浓度方差的最大值出现在井点附近,并且随着对应充填实现数的增加,质量浓度方差范围明显变小,质量浓度方差最大值也随实现数的增加而减小。

两种方法的模拟得到的质量浓度方差分布差别与均值差别的规律类似,随实现数的增加两者的区别减小,当实现数达到22 500以上时,两者的差别就很小。需要指出的是,对应充填所得的最大方差值要小于排列充填所得的最大方差值。

表2为200次相关实现数不同方法模拟得到的质量浓度方差最大值比较,在没有增加对应充填实现次数时,排列充填所得质量浓度方差最大值是小于对应充填的,但是当对应充填的实现次数增加到与排列充填相同后,质量浓度方差值就小于排列充填了,方差值变小,则所得质量浓度值变幅就变小,结果也就更符合实际。

表2 200次相关实现数不同运移时间质量浓度方差最大值

4 结 论

对产生的符合研究区统计特性的随机场多次实现,进行蒙特卡罗模拟可以充分考虑含水层水力参数的结构性和随机性,能较真实地刻画含水层中溶质运移的过程,但是对含水层非均质性的不同处理方法对蒙特卡罗模拟所得到结果是有影响的。对于本文研究的具有强烈非均质性的含水层,对应充填的蒙特卡罗模拟绝对实现数要取至2 500次,排列充填的蒙特卡罗模拟绝对实现数要取至10 000次,才能得到合理的计算结果。

[ 1 ] 薛禹群,谢春红.水文地质学的数值法[M].北京:煤炭工业出版社,1980.

[ 2 ] 薛禹群,谢春红,吴吉春.水文地质数值法存在的问题及其对策[J].地球科学进展,1996,11(5):472-474.(XUE Yuqun,XIE Chunhong,WU Jichun.Problems in hydrogeological numerical method and its countermeasures[J].Advances in Earth Science,1996,11(5):472-474.(in Chinese))

[ 3 ] 薛禹群,吴吉春.地下水数值模拟在我国:回顾与展望[J].水文地质工程地质,1997(4):21-24.(XUE Yuqun,WU Jichun.The numerical simulation of groundwater in China:retrospect and prospect[J].Hydrogeology and Engineering Geology,1997(4):21-24.(in Chinese))

[ 4 ] 薛禹群,吴吉春.面临21世纪的中国地下水模拟问题[J].水文地质工程地质,1999(5):1-3.(XUE Yuqun,WU Jichun.Problems of groundwater modeling in China facing 21st century[J].Hydrogeology and Engineering Geology,1999(5):1-3.(in Chinese))

[ 5 ] 杨金忠,蔡树英,黄冠华,等.多孔介质中水分及溶质运移的随即理论[M].北京:科学出版社,2000.

[ 6 ] 杨金忠,蔡树英,叶自桐.区域地下水溶质运移随机理论的研究与进展[J].水科学进展,1998,9(1):84-98.(YANG Jinzhong,CAI Shuying,YE Zitong.Research and development of the regional groundwater solute transport random theory[J].Advances in Water Science,1998,9(1):84-98.(in Chinese))

[ 7 ] DAGAN G.Flow and transport in porous formations[M].New York:Springer-Verlag,1989.

[ 8 ] DAGAN G.Stochastic modeling of groundwater flow by unconditional and conditional probabilities:conditional simulation and the direct problem[J].Water Resources Research,1982,18(4):813-833.

[ 9 ] NEUMAN S P,ZHANG Y K.A quasi-linear theory of non-Fickian and Fickian subsurface dispersion:theoretical analysis with applicaton to isotropic media[J].Water Resources Research,1990,26(5):887-902.

[10] WU J C,HU B X,HE C.A numerical method of moments for solute transport in a porous medium with multiscale physical and chemical heterogeneity[J].Water Resources Research,2004,40(1),W01508.

[11] AHMED S,MARSILY G D.Comparision of geostatistical methods for estimating transmissivity using data in transmissivity and specific capacity[J].Water Resources Research,1987,23(9):1717-1737.

[12] HASSAN A E,CUSHMAN J H,DELLEUR J W.A Monte Carlo assessment flow and transport perturbation models[J].Water Resources Research,1998,34(5):1143-1163.

[13] 陈彦,吴吉春.含水层渗透系数空间变异性对地下水数值模拟的影响[J].水科学进展,2005,16(4):428-487.(CHEN Yan,WU Jichun.Effect of the spatial variability of hydraulic conductivity in aquifer on the numerical simulation of groundwater[J].Advances in Water Science,2005,16(4):428-487.(in Chinese))

[14] 阎婷婷,吴剑锋.渗透系数的空间变异性对污染物运移的影响研究[J].水科学进展,2006,17(1):29-36.(YAN Tingting,WU Jianfeng.Impacts of the spatial variation of hydraulic conductivity on the transport fate of contaminant plume[J].Advances in Water Science,2006,17(1):29-36.(in Chinese))

[15] FREEZE R A.A stochastic conceptual analysis of one-dimensional groundwater flow in non-uniform homogeneous media[J].Water Resources Research,1975,11(5):725-741.

[16] SUDICKY E A.A natural gradient experiment on solute transport in a sand aquifer:spatial variability of hydraulic conductivity and its role in the dispersion process[J].Water Resources Research,1986,22(13):2069-2082.

[17] HUANG H,HASSAN A E,HU B X.Monte Carlo study of conservative transport in heterogeneous dual-porosity media[J].Journal of Hydrology,2003,275:229-241.

[18] ROBIN M J L,GUTJAHR A L,SUDICKY E A,et al.Cross-correlated random field generation with the direct Fourier transform method[J].Water Resources Research,1993,29(7):2385-2397.

[19] MCDONALD M G,HARBAUGH A W.A modular three-dimension finite-difference ground water model[R].Reston,VA,USA:USGS,1988.

[20] ZHENG C,WANG P P.A modular three-dimension multispecies transport model for simulation of advection,dispersion and chemical reactions of contaminants in ground water systems:documentation and user’s guide [R].New York:ERDC,1999.

Impacts of aquifer heterogeneity on Monte Carlo simulation of groundwater

WANG Jingjing1，FAN Zunrong2

(1.YangzhouBranchofJiangsuProvincialHydrologyandWaterResourcesBureau,Yangzhou225000,China; 2.JiangsuXingshuiConstructionEngineeringCo.,Ltd.,Nanjing210000,China)

In order to study the impacts of aquifer heterogeneity on the results of Monte Carlo simulation of groundwater, two methods, i.e., the one-to-one filling method and the arrangement filling method, were used to obtain the permeability coefficient random field, and the results of Monte Carlo simulation using the two methods were compared. The study suggests that, for the aquifer with strong heterogeneity, the necessary number of realizations of Monte Carlo simulation should be at least 2 500 for the one-to-one filling method and 10 000 for the arrangement filling method, in order to obtain reasonable results. In the case of consistent absolute number of realizations, the one-to-one filling method performs better than the arrangement filling method. However, when the number of realizations reaches 22 500 or above, there is no significant difference between the two methods.

aquifer; groundwater; heterogeneity; Monte Carlo simulation; one-to-one filling method; arrangement filling method

10.3880/j.issn.1004-6933.2017.01.010

国家自然科学基金(40272106)

王晶晶(1982—),女,工程师,硕士,主要从事地下水及水文水资源研究。E-mail:360612133@qq.com

P641

A

1004-6933(2017)01-0046-06

2016-04-26 编辑：熊水斌)