含水層非均質性對地下水蒙特卡羅模擬結果的影響

王晶晶,樊尊榮

(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 編輯：熊水斌)