地下水DNAPLs污染多相流的隨機模擬及其不確定性分析

王 涵,盧文喜*,李久輝,常振波,侯澤宇 (1.吉林大學(xué)環(huán)境與資源學(xué)院,吉林 長春 130000；2.吉林大學(xué)地下水與資源環(huán)境教育部重點實驗室,吉林 長春 130000)

石油類有機污染物難溶于水,泄漏進入地下環(huán)境后對地下水造成污染,且難以去除,嚴(yán)重威脅人類健康[1].這類有機物進入地下環(huán)境后通常以非水相流體(Non-aqueous phase liquids,NAPLs)的形式存在[2].非水相流體根據(jù)密度的不同可分為兩類[3-4],密度小于水的稱為輕非水相流體(Non-aqueous phase liquids,簡稱 LNAPLs);密度大于水的稱為重非水相流體(Dense non-aqueous phase liquids,簡稱 DNAPLs).

DNAPLs具有高密度,低水溶性和高界面張力的特性[5],比LNAPLs更難修復(fù).為了合理高效地防治地下水DNAPLs污染,常常需要開展地下水DNAPLs污染風(fēng)險評價工作.地下水DNAPLs污染風(fēng)險評價的研究通常以多相流數(shù)值模擬為基礎(chǔ).然而,多相流數(shù)值模擬模型中水文地質(zhì)參數(shù)存在很大的不確定性[6].水文地質(zhì)參數(shù)的不確定性會導(dǎo)致多相流數(shù)值模型輸出結(jié)果的不確定性.因此,地下水污染風(fēng)險評價也具有不確定性[7].對多相流數(shù)值模擬模型中水文地質(zhì)參數(shù)的不確定性進行定量分析,可以使地下水DNAPLs污染風(fēng)險評價結(jié)果更加科學(xué)、全面.

近年來,不確定性分析在地下水溶質(zhì)運移方面發(fā)展迅速.Goodrich等[8]應(yīng)用蒙特卡洛方法進行不確定性分析,借助復(fù)雜數(shù)值模擬模型來模擬從假想污染場地到地下排水渠的溶質(zhì)運移過程; Hassan[9]應(yīng)用MCMC方法(Monte Carlo Markov Chain)方法評價了Amchitka Island (Alaska) Milrow試驗場地下水模型參數(shù)的不確定性;陳彥等[10]研究了含水層滲透系數(shù)空間變異性對地下水溶質(zhì)運移數(shù)值模擬結(jié)果造成的影響;吳吉春等[11]提出了將地下水?dāng)?shù)值模擬不確定性分析結(jié)果與風(fēng)險評估相結(jié)合,在此基礎(chǔ)上,常振波等[12],李久輝等[13]分析了參數(shù)不確定性對地下水溶質(zhì)運移數(shù)值模型輸出結(jié)果的影響,并從風(fēng)險評估的角度對不確定性分析結(jié)果進行闡釋.然而,這些研究均是針對簡單的二維單相流算例,目前尚未見有關(guān)于多相流隨機模擬及不確定性分析的報道.本文將不確定性分析與地下水DNAPLs污染多相流的隨機模擬相結(jié)合,運用靈敏度分析法篩選對模型輸出結(jié)果影響較大的參數(shù),作為模型中的隨機變量;為減少反復(fù)調(diào)用多相流模擬模型產(chǎn)生的計算負(fù)荷,運用克里格方法建立多相流模擬模型的替代模型,并利用替代模型完成蒙特卡洛隨機模擬.最后,對多相流模擬模型輸出結(jié)果進行統(tǒng)計分析并完成了地下水污染風(fēng)險評價.根據(jù)污染物濃度分布函數(shù)估算單井遭受污染的風(fēng)險.繪制地下水污染風(fēng)險分區(qū)圖,從而對全區(qū)地下水遭受不同程度污染的風(fēng)險大小進行分區(qū),為地下水污染防治提供合理依據(jù).

1 研究方法

1.1 蒙特卡洛法

蒙特卡洛法是一種以概率和統(tǒng)計理論為基礎(chǔ)的具有廣泛實際應(yīng)用價值的數(shù)學(xué)方法.它的基本思想是針對待求問題,根據(jù)物理現(xiàn)象本身的統(tǒng)計規(guī)律,或人為構(gòu)造一合適的依賴隨機變量的概率模型,使某些隨機變量的統(tǒng)計值為該待求問題的解.通過進行大量統(tǒng)計實驗,得到該問題的解.

蒙特卡洛法是分析復(fù)雜數(shù)值模型不確定性最有效的方法之一[14].它假設(shè)隨機變量概率分布函數(shù)和協(xié)方差函數(shù)已知,人工產(chǎn)生多組輸入變量,然后對每組輸入變量運行數(shù)值模型,獲得多組模型計算結(jié)果[15].本文運用蒙特卡洛方法進行地下水DNAPLs污染多相流的隨機模擬.

近年來,蒙特卡羅方法越來越多地應(yīng)用于地下水?dāng)?shù)值模擬不確定性分析方面.束龍倉等[16]利用蒙特卡羅方法,以濟寧市地下水潛水含水層為例,分析了給水度空間分布的隨機性對地下水庫庫容計算結(jié)果的影響.溫忠輝等[17]運用蒙特卡羅方法,在開采量不確定的條件下,對數(shù)值模擬結(jié)果進行可靠性分析,并對比計算了在不同的給定允許降深條件下,模擬結(jié)果的可靠性.束龍倉等[18]應(yīng)用蒙特卡羅法,研究水文地質(zhì)參數(shù)的不確定性對地下水補給量計算結(jié)果以及地表水與地下水交換量的影響.陸樂等[19]利用蒙特卡羅方法對多尺度非均質(zhì)含水層中的溶質(zhì)運移進行了模擬.綜上,蒙特卡羅方法在國內(nèi)地下水模擬不確定性分析方面得到了廣泛應(yīng)用.

1.2 靈敏度分析法

靈敏度可度量一種因子變化對另一因子的影響程度.靈敏度分析的方法主要有全局靈敏度分析法和局部靈敏度分析法[20].本文利用局部靈敏度分析的方法,評價單一參數(shù)變化對模型輸出結(jié)果造成的影響.計算時通過對參數(shù)求偏導(dǎo)數(shù)反映參數(shù)變化對模型輸出結(jié)果的影響程度,公式如下:

Xk表示當(dāng)參數(shù)αk變化時對輸出結(jié)果y的影響程度,即靈敏度系數(shù).

某一參數(shù)k的靈敏度系數(shù)值,可通過保證其他所有參數(shù)不變,使該參數(shù)k的值由αk變?yōu)棣羕+Δαk,相應(yīng)的因變量值由 yi(αk)變化為 yi(αk+Δαk),用下列計算獲得其計算值,即:

由于不同參數(shù)的單位不同,導(dǎo)致靈敏度系數(shù)就無可比性,因此本文采用下列標(biāo)準(zhǔn)化無量綱形式計算靈敏度系數(shù)[21],即:

1.3 替代模型

地下水DNAPLs污染多相流的隨機模擬需要反復(fù)多次調(diào)用多相流模擬模型.為減少上述計算負(fù)荷,應(yīng)用克里格方法,建立地下水DNAPLs污染多相流模擬模型的替代模型.在蒙特卡洛模擬過程中可以直接調(diào)用模擬模型的替代模型.

克里格方法是常用的地質(zhì)學(xué)統(tǒng)計方法之一[22],它從變量的相關(guān)性和變異性出發(fā)[23],在有限的研究區(qū)域內(nèi),對區(qū)域中變量取值進行誤差最小的優(yōu)化估計.克里格方法利用協(xié)方差的變化來表達空間的變化,現(xiàn)被延伸為一種建立替代模型的方法,是一種黑箱模型[24].建立替代模型要完成以下幾個步驟:

式中: R ( xi,xj)(i = 1,2， … m; j = 1,2， … m )代表采樣點xi和點 xj的關(guān)聯(lián)函數(shù).關(guān)聯(lián)函數(shù)模型包括很多種,本文采用應(yīng)用最為廣泛的高斯模型,模型形式如下:

式中:θk為待定系數(shù),xki是第i個樣本的k維坐標(biāo).

根據(jù)克里格模型,在預(yù)測點 x點處的響應(yīng)值 y(x)的預(yù)測估計值為

式中:f為基函數(shù),為方便起見,本文選擇常數(shù)型, f為一常數(shù)列向量;r( x)為點x與n個訓(xùn)練樣本采樣點(x1,x2,… xn)之間的相關(guān)向量;y為與n個采樣點對應(yīng)的響應(yīng)值,為n×1的向量;β為線性回歸部分的待定參數(shù),可以通過最優(yōu)線性無偏估計求得:

R為n個采樣點的相關(guān)矩陣:

方差σ2的估計值為:

替代模型的建立可以通過求解上面的非線性無約束優(yōu)化問題來實現(xiàn).待定參數(shù)θk可以通過一個無約束優(yōu)化問題求得:

2 應(yīng)用研究

2.1 模型建立

本文針對假想算例展開研究.模擬研究場地為140m×60m×20m的三維區(qū)域,計算目的層為松散巖類孔隙潛水含水層.含水層介質(zhì)為均質(zhì)各向同性介質(zhì).地下水流方向由左向右.左右邊界為定壓邊界(左邊界壓力為 102.972kPa,右邊界壓力為 101.325kPa),其余邊界為隔水邊界.在研究區(qū)內(nèi),有一個硝基苯泄漏點,泄漏量為 8m3/d,泄漏時間為 10d.硝基苯濃度為70%(體積分?jǐn)?shù)).模擬時間持續(xù) 300d,模型中各參數(shù)取值見表 1.潛水含水層污染物初始濃度為 0.研究區(qū)內(nèi)污染源與觀測井的平面位置見圖1.

在研究區(qū)水文地質(zhì)概念模型的基礎(chǔ)上建立地下水DANPLs污染多相流數(shù)值模擬模型.多相流數(shù)值模擬模型由偏微分方程和定解條件兩部分組成.偏微分方程以質(zhì)量守恒定律為理論基礎(chǔ)[25],公式如下:

圖1 研究區(qū)污染源與觀測井位置示意Fig.1 Location of the pollution source and observation wells in study area

表1 模型參數(shù)取值情況Table 1 Values of model parameters

式中:k代表不同組分的標(biāo)號,k值取1,2時分別代表水,油;l代表不同相標(biāo)號,l取1,2時分別代表水相、油相.φ代表孔隙度;Ck代表k組分的總濃度(體積分?jǐn)?shù));Ckl代表 k組分在 l相中的濃度(體積分?jǐn)?shù));ρk代表 k組分的密度(量綱為 ML-3);Sl代表相 l的飽和度;代表 k組分在 l相的彌散系數(shù)張量(量綱為 L2/T);Rk代表 k組分總的源匯項;→代表 l相的滲透流速(量綱為 L/T);krl代表相 l的相對滲透率;代表固有滲透率(量綱為 L2);ρl代表相 l的密度(量綱為 ML-3);g代表重力加速度(量綱為LT-2);z代表垂向距離(量綱為 L);Pl代表相 l的壓力[量綱為 M/(L?T2)].

定解條件包括初始條件和邊界條件兩部分.初始條件如下式所示:

式中:Ckl0(x, y, z)代表k組分在l相中的初始濃度(體積分?jǐn)?shù)); Sl0(x, y, z)代表 l相的初始飽和度;Pl0(x, y, z)代表l相的初始壓力分布(量綱為ML-1T-2).邊界條件如下式所示:

圖2 研究區(qū)剖分示意Fig.2 Subdivision schematic diagram of study area

2.2 靈敏度分析法篩選隨機參數(shù)

靈敏度分析可以定量評價參數(shù)不確定性對模擬模型輸出結(jié)果造成的影響.利用靈敏度分析方法篩選對模擬模型輸出結(jié)果影響較大的參數(shù),作為模型中的隨機變量,既減少計算負(fù)荷,又提高研究精度.

表2 各井各參數(shù)靈敏度計算結(jié)果Table 2 Sensitivity results of each parameter of each well

進行靈敏度分析時,將模型中各參數(shù)取波動范圍內(nèi)均值,見表3,輸出一組污染物濃度值.再分別將模型中各參數(shù)加減 20%,10%,同時保證其他參數(shù)不變,分別輸出污染物濃度值.利用公式(3)計算各井各參數(shù)靈敏度系數(shù)值,見表 2,并篩選出兩個靈敏度較大的參數(shù)作為模型中的隨機變量.

圖3 1號井參數(shù)靈敏度分析Fig.3 Sensitivity analysis of No.1 well

圖4 2號井參數(shù)靈敏度分析Fig.4 Sensitivity analysis of No.2 well

由圖3,圖4,圖5可知,對于1、2、3號井,靈敏度較大的兩個參數(shù)均為孔隙度和滲透率.將孔隙度和滲透率作為隨機模擬過程中的隨機變量,將其他參數(shù)作為常數(shù).

圖5 3號井參數(shù)靈敏度分析Fig.5 Sensitivity analysis of No.3 well

2.3 拉丁超立方抽樣

本文替代模型的輸入變量為隨機模擬過程中的隨機變量,即孔隙度和滲透率,輸出變量為污染物濃度.建立模擬模型的替代模型需要一定數(shù)量的訓(xùn)練樣本和檢驗樣本,為使樣本具有代表性,應(yīng)用拉丁超立方抽樣方法在孔隙度、滲透率的合理變化范圍內(nèi)進行抽樣.

拉丁超立方抽樣是一種基于分層采樣的抽樣方法,具有較高的采樣效率,能覆蓋到所有采樣區(qū)間.它對變量按照假定的概率密度等概率劃分,在每個區(qū)間內(nèi)抽取等數(shù)量的參數(shù)組成參數(shù)樣本.根據(jù)前人經(jīng)驗可知,孔隙度服從正態(tài)分布,滲透率服從對數(shù)正態(tài)分布[26].

運用 MATLAB軟件編寫拉丁超立方抽樣程序,通過兩次抽樣分別獲得40組訓(xùn)練樣本,和10組檢驗樣本.

表3 參數(shù)概率分布及取值情況Table 3 Values and probability distribution of model parameters

2.4 替代模型的建立及基于替代模型的隨機模擬

將50組樣本輸入地下水DNAPLs污染多相流模擬模型中,利用UTCHEM軟件進行求解,獲得50組輸入-輸出樣本數(shù)據(jù)集.利用訓(xùn)練樣本建立替代模型,并利用檢驗樣本檢驗替代模型對模擬模型的逼近精度.

由圖6可知,替代模型與模擬模型輸出結(jié)果的相對誤差均小于0.80%,誤差較小,故認(rèn)為替代模型精度較高,可以用來代替模擬模型.

運用蒙特卡洛法進行地下水DNAPLs污染的隨機模擬.首先,將孔隙度和滲透率作為隨機變量,利用拉丁超立方抽樣獲得1000組參數(shù)組合;然后,將1000組參數(shù)組合輸入替代模型,輸出1、2、3號井1000組污染物濃度數(shù)據(jù);最后,對1、2、3號井污染物濃度數(shù)據(jù)進行統(tǒng)計分析.

圖7 各井污染物濃度累計頻次Fig.7 Pollutant concentration histogram of each observation well

3 結(jié)果與討論

3.1 地下水污染情況分析

3.1.1 單井污染物濃度統(tǒng)計分析 對3口觀測井各1000組污染物濃度數(shù)據(jù)的各項統(tǒng)計指標(biāo)進行分析,繪制3口觀測井污染物濃度累計頻次直方圖(圖7).1號井污染物濃度在 20～40mg/L范圍內(nèi)的概率高達83.5%,可以預(yù)測 300d后 1號井污染物濃度在20～40mg/L范圍內(nèi)的可能性最高.同理可以預(yù)測,300d后 2、3號井污染物濃度 50～90mg/L和 10～30mg/L范圍內(nèi)的可能性最高.由表4可以看出,2號井污染物濃度數(shù)據(jù)的變異系數(shù)較小,分布較集中.3號井污染物濃度數(shù)據(jù)的變異系數(shù)較大,分布較分散.這說明,在同一時間條件下,污染物濃度的不確定性大小主要取決于觀測井和污染源的距離.距離污染源越近,污染質(zhì)運移的距離越短,受含水層各參數(shù)的影響越小,所以,污染物濃度的不確定性就越小.

在同一空間位置的條件下,以3號井為例,經(jīng)計算,第 100,200,300d污染物濃度數(shù)據(jù)的變異系數(shù)分別為0.12,0.25,0.31.這說明在相同位置處,污染物濃度不確定性的變化速率不與時間成正比,增速逐漸降低.其他位置同理.

表4 各觀測井污染物濃度輸出值統(tǒng)計指標(biāo)Table 4 Statistical indexes of pollutant concentration of observation well

圖8 研究區(qū)地下水污染程度分區(qū)Fig.8 Contour map of contaminant concentration in study area

3.1.2 全區(qū)地下水污染情況分析 為劃分全區(qū)地下水遭受不同程度污染的區(qū)域,以第 10層為例,計算每個剖分網(wǎng)格各 1000組污染物濃度的均值,繪制等值線圖(圖8).假設(shè)污染物濃度超過20mg/L為輕度污染區(qū)域,超過40mg/L為中度污染區(qū)域,超過60mg/L為重度污染區(qū)域.由圖8可知,300d后研究區(qū)輕度污染區(qū)域面積約為2715m2,占全區(qū)總面積的32.3%.嚴(yán)重污染區(qū)域面積約為825m2,占全區(qū)總面積的18.5%.

圖9 各井污染物濃度分布函數(shù)曲線Fig.9 Pollutant concentration distribution function of each observation well

3.2 地下水污染風(fēng)險評價

3.2.1 單井污染風(fēng)險評價 利用 MATLAB軟件將各井污染物濃度數(shù)據(jù)作為隨機變量,繪制污染物濃度分布函數(shù)曲線(圖9).假設(shè)污染物濃度超過20mg/L為輕度污染,超過40mg/L為中度污染,超過60mg/L為重度污染.由圖9可知1,2,3號井為輕度污染的風(fēng)險分別

為 92%,100%,40%,為中度污染的風(fēng)險分別為20%,99%,0%;為重度污染的風(fēng)險分別為 0%,51%,0%.由風(fēng)險評價結(jié)果可知,3號井水質(zhì)較好,該井附近適宜人類居住.

3.2.2 全區(qū)地下水污染風(fēng)險評價 為了對全區(qū)地下水 遭受不同程度污染的風(fēng)險大小進行分區(qū),以第 10層為例,計算1000個剖分網(wǎng)格處地下水輕度污染、中度污染以及重度污染的風(fēng)險,繪制地下水遭受不同程度污染的風(fēng)險分區(qū)圖.由圖 10可知,在地下水遭受輕度污染的條件下,當(dāng)污染風(fēng)險依次為30%,60%,90%時,風(fēng)險分區(qū)的面積呈現(xiàn)逐漸減小的趨勢.由圖11、12可知,中度污染和重度污染的結(jié)論同上.比較圖10、11、12可知,在污染風(fēng)險一定(30%,60%或90%)的條件下,當(dāng)?shù)叵滤謩e遭受輕度、中度、重度污染時,風(fēng)險分區(qū)的面積逐漸減小.

綜上,在污染程度相同的條件下,污染風(fēng)險越大,風(fēng)險分區(qū)的面積越小;在污染風(fēng)險相同的條件下,污染程度越重,風(fēng)險分區(qū)的面積越小.地下水污染風(fēng)險分區(qū)圖可以幫助決策者分析出地下水污染的高風(fēng)險區(qū),估計高風(fēng)險區(qū)的位置和面積,合理開展地下水污染風(fēng)險評價工作,為地下水污染防治提供更加科學(xué)可靠的依據(jù).相關(guān)人員也可以根據(jù)地下水污染風(fēng)險分區(qū)圖,選擇人類宜居地點.

圖10 地下水輕度污染風(fēng)險分區(qū)Fig.10 Contour map of groundwater mild pollution risk

圖11 地下水中度污染風(fēng)險分區(qū)Fig.11 Contour map of groundwater moderate pollution risk

圖12 地下水重度污染風(fēng)險分區(qū)Fig.12 Contour map of groundwater serious pollution risk

4 結(jié)論

4.1 靈敏度分析結(jié)果顯示:地下水DNAPLs污染多相流模擬模型中孔隙度和滲透率的靈敏度較大.利用靈敏度分析法可以篩選對模型輸出結(jié)果影響較大的參數(shù),減小不確定性分析的工作負(fù)荷.

4.2 本文運用克里格方法建立了地下水 DNAPLs污染多相流模擬模型的替代模型.既能大幅度減少不確定性分析過程中的計算負(fù)荷,還可以保持較高的精度.

4.3 本文分析了水文地質(zhì)參數(shù)不確定性對多相流數(shù)值模擬模型輸出結(jié)果的影響.利用蒙特卡洛方法進行地下水DNAPLs污染多相流的隨機模擬,并對隨機模擬結(jié)果進行統(tǒng)計分析.利用污染物濃度分布函數(shù)估算單井遭受污染的風(fēng)險;利用地下水污染風(fēng)險分區(qū)圖分析地下水污染的高風(fēng)險區(qū),為地下水污染的防治提供合理依據(jù).