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

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

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

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

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

1 研究方法

1.1 蒙特卡洛法

蒙特卡洛法是一種以概率和統計理論為基礎的具有廣泛實際應用價值的數學方法.它的基本思想是針對待求問題,根據物理現象本身的統計規律,或人為構造一合適的依賴隨機變量的概率模型,使某些隨機變量的統計值為該待求問題的解.通過進行大量統計實驗,得到該問題的解.

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

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

1.2 靈敏度分析法

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

Xk表示當參數αk變化時對輸出結果y的影響程度,即靈敏度系數.

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

由于不同參數的單位不同,導致靈敏度系數就無可比性,因此本文采用下列標準化無量綱形式計算靈敏度系數[21],即:

1.3 替代模型

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

克里格方法是常用的地質學統計方法之一[22],它從變量的相關性和變異性出發[23],在有限的研究區域內,對區域中變量取值進行誤差最小的優化估計.克里格方法利用協方差的變化來表達空間的變化,現被延伸為一種建立替代模型的方法,是一種黑箱模型[24].建立替代模型要完成以下幾個步驟:

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

式中:θk為待定系數,xki是第i個樣本的k維坐標.

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

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

R為n個采樣點的相關矩陣:

方差σ2的估計值為:

替代模型的建立可以通過求解上面的非線性無約束優化問題來實現.待定參數θk可以通過一個無約束優化問題求得:

2 應用研究

2.1 模型建立

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

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

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

表1 模型參數取值情況Table 1 Values of model parameters

式中:k代表不同組分的標號,k值取1,2時分別代表水,油;l代表不同相標號,l取1,2時分別代表水相、油相.φ代表孔隙度;Ck代表k組分的總濃度(體積分數);Ckl代表 k組分在 l相中的濃度(體積分數);ρk代表 k組分的密度(量綱為 ML-3);Sl代表相 l的飽和度;代表 k組分在 l相的彌散系數張量(量綱為 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相中的初始濃度(體積分數); Sl0(x, y, z)代表 l相的初始飽和度;Pl0(x, y, z)代表l相的初始壓力分布(量綱為ML-1T-2).邊界條件如下式所示:

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

2.2 靈敏度分析法篩選隨機參數

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

表2 各井各參數靈敏度計算結果Table 2 Sensitivity results of each parameter of each well

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

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

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

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

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

2.3 拉丁超立方抽樣

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

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

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

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

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

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

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

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

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

3 結果與討論

3.1 地下水污染情況分析

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

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

表4 各觀測井污染物濃度輸出值統計指標Table 4 Statistical indexes of pollutant concentration of observation well

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

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

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

3.2 地下水污染風險評價

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

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

3.2.2 全區地下水污染風險評價 為了對全區地下水 遭受不同程度污染的風險大小進行分區,以第 10層為例,計算1000個剖分網格處地下水輕度污染、中度污染以及重度污染的風險,繪制地下水遭受不同程度污染的風險分區圖.由圖 10可知,在地下水遭受輕度污染的條件下,當污染風險依次為30%,60%,90%時,風險分區的面積呈現逐漸減小的趨勢.由圖11、12可知,中度污染和重度污染的結論同上.比較圖10、11、12可知,在污染風險一定(30%,60%或90%)的條件下,當地下水分別遭受輕度、中度、重度污染時,風險分區的面積逐漸減小.

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

圖10 地下水輕度污染風險分區Fig.10 Contour map of groundwater mild pollution risk

圖11 地下水中度污染風險分區Fig.11 Contour map of groundwater moderate pollution risk

圖12 地下水重度污染風險分區Fig.12 Contour map of groundwater serious pollution risk

4 結論

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

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

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