王晶晶,樊尊榮

(1.江蘇省水文水資源勘測(cè)局揚(yáng)州分局,江蘇 揚(yáng)州 225000; 2.江蘇興水建設(shè)工程有限公司,江蘇 南京 210001)

?

含水層非均質(zhì)性對(duì)地下水蒙特卡羅模擬結(jié)果的影響

王晶晶1,樊尊榮2

(1.江蘇省水文水資源勘測(cè)局揚(yáng)州分局,江蘇 揚(yáng)州 225000; 2.江蘇興水建設(shè)工程有限公司,江蘇 南京 210001)

為研究含水層非均質(zhì)性對(duì)地下水蒙特卡羅模擬結(jié)果的影響,對(duì)非均質(zhì)性分區(qū)滲透系數(shù)最終隨機(jī)場(chǎng)的產(chǎn)生采用對(duì)應(yīng)充填和排列充填兩種處理方法,并對(duì)兩種處理方法的蒙特卡羅模擬結(jié)果進(jìn)行了比較。結(jié)果表明:對(duì)于具有強(qiáng)烈非均質(zhì)性的含水層,蒙特卡羅模擬的實(shí)現(xiàn)數(shù)對(duì)于對(duì)應(yīng)充填至少要取至2 500次,排列充填至少要取至10 000次,所得結(jié)果才是合理的;在絕對(duì)實(shí)現(xiàn)數(shù)相同的情況下,對(duì)應(yīng)充填處理結(jié)果優(yōu)于排列充填,但當(dāng)實(shí)現(xiàn)數(shù)達(dá)到22 500次以上時(shí),兩種處理方法所得結(jié)果就無(wú)明顯差別。

含水層;地下水;非均質(zhì)性;蒙特卡羅模擬;對(duì)應(yīng)充填;排列充填

地下水?dāng)?shù)值模擬是研究地下水最主要的方法和手段[1-2],確定性方法是進(jìn)行地下水?dāng)?shù)值模擬的傳統(tǒng)方法,目前已比較成熟,在解決實(shí)際問(wèn)題上取得了很好的成果[3-4]。對(duì)一個(gè)實(shí)際問(wèn)題建立模型最困難的往往是水文地質(zhì)參數(shù)的確定,對(duì)于水文地質(zhì)參數(shù)的處理,傳統(tǒng)的確定性方法均為在研究區(qū)內(nèi)根據(jù)實(shí)際的水文地質(zhì)條件進(jìn)行參數(shù)分區(qū),并且認(rèn)為同一參數(shù)分區(qū)內(nèi)的水文地質(zhì)參數(shù)是相同的。但由于天然土壤和含水層沉積過(guò)程的隨機(jī)性,含水層的水文地質(zhì)特征以及土壤的結(jié)構(gòu)、構(gòu)造和土壤各種礦物的組成等具有非常明顯的空間變異性,含水介質(zhì)的空間變異性是通過(guò)水文地質(zhì)參數(shù)的空間變異性表現(xiàn)出來(lái)的。

眾多研究成果表明,水文地質(zhì)參數(shù)尤其是滲透系數(shù)的空間變異性是影響溶質(zhì)運(yùn)移的決定性因素。目前隨機(jī)方法已成為研究非均質(zhì)含水層中地下水流和溶質(zhì)運(yùn)移問(wèn)題的重要手段[5-10]。地下水隨機(jī)模擬方法主要有矩方程法和蒙特卡羅法。矩方程法通過(guò)求解有關(guān)均值和協(xié)方差的隨機(jī)偏微分方程組獲得隨機(jī)問(wèn)題的解,而蒙特卡羅方法則是通過(guò)平均一系列反映含水層實(shí)際性質(zhì)的確定性問(wèn)題來(lái)模擬隨機(jī)過(guò)程的一種計(jì)算機(jī)模擬方法。采用蒙特卡羅方法進(jìn)行地下水流和溶質(zhì)運(yùn)移的數(shù)值模擬具有易于理解、便于操作的優(yōu)點(diǎn),已被越來(lái)越多的研究人員接受[11-12]。陳彥等[13]用蒙特卡羅方法研究了含水層滲透系數(shù)的空間變異性對(duì)地下水?dāng)?shù)值模擬結(jié)果的影響;閻婷婷等[14]同樣用蒙特卡羅方法研究了滲透系數(shù)的均值和方差對(duì)污染羽的二維空間分布的影響。上述研究涉及的含水層雖然有空間變異性,但相對(duì)還是比較均勻的(lnK的方差一般小于0.5,采用滲透系數(shù)符合對(duì)數(shù)正態(tài)分布的假設(shè)[15])。野外實(shí)際含水層的非均質(zhì)性是非常復(fù)雜的,非均質(zhì)性分區(qū)處理一般不可避免,此外,在砂層中夾有較大的亞砂土或亞黏土透鏡體也是很常見(jiàn)的,當(dāng)研究尺度不是很大時(shí),不能簡(jiǎn)單做統(tǒng)計(jì)平均處理,必須按非均質(zhì)性分區(qū)處理。本文研究含水層非均質(zhì)性分區(qū)對(duì)地下水蒙特卡羅模擬結(jié)果的影響,采用對(duì)應(yīng)充填和排列充填兩種方法在不同實(shí)現(xiàn)數(shù)下產(chǎn)生非均質(zhì)性分區(qū)滲透系數(shù)最終隨機(jī)場(chǎng),并對(duì)兩種方法蒙特卡羅模擬的結(jié)果進(jìn)行比較。本文研究不考慮研究區(qū)剖分單元的密度對(duì)模擬結(jié)果的影響,因此,所有模擬計(jì)算對(duì)研究區(qū)均采用同一種剖分密度。

1 研究區(qū)概況

研究區(qū)是在Borden試驗(yàn)場(chǎng)局部含水層的基礎(chǔ)上進(jìn)行改造得到的,為了體現(xiàn)含水層的非均質(zhì)性,在原Borden試驗(yàn)場(chǎng)中加入了一個(gè)矩形高透水性小區(qū)域(滲透系數(shù)為K2,其他區(qū)域滲透系數(shù)為K1),以模擬分析三維的非均質(zhì)各向異性承壓含水層中的溶質(zhì)運(yùn)移問(wèn)題。含水層的平面圖如圖1所示,整個(gè)研究區(qū)東西向距離為19 m,南北向距離為12 m,含水層厚為1.8 m。其中有一個(gè)相對(duì)高透水性的小區(qū),其平面上東西向長(zhǎng)為4.5 m,南北向?qū)挒?.5 m,垂向上貫穿整個(gè)含水層厚度。

圖1 研究區(qū)平面示意圖(單位:m)

對(duì)于水文地質(zhì)參數(shù)的確定,除了相對(duì)高透水性小區(qū)的滲透系數(shù)K2外,其他水文地質(zhì)參數(shù)均采用Borden試驗(yàn)場(chǎng)局部含水層的數(shù)據(jù)[16-17]:含水層在水平方向的相關(guān)長(zhǎng)度為2.8 m,垂向的相關(guān)長(zhǎng)度為0.12 m;含水層孔隙度為0.1,儲(chǔ)水率為0.25 mL/ m。本文研究基于滲透系數(shù)服從對(duì)數(shù)正態(tài)分布的假設(shè),在Borden含水層資料的基礎(chǔ)上突出大小區(qū)的透水性不同,在結(jié)果合理的前提下,假設(shè)lnK2的均值為4.605 170,方差為2.302 585,lnK1的均值為1.609 438,方差為0.3 147 165(即原Borden試驗(yàn)場(chǎng)的數(shù)值)。

含水層的左右兩邊邊界設(shè)為定水頭邊界,水頭分別為20.0 m和19.8 m，水流方向自左向右。總的模擬時(shí)間為9 d。假設(shè)圖1中Q點(diǎn)為一口不完整注水井,井深達(dá)到地面以下0.5 m,注水速率隨時(shí)間變化,在前3天為0.5 L/d,第4～6天為0.8 L/d,第7～9天為0.7 L/d。注入的水中,溶質(zhì)的質(zhì)量濃度為100 g/L。由于主要分析含水層的非均質(zhì)性分區(qū)對(duì)地下水蒙特卡羅模擬的影響,故僅考慮滲透系數(shù)的空間變異性,不考慮孔隙度的空間變異性。

2 蒙特卡羅模擬

2.1 滲透系數(shù)場(chǎng)的產(chǎn)生

采用Robin等[18]提出的直接傅立葉變換方法來(lái)生成單一的滲透系數(shù)隨機(jī)場(chǎng),對(duì)于本文的高度非均質(zhì)性研究區(qū)(即滲透系數(shù)變化大)必須進(jìn)行非均質(zhì)分區(qū),但如何產(chǎn)生最終的滲透系數(shù)隨機(jī)場(chǎng)目前還沒(méi)有確定的方法,本文分別采用兩種可能的方法得到最終滲透系數(shù)隨機(jī)場(chǎng)來(lái)進(jìn)行比較。

首先將研究區(qū)按非均質(zhì)性進(jìn)行分區(qū),外圍的大矩形區(qū)定為Ⅰ區(qū),內(nèi)部的小矩形區(qū)定為Ⅱ區(qū);然后分別按前述數(shù)據(jù)(lnK1的均值為1.609 438,方差為0.3 147 165;lnK2的均值為4.605 170,方差為2.302 585)產(chǎn)生Ⅰ和Ⅱ區(qū)滲透系數(shù)隨機(jī)場(chǎng);最后采用兩種可能的拼接方法得到最終的研究區(qū)滲透系數(shù)隨機(jī)場(chǎng)。

將滲透系數(shù)隨機(jī)場(chǎng)Ⅱ充填進(jìn)滲透系數(shù)隨機(jī)場(chǎng)Ⅰ是得到最終研究區(qū)滲透系數(shù)隨機(jī)場(chǎng)的基礎(chǔ),而本文采用的兩種拼接方法的區(qū)別也就在于充填方法的不同,下面以3次實(shí)現(xiàn)數(shù)為例說(shuō)明這兩種拼接方法。

假設(shè)已經(jīng)分別得到了隨機(jī)場(chǎng)Ⅰ和Ⅱ的3個(gè)不同實(shí)現(xiàn)(K1的隨機(jī)場(chǎng)1、K1的隨機(jī)場(chǎng)2、K1的隨機(jī)場(chǎng)3和K2的隨機(jī)場(chǎng)1、K2的隨機(jī)場(chǎng)2、K2的隨機(jī)場(chǎng)3),然后通過(guò)以下兩種充填拼接方法得到最終的研究區(qū)滲透系數(shù)K的隨機(jī)場(chǎng):

a. 對(duì)應(yīng)充填。將K2的3個(gè)隨機(jī)場(chǎng)按照實(shí)現(xiàn)數(shù)一一對(duì)應(yīng)充填到K1隨機(jī)場(chǎng)中的相應(yīng)位置,即

由此得到3個(gè)最終的研究區(qū)滲透系數(shù)K的隨機(jī)場(chǎng)。

b. 排列充填。將K2的3個(gè)隨機(jī)場(chǎng)按照實(shí)現(xiàn)數(shù)用排列的方法充填到K1的隨機(jī)場(chǎng)中,即

由此得到9個(gè)最終的研究區(qū)滲透系數(shù)K的隨機(jī)場(chǎng)。

2.2 蒙特卡羅模擬方法和步驟

蒙特卡羅方法是一種基于計(jì)算機(jī)模擬的統(tǒng)計(jì)試驗(yàn)方法,首先利用偽隨機(jī)數(shù)生成足夠多組服從給定分布規(guī)律的隨機(jī)變量(如滲透系數(shù)),然后對(duì)每組輸入變量分別進(jìn)行數(shù)值模擬。本文研究含水層非均質(zhì)性對(duì)蒙特卡羅模擬結(jié)果的影響大致分為以下5個(gè)步驟:①利用直接傅立葉變換方法分別在Ⅰ區(qū)和Ⅱ區(qū)空間離散的網(wǎng)格節(jié)點(diǎn)上隨機(jī)產(chǎn)生符合給定分布特征的滲透系數(shù)場(chǎng),即完成滲透系數(shù)場(chǎng)的1次實(shí)現(xiàn);②用兩種不同充填拼接方法得到最終的研究區(qū)滲透系數(shù)場(chǎng)的1次實(shí)現(xiàn);③根據(jù)生成的最終滲透系數(shù)場(chǎng),運(yùn)行對(duì)應(yīng)條件下的地下水流模型和溶質(zhì)運(yùn)移模型,以獲得整個(gè)研究區(qū)的地下水流場(chǎng)和溶質(zhì)的濃度分布場(chǎng);④由確定的實(shí)現(xiàn)次數(shù)重復(fù)前3步;⑤最后對(duì)所有實(shí)現(xiàn)的結(jié)果進(jìn)行統(tǒng)計(jì)。

將研究區(qū)離散為25行、39列、36層,其中小矩形區(qū)占據(jù)了第7～18行、第10～19列,垂向上為36層。水平面上沿x和y方向網(wǎng)格的間距均為0.5 m,垂向上網(wǎng)格間距為0.05 m。使用地下水流模擬模型MODFLOW[19]模擬地下水的流動(dòng),得到地下水流場(chǎng),然后利用溶質(zhì)運(yùn)移模型MT3DMS[20]追蹤溶質(zhì)質(zhì)點(diǎn)隨水流的運(yùn)動(dòng),最終得到研究區(qū)的濃度分布場(chǎng)。為了充分比較兩種充填方法的結(jié)果,在采用蒙特卡羅方法進(jìn)行隨機(jī)模擬時(shí),對(duì)應(yīng)充填分別嘗試了20次、50次、100次、150次和200次實(shí)現(xiàn),相應(yīng)的排列充填實(shí)現(xiàn)數(shù)分別為400次、2 500次、10 000次、22 500次和40 000次;另外為了更充分地比較模擬結(jié)果,還用對(duì)應(yīng)充填的方法計(jì)算了400次、2 500次、10 000次、22 500次和40 000次實(shí)現(xiàn)的結(jié)果。

3 模擬結(jié)果及分析

3.1 水頭場(chǎng)的比較

通過(guò)地下水流模型MODFLOW求解可以得到第3天、第6天和第9天的地下水流場(chǎng)。圖2給出了150次相關(guān)實(shí)現(xiàn)數(shù)第9天的水頭均值等值線不同方法模擬結(jié)果,可以看出兩種方法以及不同實(shí)現(xiàn)數(shù)間差別不明顯。這是由于本次模擬中只有單口井,注水量很小,且模擬時(shí)間短(只有9 d),所以模擬過(guò)程中研究區(qū)水頭的變化都不大。

圖2 150次相關(guān)實(shí)現(xiàn)數(shù)第9天的水頭均值等值線對(duì)比(單位:m)

此外,從圖2還可以看出,中間相對(duì)高透水性區(qū)域水頭均值等值線排列疏松,即水力坡度下降較慢;兩邊水頭均值等值線相對(duì)較密,水力坡度下降快,即為相對(duì)滲透系數(shù)小的區(qū)域,因此模擬得到的地下水流場(chǎng)是合理的。

3.2 水頭方差的比較

圖3 150次相關(guān)實(shí)現(xiàn)數(shù)第9天的水頭方差等值線對(duì)比(單位:m2)

用蒙特卡羅方法計(jì)算出的水頭和濃度都是存在方差的。由圖3可以看出兩種處理方法得到的水頭方差的差別主要表現(xiàn)在中間相對(duì)高透水性小區(qū)域中,而且較大的方差也是出現(xiàn)在小區(qū)域中,最大方差值出現(xiàn)在井點(diǎn)附近。總體而言,排列充填的方差范圍要比對(duì)應(yīng)充填的方差范圍小,而且隨實(shí)現(xiàn)數(shù)的增加,方差范圍逐漸縮小。從圖3可以看出,當(dāng)對(duì)應(yīng)充填的實(shí)現(xiàn)數(shù)增加到400次以后,所得同樣大小的方差范圍比相同實(shí)現(xiàn)數(shù)的排列充填的范圍小,也就是說(shuō)所得水頭均值變幅更小,即更準(zhǔn)確。濃度方差亦有相同規(guī)律,具體見(jiàn)下文分析。

3.3 質(zhì)量濃度均值的比較

使用特征線法追蹤溶質(zhì)質(zhì)點(diǎn)隨水流的運(yùn)動(dòng),可以獲得研究區(qū)溶質(zhì)在不同模擬時(shí)間的質(zhì)量濃度分布場(chǎng),為簡(jiǎn)便起見(jiàn),本文僅對(duì)最頂層的質(zhì)量濃度均值等值線進(jìn)行比較。

為了分析計(jì)算結(jié)果的可靠性,仍按傳統(tǒng)方法將同一參數(shù)分區(qū)里的滲透系數(shù)取為同一數(shù)值,即將滲透系數(shù)按均值處理進(jìn)行一次地下水流與溶質(zhì)運(yùn)移的模擬(以下簡(jiǎn)稱為按均值處理),然后與蒙特卡羅方法進(jìn)行比較。根據(jù)蒙特卡羅模擬所得的均值和方差結(jié)果可以繪制出蒙特卡羅計(jì)算結(jié)果的波動(dòng)范圍,這一波動(dòng)范圍代表了各種可能性的計(jì)算結(jié)果,也包括按均值處理的情況;因此,如果蒙特卡羅的計(jì)算結(jié)果是合理的,那么按均值處理的結(jié)果應(yīng)包含在蒙特卡羅計(jì)算結(jié)果的波動(dòng)范圍之中。將兩種方法不同實(shí)現(xiàn)數(shù)的蒙特卡羅計(jì)算結(jié)果波動(dòng)范圍與按均值處理的結(jié)果進(jìn)行對(duì)比分析表明,只有當(dāng)對(duì)應(yīng)充填的絕對(duì)實(shí)現(xiàn)數(shù)取到2 500次以上,排列充填的絕對(duì)實(shí)現(xiàn)數(shù)取到10 000次以上時(shí),按均值處理的結(jié)果才會(huì)完全落在蒙特卡羅計(jì)算結(jié)果的波動(dòng)范圍內(nèi)(圖4)。由此可見(jiàn),對(duì)于這種具有強(qiáng)烈非均質(zhì)性的含水層進(jìn)行蒙特卡羅模擬時(shí),以往經(jīng)常采用的幾百次的實(shí)現(xiàn)數(shù)是遠(yuǎn)遠(yuǎn)不夠的,對(duì)應(yīng)充填實(shí)現(xiàn)數(shù)要取到2 500次以上,而排列充填實(shí)現(xiàn)數(shù)要取到10 000次以上時(shí)所得的結(jié)果才是合理的。

圖4 蒙特卡羅模擬與按均值處理運(yùn)移5 d的質(zhì)量濃度均值對(duì)比(單位:g/L)

另外對(duì)比蒙特卡羅模擬的結(jié)果和按均值處理的結(jié)果,可以發(fā)現(xiàn)按均值處理的結(jié)果質(zhì)量濃度均值等值線主要為長(zhǎng)條形,質(zhì)點(diǎn)在縱向上隨水流運(yùn)移明顯,而蒙特卡羅模擬得到的質(zhì)量濃度均值等值線則在橫向上亦有明顯的分布,且在縱向上的擴(kuò)散要比按均值處理慢,由于水流方向是沿縱向自左向右,可見(jiàn)蒙特卡羅方法的計(jì)算結(jié)果既有縱向上因?qū)α饕鸬娜苜|(zhì)擴(kuò)散,又較好地體現(xiàn)了溶質(zhì)質(zhì)點(diǎn)在橫向上的彌散現(xiàn)象,蒙特卡羅方法能更好地體現(xiàn)含水層介質(zhì)的空間變異性。

比較不同實(shí)現(xiàn)數(shù)不同運(yùn)移時(shí)間下兩種處理方法的質(zhì)量濃度均值等值線(圖5),可以看出對(duì)應(yīng)充填方法和排列充填方法得到的質(zhì)量濃度場(chǎng)存在著明顯的差別,且隨著運(yùn)移時(shí)間的增長(zhǎng),區(qū)別也愈加明顯,但是隨實(shí)現(xiàn)數(shù)的增加,二者區(qū)別則相對(duì)減小。排列充填在實(shí)現(xiàn)數(shù)很少(100次以內(nèi))的情況下得出的結(jié)果顯然是不合理的,而當(dāng)對(duì)應(yīng)充填的實(shí)現(xiàn)數(shù)增加至與排列充填的絕對(duì)實(shí)現(xiàn)數(shù)相同時(shí),總體上溶質(zhì)在縱向上的擴(kuò)散排列充填要比對(duì)應(yīng)充填快。由圖5可以看出,在實(shí)現(xiàn)數(shù)低于10 000次的情況下,兩種方法得到的質(zhì)量濃度均值等值線有差別,主要表現(xiàn)在縱向上的擴(kuò)散快慢,而當(dāng)實(shí)現(xiàn)數(shù)達(dá)到22 500次以上時(shí),兩種方法所得結(jié)果就無(wú)明顯差別。

圖5 運(yùn)移9 d質(zhì)量濃度均值等值線對(duì)比(單位:g/L)

當(dāng)對(duì)應(yīng)充填的實(shí)現(xiàn)數(shù)為20次、50次、100次、150次和200次時(shí),排列充填的絕對(duì)實(shí)現(xiàn)數(shù)為400次、2 500次、10 000次、22 500次和40 000次,此時(shí)兩種方法所產(chǎn)生的基礎(chǔ)K1和K2的隨機(jī)場(chǎng)是相同的,但是對(duì)于排列充填而言等于是增加了蒙特卡羅模擬實(shí)現(xiàn)的次數(shù)。根據(jù)隨機(jī)理論,增加實(shí)現(xiàn)數(shù)最終結(jié)果的波動(dòng)會(huì)變小,排列充填得到的結(jié)果會(huì)更準(zhǔn)確。當(dāng)對(duì)應(yīng)充填的實(shí)現(xiàn)數(shù)增加到2 500次之后,所得的結(jié)果比相同絕對(duì)實(shí)現(xiàn)數(shù)下排列充填得到的結(jié)果更準(zhǔn)確,因?yàn)榇藭r(shí)排列充填的絕對(duì)實(shí)現(xiàn)數(shù)雖然和對(duì)應(yīng)充填相同,但是其并未增加基礎(chǔ)K1和K2的隨機(jī)場(chǎng)數(shù)目,而對(duì)應(yīng)充填在增加了實(shí)現(xiàn)數(shù)的同時(shí)還增加了產(chǎn)生的K1和K2的隨機(jī)場(chǎng)數(shù)目,這樣得到的最終K的隨機(jī)場(chǎng)更為符合實(shí)際情況,如表1所示。隨實(shí)現(xiàn)數(shù)的增加,所得到的K1和K2隨機(jī)場(chǎng)的均值和方差越接近給定的均值和方差,誤差越小,得到的最終滲透系數(shù)場(chǎng)也就更為準(zhǔn)確,因而所得到的質(zhì)量濃度結(jié)果亦更為準(zhǔn)確(前述水頭方差隨實(shí)現(xiàn)數(shù)增加而變小原因同此)。需要注意的是當(dāng)絕對(duì)實(shí)現(xiàn)數(shù)達(dá)到22 500次以上時(shí),兩種方法所得結(jié)果就無(wú)明顯差別。

表1 不同實(shí)現(xiàn)數(shù)時(shí)lnK1和lnK2的均值與方差

3.4 質(zhì)量濃度方差的比較

模擬結(jié)果表明,質(zhì)量濃度方差的分布規(guī)律與質(zhì)量濃度均值的分布規(guī)律相同。當(dāng)對(duì)應(yīng)充填的實(shí)現(xiàn)數(shù)為20次、50次和100次時(shí),質(zhì)量濃度方差對(duì)應(yīng)于質(zhì)量濃度值也存在一個(gè)狹長(zhǎng)帶,隨實(shí)現(xiàn)數(shù)的增加,該狹長(zhǎng)帶也就逐漸縮短并消失。質(zhì)量濃度方差的最大值出現(xiàn)在井點(diǎn)附近,并且隨著對(duì)應(yīng)充填實(shí)現(xiàn)數(shù)的增加,質(zhì)量濃度方差范圍明顯變小,質(zhì)量濃度方差最大值也隨實(shí)現(xiàn)數(shù)的增加而減小。

兩種方法的模擬得到的質(zhì)量濃度方差分布差別與均值差別的規(guī)律類(lèi)似,隨實(shí)現(xiàn)數(shù)的增加兩者的區(qū)別減小,當(dāng)實(shí)現(xiàn)數(shù)達(dá)到22 500以上時(shí),兩者的差別就很小。需要指出的是,對(duì)應(yīng)充填所得的最大方差值要小于排列充填所得的最大方差值。

表2為200次相關(guān)實(shí)現(xiàn)數(shù)不同方法模擬得到的質(zhì)量濃度方差最大值比較,在沒(méi)有增加對(duì)應(yīng)充填實(shí)現(xiàn)次數(shù)時(shí),排列充填所得質(zhì)量濃度方差最大值是小于對(duì)應(yīng)充填的,但是當(dāng)對(duì)應(yīng)充填的實(shí)現(xiàn)次數(shù)增加到與排列充填相同后,質(zhì)量濃度方差值就小于排列充填了,方差值變小,則所得質(zhì)量濃度值變幅就變小,結(jié)果也就更符合實(shí)際。

表2 200次相關(guān)實(shí)現(xiàn)數(shù)不同運(yùn)移時(shí)間質(zhì)量濃度方差最大值

4 結(jié) 論

對(duì)產(chǎn)生的符合研究區(qū)統(tǒng)計(jì)特性的隨機(jī)場(chǎng)多次實(shí)現(xiàn),進(jìn)行蒙特卡羅模擬可以充分考慮含水層水力參數(shù)的結(jié)構(gòu)性和隨機(jī)性,能較真實(shí)地刻畫(huà)含水層中溶質(zhì)運(yùn)移的過(guò)程,但是對(duì)含水層非均質(zhì)性的不同處理方法對(duì)蒙特卡羅模擬所得到結(jié)果是有影響的。對(duì)于本文研究的具有強(qiáng)烈非均質(zhì)性的含水層,對(duì)應(yīng)充填的蒙特卡羅模擬絕對(duì)實(shí)現(xiàn)數(shù)要取至2 500次,排列充填的蒙特卡羅模擬絕對(duì)實(shí)現(xiàn)數(shù)要取至10 000次,才能得到合理的計(jì)算結(jié)果。

[ 1 ] 薛禹群,謝春紅.水文地質(zhì)學(xué)的數(shù)值法[M].北京:煤炭工業(yè)出版社,1980.

[ 2 ] 薛禹群,謝春紅,吳吉春.水文地質(zhì)數(shù)值法存在的問(wèn)題及其對(duì)策[J].地球科學(xué)進(jìn)展,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 ] 薛禹群,吳吉春.地下水?dāng)?shù)值模擬在我國(guó):回顧與展望[J].水文地質(zhì)工程地質(zhì),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ì)的中國(guó)地下水模擬問(wèn)題[J].水文地質(zhì)工程地質(zhì),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 ] 楊金忠,蔡樹(shù)英,黃冠華,等.多孔介質(zhì)中水分及溶質(zhì)運(yùn)移的隨即理論[M].北京:科學(xué)出版社,2000.

[ 6 ] 楊金忠,蔡樹(shù)英,葉自桐.區(qū)域地下水溶質(zhì)運(yùn)移隨機(jī)理論的研究與進(jìn)展[J].水科學(xué)進(jìn)展,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] 陳彥,吳吉春.含水層滲透系數(shù)空間變異性對(duì)地下水?dāng)?shù)值模擬的影響[J].水科學(xué)進(jìn)展,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] 閻婷婷,吳劍鋒.滲透系數(shù)的空間變異性對(duì)污染物運(yùn)移的影響研究[J].水科學(xué)進(jìn)展,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，F(xiàn)AN 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

國(guó)家自然科學(xué)基金(40272106)

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

P641

A

1004-6933(2017)01-0046-06

2016-04-26 編輯：熊水斌)