李登波， 鄧 芳， 唐玉川， 張建民

(1.中國電建集團中南勘測設計研究院有限公司， 湖南 長沙 410014； 2.中國長江三峽集團有限公司， 北京 100038；3.四川大學 水力學與山區(qū)河流開發(fā)保護國家重點實驗室， 四川 成都 610065)

1 研究背景

地下水溶質運移受到多方面的影響，其中水文地質參數(shù)的不確定性是影響模型準確性的主要原因[1]。邊界條件或地質地貌改變都可能導致該研究區(qū)的水文地質參數(shù)發(fā)生變化，另外因含水介質形成條件的復雜性，導致大多數(shù)水文地質參數(shù)，比如滲透系數(shù)、孔隙度等具有空間變異性[2-4]。對于非均質介質中的溶質運移問題，水文地質參數(shù)的不確定性對于溶質運移有著重要的影響。Freeze的研究表明了滲透系數(shù)服從對數(shù)正態(tài)分布，后來這一結論也被一些學者的研究證明[5]。楊助澤[6]、李少龍等[7]和宋浩等[8]通過建立隨機數(shù)學模型得到了水文地質參數(shù)的空間變異性對水頭分布的影響。蔣立群等[9]和施小清等[10]通過隨機模擬法研究了水文地質參數(shù)的空間變異性對氯化物運移結果的影響。王建娥等[11]基于Choleshy分解進行了堆石材料的參數(shù)隨機場模擬，通過隨機有限元模型進行計算，得到了隨機有限元模型計算結果比常規(guī)有限元計算結果更接近實測值的結論。綜上所述，隨機模擬法可以較好地考慮到水文地質參數(shù)不確定性的影響，從而對地下水溶質運移有更為準確的描述。

本文分別考慮了水文地質參數(shù)數(shù)值不確定性和空間分布隨機性兩種因素對于氯化物運移的影響，定量比較了氯化物濃度和污染范圍的不確定性。在水文地質參數(shù)空間分布隨機性方面，基于3種空間插值技術對滲透系數(shù)的空間分布進行隨機模擬計算。通過不確定性分析可以掌握該地區(qū)氯化物濃度的變化范圍，在劃分不同濃度置信區(qū)間的基礎上進行風險分析，以期為地下水污染防治提供數(shù)據(jù)參考。

2 數(shù)據(jù)來源與研究方法

2.1 研究區(qū)概況

某垃圾填埋場位于非承壓含水層，含水層在該區(qū)域西南部最厚，在該區(qū)域的北端變薄至9.5 m，該含水層的厚度為15～30 m不等，為冰川河流砂含水層，含水層埋深較淺，底部由大面積的黏土質和砂質淤泥層構成相對隔水屏障。含水層主要由細粒、中粒和粗粒砂組成復雜地層，由于透鏡體的存在使得非承壓含水層局部是非均勻的，研究區(qū)內(nèi)的各水文地質參數(shù)存在一定的變化范圍，表1為研究區(qū)水文地質參數(shù)變化范圍表，參數(shù)的變化在其均值的20%范圍內(nèi)。

表1 研究區(qū)水文地質參數(shù)變化范圍表

根據(jù)已有資料建立該地區(qū)的地下水數(shù)值模型，如圖1所示。其中網(wǎng)格剖分采用結構化網(wǎng)格，網(wǎng)格總體數(shù)量為119 240個，最大網(wǎng)格尺寸為280 m。含水層補給來源為降水入滲，排泄為側向流出。上、下游為給定水頭邊界條件，每個月的水位隨著季節(jié)變化但是變化不大，上游邊界平均水頭為222.2 m，下游邊界平均水頭為218.5 m，左、右邊界條件為零流量邊界條件。研究區(qū)垃圾滲濾液入滲區(qū)位置如圖1所示，假設在某時刻垃圾場滲濾液開始滲漏并影響含水層，將垃圾滲濾液濃度概化為一個連續(xù)變化函數(shù)，濃度變化為先增大后逐漸減小最終保持穩(wěn)定，其在第30 a時濃度達到最大值。

圖1 研究區(qū)模型建立及其邊界條件

2.2 數(shù)據(jù)來源

本文的研究區(qū)位于加拿大多倫多以北約100 km處的某軍事基地，該填埋場啟用于1940年，一直使用至1976年。在1974-1980年期間，該地區(qū)逐漸建立起大量立式壓力計用于監(jiān)測該地區(qū)地下水位值，文獻[12]給出了該地區(qū)1983年的實測水位場分布圖。另外在含水層范圍內(nèi)有80臺多級監(jiān)測儀，從820個采樣點中獲取了氯化物的空間濃度值，并繪制了垃圾填埋場運行第43 a(1983年)時的氯化物濃度實測分布圖。根據(jù)Macfarlane等[13]和Cherry等[14]的研究成果，得到了該地區(qū)的滲透系數(shù)和孔隙度等水文參數(shù)的實測值，Sudicky[15]通過研究得到了該地區(qū)的彌散系數(shù)的實測值，本文基于這些文獻的實測資料值，對該地區(qū)氯化物的運移情況進行了不確定性分析。

2.3 研究方法

2.3.1 計算模型

(1)地下水流運動模型?？紤]到該研究區(qū)為非均質各向異性三維含水層，其控制方程可以用如下公式表示：

(1)

式中：K為滲透系數(shù)，m/d；μs為貯水率，m-1；H(x,y,z)為邊界Г1上的水頭，m；φ(x,y,z)為邊界Г1上的已知函數(shù)；H0(x,y,z)為t=0時邊界D上的初始水頭，m；Ω為研究區(qū)域； Г1、D為邊界條件。

(2)地下水對流-彌散方程。 該垃圾填埋場氯化物滲濾液的運移，其控制方程可以用如下公式表示：

(2)

式中：Dxx為縱向彌散系數(shù)，m2/s；Dyy、Dzz為橫向彌散系數(shù)，m2/s；ux、uy、uz為實際平均流速在坐標軸的3個分量，m/s；c(x,y,z,t)為邊界Г1上的濃度，kg/m3；φ(x,y,z,t)為邊界Г1上的已知函數(shù)；c0(x,y,z)為t=0時邊界D上的初始濃度，kg/m3。

隨機模型即在上述常規(guī)數(shù)值模型的基礎上考慮了參數(shù)的數(shù)值不確定性和空間分布隨機性，將得到的不同參數(shù)組合或者滿足一定空間分布狀態(tài)的輸入?yún)?shù)代入到模型中進行多次蒙特卡洛計算，最后進行計算結果的統(tǒng)計規(guī)律分析。

2.3.2 插值方法

(1)反距離權重插值。反距離權重插值(inverse distance weighted， IDW)使用一組采樣點的線性權重組合來確定像元值的方法。現(xiàn)有n個離散點Z(x1,y1,z1)，Z(x2,y2,z2)，…，Z(xn,yn,zn)，需要對插值點(x,y,z)進行插值預測[16]：

(3)

式中：F(x,y,z)為插值點(x,y,z)處的插值函數(shù)；wi為各個散點的權值函數(shù)；Qi為插值點的梯度超平面。

(2)克里金插值?？死锝鸱?Kriging)最早是由法國著名統(tǒng)計學家Matheron和南非礦山地質工程師Krige提出的，它是一種數(shù)學內(nèi)插及外推方法[17]?？死锝鹉Ｐ椭凶宰兞颗c響應值之間的關系可以用下式表示[18-19]：

(4)

式中：Z*(x)為x處克里金法響應值；λi為滿足無偏和估計方差最小兩個條件的待定權重系數(shù)；Z(xi)為第i個位置處的實際觀測值。

(3)自然鄰域插值。自然鄰域插值法(natural neighbor， NN)可找到距離查詢點最近的輸入樣本子集，并基于區(qū)域大小按照比例對這些樣本根據(jù)權重進行插值。當求某一位置點的預測值時，按下式進行計算[18]：

(5)

式中：f(x)為x處的預測值；Ni(x)為x點關于其自然鄰域節(jié)點的權重；fi為已知點函數(shù)值。

2.3.3 溶質運移結果的不確定性評價 為了定量描述氯化物運移結果的不確定程度，引入平均偏差值AEAE和平均方差值AESD兩個參數(shù)共同來描述計算結果的不確定性程度[20-21]。AEAE用于描述模擬值較實測值的偏差程度；AESD用于說明計算結果距離模擬值均值的波動程度，其表達式如下：

(6)

(7)

3 結果與分析

3.1 水文地質參數(shù)數(shù)值不確定性對模擬結果的影響及分析

通過參數(shù)靈敏度分析計算可知影響氯化物濃度最大的4個參數(shù)分別為水平滲透系數(shù)Kx、孔隙度n、垃圾填埋場處的入滲強度R和縱向離散度αL，因而取該4個水文地質參數(shù)進行數(shù)值不確定性分析。假設該4個水文地質參數(shù)在研究區(qū)內(nèi)相互獨立且均勻分布，每次隨機模擬時在每個參數(shù)取值范圍內(nèi)隨機抽取1個值并將其與其他3個參數(shù)進行組合，剩余水文地質參數(shù)計算時采用平均值，由此得到了1組輸入?yún)?shù)樣本值；將其代入隨機數(shù)值模型中進行計算得到1次蒙特卡洛隨機模擬結果。進行100次蒙特卡洛隨機模擬，統(tǒng)計100次隨機模擬結果值，得到了垃圾填埋場運行第43 a(1983年)4個監(jiān)測點W2、W3、W4和W6處(位置分布見圖1)的氯化物濃度統(tǒng)計參數(shù)值，表2為該4個監(jiān)測點濃度分布參數(shù)統(tǒng)計表。其中，最優(yōu)值是將濃度范圍等分為20個區(qū)間，取相對頻率最大的區(qū)間濃度中點值。

表2 4個監(jiān)測點濃度分布參數(shù)統(tǒng)計表(1983年)

由表2可知，統(tǒng)計第43 a監(jiān)測點W2的氯化物濃度分布區(qū)間為251.685～366.257 mg/L，氯化物濃度集中分布在316.000 mg/L，其平均濃度為315.847 mg/L；氯化物濃度分布偏度系數(shù)為-0.209，稍有負偏；峰度系數(shù)為-0.228，即氯化物濃度分布范圍較正態(tài)分布模型稍大，直方圖呈“矮胖”分布；對于監(jiān)測點W3，氯化物濃度范圍為111.859～196.938 mg/L，氯化物濃度集中分布在157.500 mg/L，其平均濃度為156.881 mg/L；其濃度分布偏度系數(shù)為0.068，基本呈對稱分布；峰度系數(shù)為-0.628，即氯化物濃度分布范圍較正態(tài)分布模型偏大，有“肥尾”現(xiàn)象；監(jiān)測點W4氯化物濃度范圍為63.352～118.379 mg/L，氯化物濃度集中分布在93.000 mg/L，其平均濃度為93.129 mg/L；氯化物濃度分布偏度系數(shù)為0.054，稍有“正偏”；峰度系數(shù)為-0.523，即氯化物濃度分布范圍較正態(tài)分布模型大，直方圖呈“矮胖”分布；監(jiān)測點W6氯化物濃度范圍為22.026～73.807 mg/L，其中氯化物濃度主要分布在48.500 mg/L，其平均濃度為47.049 mg/L；氯化物濃度分布偏度系數(shù)為0.251，呈“正偏”分布；峰度系數(shù)為-0.584，即氯化物濃度分布范圍較正態(tài)分布模型稍大，直方圖呈“矮胖”分布。

圖2為實測水位平均值與多次蒙特卡洛隨機模擬計算的平均值對比圖。由圖2可以看出，考慮參數(shù)數(shù)值不確定性的隨機模型水頭分布與實測水頭平均值符合較好，具有良好的適用性。

圖2 研究區(qū)實測水位與隨機模擬水位分布(單位：m)

繪制了數(shù)值不確定性情況下垃圾填埋場運行第43 a多次隨機模擬平均氯化物濃度分布圖，圖3為圖1中垃圾填埋場A-A斷面的剖切圖。由圖3可見，考慮數(shù)值不確定性計算得到的50 mg/L氯化物等濃度曲線在水平及垂直方向上的范圍與實測值范圍相同，400 mg/L等濃度曲線的范圍要大于實測值的濃度范圍。

圖3 研究區(qū)A-A斷面數(shù)值模擬的平均氯化物濃度與實測值對比

考慮到氯化物監(jiān)測點濃度的不確定性，本研究對監(jiān)測點濃度的不同置信水平濃度區(qū)間進行了劃分，表3為第43 a監(jiān)測點的百分位點濃度表。其中90%的置信水平濃度區(qū)間是監(jiān)測點5%～95%之間的濃度，其他置信水平濃度區(qū)間以此類推。由表3可見監(jiān)測點W290%的置信水平濃度區(qū)間為[265.211 mg/L，360.184 mg/L]，監(jiān)測點W390%的置信水平濃度區(qū)間為[123.320 mg/L，195.231 mg/L]，監(jiān)測點W490%的置信水平濃度區(qū)間為[72.580 mg/L，114.238 mg/L]，監(jiān)測點W690%的置信水平濃度區(qū)間為[24.633 mg/L，68.833 mg/L]，其他年份的置信水平濃度區(qū)間同理可得。隨著監(jiān)測點距離垃圾入滲區(qū)距離的增加，同等條件下氯化物濃度變化范圍逐漸減小，氯化物濃度的不確定性逐漸減小。

表3 研究區(qū)4個監(jiān)測點第43 a的氯化物百分位點濃度 mg/L

為了研究60 a計算時間段內(nèi)水文地質參數(shù)數(shù)值不確定性對垃圾填埋場氯化物濃度不確定性影響的變化情況，計算了整個研究區(qū)垃圾填埋場運行60 a內(nèi)即1940-2000年各時間節(jié)點的AEAE和AESD兩個參數(shù)值，其結果如表4所示。

表4 研究區(qū)各時間節(jié)點的AEAE和AESD計算結果 mg/L

表4中的數(shù)值結果顯示，研究區(qū)內(nèi)氯化物濃度變化幅度會隨著時間的改變而發(fā)生變化，氯化物濃度變化幅度越大，表明其不確定性也越大。從第10 a到第43 a，研究區(qū)內(nèi)AEAE和AESD不斷增加，表明其不確定性也逐漸增加，到第43 a達到最大值，AEAE和AESD分別為15.187和18.600；而后到第60 a逐漸減小，其不確定性也逐漸減小。就60 a整個計算時間段內(nèi)變化趨勢而言，數(shù)值不確定性對垃圾填埋場氯化物濃度的影響在第43 a最大。分析出現(xiàn)在第43 a的原因是污染源濃度在第30 a達到最大值，而由于地下水氯化物運移的滯后效應，所以氯化物濃度不確定性的最大值沒有與污染源濃度最大值出現(xiàn)在同一年。

由于水平滲透系數(shù)Kx、孔隙度n、垃圾填埋場處的入滲強度R、縱向離散度αL對氯化物滲濾液輸入量有影響，當該4個參數(shù)存在不確定性時，導致氯化物滲濾液入滲量存在不確定性，隨著氯化物濃度逐漸增大，其入滲量的變化范圍也逐漸增加，由于污染物的滯后效應，研究區(qū)氯化物濃度的AEAE和AESD兩個參數(shù)值逐漸增加并在第43 a達到極值，隨著氯化物入滲量變化范圍的減小，兩個參數(shù)值逐漸減小。

3.2 水文地質參數(shù)空間分布隨機性對模擬結果的影響及分析

比較計算結果可知，3種空間插值手段均值差異不大，因此本文以IDW法計算得到的平均濃度為例。圖4為考慮水平滲透系數(shù)空間分布隨機性條件下，隨機模擬平均值與實測氯化物濃度的對比。由圖4可見，考慮滲透系數(shù)空間分布隨機性模擬的濃度平均值較參數(shù)不確定性計算的平均值在極值濃度的分布范圍上有所減小，極值濃度的位置也與實測值位置更為吻合。

圖4 研究區(qū)A-A斷面數(shù)值模擬的平均氯化物濃度與實測值對比(考慮水平滲透系數(shù)空間分布隨機性)

仍選取W2、W3、W4和W6共4個監(jiān)測點，比較了水平滲透系數(shù)在空間分布隨機性的情況下對氯化物濃度的影響，結果見圖5。由圖5首先可以看到，相比于確定性模型，由于受到水平滲透系數(shù)的空間分布隨機性的影響，3種空間插值方法下模擬的各監(jiān)測點氯化物濃度不再是一個確定的值，而是一個變化的濃度范圍，該范圍包含所有可能解。

圖5 3種空間插值方法下4個監(jiān)測點氯化物濃度的變化范圍和標準差(第43 a)

不同監(jiān)測點的濃度變化范圍和標準差隨著距垃圾滲濾液入滲區(qū)距離的增大而減小，表明氯化物濃度的不確定性逐漸減小，雖然垃圾滲濾液入滲濃度為確定值，但受到水平滲透系數(shù)的空間分布隨機性的影響，研究區(qū)監(jiān)測點濃度為一個變化范圍，表明水文地質參數(shù)的空間隨機分布對于濃度影響較大。其次對比空間隨機性的3種空間插值方法，3種方法對各監(jiān)測點氯化物濃度均值的影響均較小，但是對于氯化物濃度的變化范圍和標準差，Kriging方法計算的結果最小，其次為NN方法，而IDW方法計算的結果最大。這說明3種隨機模擬方法中，計算結果的不確定性程度最小的是Kriging法，其次是NN法，而IDW方法不確定性程度最大，表明在數(shù)值模擬中采用不同的空間插值方法實現(xiàn)參數(shù)空間分布隨機性最終會使得計算結果的不確定性程度有所差別。

根據(jù)3種隨機模擬方法下氯化物濃度的不確定性分析，進行了監(jiān)測點濃度的置信區(qū)間劃分，圖6為第43 a研究區(qū)4個監(jiān)測點濃度置信區(qū)間劃分。其中90%的置信區(qū)間是監(jiān)測點百分位5%～95%之間的濃度，其他置信區(qū)間依此類推。就同一種空間插值方法，80%置信濃度區(qū)間的長度小于90%置信濃度區(qū)間的長度，隨著監(jiān)測點遠離垃圾滲濾液入滲點，它們之間的差距進一步減小，對比3種隨機模擬方法的90%置信濃度區(qū)間可知，Kriging法計算的置信濃度區(qū)間范圍最小，NN法次之，而IDW法置信濃度區(qū)間范圍最大。

圖6 研究區(qū)4個監(jiān)測點濃度的置信區(qū)間劃分(第43 a)

為了研究60 a計算時間段內(nèi)水文地質參數(shù)空間分布隨機性對垃圾填埋場氯化物濃度不確定性影響的變化情況，本研究針對AEAE和AESD兩個參數(shù)進行了分析計算。表5中列出了采用3種不同空間插值方法時，研究區(qū)AEAE和AESD系數(shù)的計算結果。由表5中的計算結果首先得到的結論是：水平滲透系數(shù)的空間隨機分布場對氯化物濃度計算結果有影響，且隨著時間的推移，其影響程度也隨之改變，呈現(xiàn)出先增加后減小的趨勢。

表5 研究區(qū)AEAE和AESD參數(shù)的3種不同空間插值法計算結果 mg/L

具體而言，3種插值方法在第10 a至第43 aAEAE和AESD值逐漸增加，第43 a達到最大，表明該時間節(jié)點氯化物濃度變化幅度最大，其不確定性達到最大；在43 a至60 a逐漸穩(wěn)步減小，表明該時段氯化物濃度變化幅度逐漸減小，其不確定性逐漸減小。對比3種空間插值方法，對于同一年份，Kriging法計算的AEAE和AESD兩個參數(shù)均較小，如在第43 a，Kriging法的AEAE和AESD分別為5.469和6.656 mg/L，其值為IDW法的45.03%和45.43%，NN法的59.93%和60.40%。說明Kriging方法得到的水平系數(shù)空間分布場，其計算的氯化物濃度不確定性最小，NN法次之，IDW法的不確定性最大。在AEAE和AESD最大值出現(xiàn)的時間上，3種空間插值方法的隨機模擬最大值均出現(xiàn)在第43 a。

將水文地質參數(shù)空間分布隨機性的3種空間插值方法和水文地質參數(shù)數(shù)值不確定性計算的AEAE和AESD進行比較，繪制了該4種方法的AEAE和AESD計算結果，如圖7所示。首先，氯化物濃度不確定性最大的排序為參數(shù)數(shù)值不確定性>空間分布隨機性。其次，對比60 a內(nèi)3種插值方法AEAE和AESD變化趨勢的空間分布隨機性，其排序為：反距離插值(IDW)法隨機模擬>自然鄰域(NN)法隨機模擬>克里格插值(Kriging)法隨機模擬。在AEAE和AESD最大值出現(xiàn)的時間上，在計算的60 a范圍內(nèi)參數(shù)數(shù)值不確定性和空間分布隨機性最大值均出現(xiàn)在第43 a，而污染源濃度在第30 a達到最大值，由于地下水氯化物運移的滯后效應，研究區(qū)內(nèi)垃圾滲濾液入滲量也滯后出現(xiàn)。無論是參數(shù)數(shù)值不確定性還是空間分布隨機性主要是由垃圾滲濾液入滲量決定最終計算結果，其不確定性程度大小與垃圾滲濾液入滲量成正比。

圖7 4種方法計算的AEAE和AESD結果對比

4 結 論

(1)不確定性模型隨機解的均值體現(xiàn)了地下水系統(tǒng)的結構性，它的極值和標準差又體現(xiàn)了地下水系統(tǒng)的隨機性。相比于確定性模型，不確定性模型更能體現(xiàn)出其在實際應用中的優(yōu)勢。

(2)地下水溶質運移受到參數(shù)數(shù)值不確定性和空間分布隨機性影響，在輸入?yún)?shù)±20%變化范圍內(nèi)，本案例中參數(shù)數(shù)值不確定性影響較空間分布隨機性影響程度大。

(3)在60 a的計算時間內(nèi)3種空間插值方法中，Kriging法得到的水平系數(shù)空間分布場，其計算的氯化物濃度不確定性最小，以43 a為例，AEAE和AESD參數(shù)值僅為同時間IDW法的45.03%和45.43%，NN法的59.93%和60.40%。

(4)在60 a的計算時間內(nèi)，參數(shù)數(shù)值不確定性和空間分布隨機性最大值都滯后出現(xiàn)(第43 a)，說明無論是參數(shù)數(shù)值不確定性還是空間分布隨機性均是通過影響垃圾滲濾液的入滲量來影響最終計算結果的。