陳夢佳，吳劍鋒，孫曉敏，林 錦，吳吉春

(1.南京大學(xué)地球科學(xué)與工程學(xué)院水科學(xué)系/表生地球化學(xué)教育部重點實驗室，江蘇 南京 210023；2.南京水利科學(xué)研究院，江蘇 南京 210029)

非水相液體(non-aqueous phase liquids, NAPLs)在工業(yè)生產(chǎn)中應(yīng)用廣泛，主要包括煤焦油、氯化溶劑等，是目前地下水污染的主要污染源之一。由于NAPLs在水中的溶解度相對較低，一旦發(fā)生泄漏，就會形成一個獨立自由相[1]。在飽和帶中，NAPLs常以不連續(xù)的形式分布在污染源區(qū)，緩慢地溶解于地下水并隨之流動擴散，逐漸形成大面積的污染羽[2-4]。

NAPLs運移模型可為制定污染修復(fù)方案提供關(guān)鍵技術(shù)支撐。目前，UTCHEM是一款廣泛用于地下水NAPLs運移的多相流數(shù)值模擬軟件。在復(fù)雜的地下水污染運移數(shù)值模擬過程中，場地面積、網(wǎng)格大小、運移持續(xù)時間等都對計算時間和模擬精度有顯著影響[5]。一般說來，為了盡可能節(jié)省計算成本，在地下水?dāng)?shù)值模擬中尤其涉及到蒙特卡洛隨機模擬時，模型的參數(shù)網(wǎng)格尺度一般大于實際測量的參數(shù)尺度[6-8]。為此，需要在保證模擬精度的情況下，將小尺度的參數(shù)轉(zhuǎn)換為大尺度的等效參數(shù)，亦即需要尺度提升。

隨著社會的發(fā)展，提出的問題變得越來越復(fù)雜，范圍更大、時間段更長、水文地質(zhì)條件更為復(fù)雜。對于這些問題的計算，傳統(tǒng)的地下水?dāng)?shù)值模擬需要巨大的計算成本來保證解的精度。研究表明，對參數(shù)進(jìn)行尺度提升能夠非常有效地減少計算時間[9-11]。由于滲透系數(shù)的確定在地下污染物的運移模擬中起著至關(guān)重要的作用，所以尺度提升的重點是將小尺度的滲透系數(shù)進(jìn)行轉(zhuǎn)換，以得到數(shù)值模型尺度上的等效參數(shù)，使得大尺度模型保持小尺度上的污染羽特征[12-13]。常用的尺度提升法有簡單平均法(simple averaging)、流管法(stream-tube)、簡單拉普拉斯法(simple Laplacian)等[14-16]。拉普拉斯-外殼法(Laplacian with skin)是對簡單拉普拉斯法的改進(jìn)，是一種基于求解流動方程的滲透系數(shù)尺度提升方法，而非依靠經(jīng)驗，因此能夠運用于各類場地。本文利用拉普拉斯-外殼法對滲透系數(shù)場進(jìn)行尺度提升并建立大尺度地下水典型NAPLs數(shù)值模型，與尺度提升前的小尺度模型及基于算術(shù)平均的尺度提升模擬結(jié)果進(jìn)行對比，評價拉普拉斯-外殼法的提升效果。

1 研究方法

與傳統(tǒng)的基于算術(shù)平均的尺度提升方法相比，簡單拉普拉斯尺度提升法基于求解流動方程，因此在理論與實際領(lǐng)域都有廣泛的應(yīng)用[17-20]。但是，簡單拉普拉斯尺度提升方法假設(shè)所得的等效滲透系數(shù)為對角張量，一旦當(dāng)小尺度的滲透系數(shù)產(chǎn)生的總流量不平行于參考軸，則通過這種尺度提升方法得到的模擬結(jié)果無法準(zhǔn)確描述小尺度模型中的運移行為[21-22]。因此，針對這一缺陷，Gómez-Hernández[23]假設(shè)等效滲透系數(shù)為全張量，提出了拉普拉斯-外殼法；李良平等[24]將拉普拉斯-外殼法應(yīng)用到綜合實例中，論證了拉普拉斯-外殼法相較于其他尺度提升方法的優(yōu)勢。

滲透系數(shù)尺度提升的目的是將小尺度的滲透系數(shù)轉(zhuǎn)換為大尺度的等效滲透系數(shù)，使其能運用于數(shù)值模型中，因此轉(zhuǎn)換需要滿2個條件：

(1)大尺度網(wǎng)格(block)中的水頭應(yīng)等于小尺度網(wǎng)格(cell)中水頭的平均：

(1)

式中：V——大尺度網(wǎng)格范圍；

nV——大尺度網(wǎng)格內(nèi)小尺度網(wǎng)格的數(shù)目；

hV(I,J)——大尺度網(wǎng)格內(nèi)水頭/m；

h(i,j)——小尺度網(wǎng)格內(nèi)水頭/m。

(2)大尺度網(wǎng)格中的流量應(yīng)等于小尺度網(wǎng)格中流量的平均：

(2)

式中：qV(I,J)——大尺度網(wǎng)格內(nèi)流量/(m3·s-1)；

q(i,j)——小尺度網(wǎng)格內(nèi)流量/(m3·s-1)。

拉普拉斯-外殼法的一個重要特點是使用外殼區(qū)域估計實際的大尺度網(wǎng)格的邊界條件，而不是拉普拉斯法中簡單的人為定義邊界條件，因此小尺度滲透系數(shù)場的區(qū)域大小要比實際含水層模型稍大，以便使模型邊界附近的網(wǎng)格也含有外殼。大尺度模型中外殼選取、網(wǎng)格內(nèi)水頭與流量定義見圖1。

對于小尺度網(wǎng)格的多套邊界條件，計算小尺度網(wǎng)格中的流動：

(3)

將等效滲透系數(shù)定義為：

(4)

運用計算得到大尺度網(wǎng)格中的參數(shù)，建立超定系統(tǒng)方程：

(5)

式中：Kxx，Kxy，Kyy——未知的等效滲透系數(shù)全張量的成分。

圖1 尺度提升示意圖Fig.1 Schematic diagram showing the scale-up mechanism

利用最小二乘法[25]解方程組，可獲得尺度提升后各計算網(wǎng)格的等效滲透系數(shù)，由此可通過UTCHEM代碼進(jìn)一步計算地下水中NAPLs污染羽隨時間變化的空間分布。

2 算例應(yīng)用

2.1 算例概述

研究對象為二維非均質(zhì)承壓含水層，含水層面積為15 085 m2，厚度為5 m，地下水流方向從左向右，采用序貫高斯模擬生成二維非均質(zhì)滲透率場[26-28]，滲透率場均值為10-10m2。含水層的有效孔隙度為0.35，水力梯度為0.001，縱向彌散度為1 m，橫向彌散度為0.1 m。模型中涉及的液相參數(shù)：水的密度為1.0 g/cm3，PCE的密度為1.63 g/cm3，水的黏滯性為0.001 Pa·s，PCE的黏滯性為0.000 89 Pa·s，PCE與水的界面張力為0.045 N/m，PCE在水中的溶解度為240.0 mg/L。由于某種人為原因，造成PCE泄漏，并穿過弱透水層進(jìn)入含水層。污染源分別有以下兩種情形：

Case 1：單一污染源(單源)，PCE泄漏點(圖2a)，泄漏量為1.0 m3/d。運移過程在時間上分為兩個階段：0～365 d的PCE泄露過程，366～730 d的PCE自然運移過程。

Case 2：雙重污染源(二源)，PCE泄漏點(圖2b)，泄漏量均為0.7 m3/d。運移過程與Case 1相同。

2.2 模型建立及尺度提升方案

根據(jù)參數(shù)先建立小尺度數(shù)值模型，橫向劃分為150個網(wǎng)格，縱向劃分為85個網(wǎng)格，垂向為一層，則小尺度離散為12 750個單元格，每個單元格大小為x=1.0 m，y=1.0 m，這與微水試驗測量的滲透系數(shù)尺度相似。平行于流動方向的邊界為隔水邊界，垂直于流動方向的邊界為定水頭邊界。

大尺度離散采用非均勻離散，尺度提升方法運用算術(shù)平均法和拉普拉斯-外殼法。拉普拉斯-外殼法采用的是網(wǎng)格中心尺度提升，x軸方向外殼的大小取小尺度的20網(wǎng)格，y軸方向取小尺度的10網(wǎng)格。

圖2 污染源模型的非均勻剖分Fig.2 Non-uniform block discretized model with polluted sources

對單一污染源污染點位置進(jìn)行非均勻剖分(圖2a)。算術(shù)平均尺度提升法將原始的滲透系數(shù)網(wǎng)格15 085(單元格大小11)提升到3 421，網(wǎng)格數(shù)從小尺度單元數(shù)12 750提升到大尺度單元數(shù)714。拉普拉斯-外殼法由于外殼不進(jìn)行計算，因此將滲透系數(shù)網(wǎng)格提升到2 617，單元數(shù)為442。

對雙重污染源兩個污染點位置分別進(jìn)行網(wǎng)格加密(圖2b)，算術(shù)平均法將原始的滲透系數(shù)網(wǎng)格15 085(單元格大小11)尺度提升到3 825，從小尺度單元數(shù)12 750提升到大尺度單元數(shù)950，拉普拉斯-外殼法將網(wǎng)格提升到3 021，單元數(shù)為630。

3 結(jié)果與分析

3.1 單一污染源

圖3為均質(zhì)滲透系數(shù)場下，單源泄漏的PCE運行了2 a(730 d)后，在不同尺度下UTCHEM計算出的PCE運移結(jié)果，表1為污染羽的空間矩計算結(jié)果。由表1可知，在均質(zhì)情況下，大尺度模型和小尺度模型得到污染羽的二階矩存在一定誤差，但其零階矩和一階矩誤差很小，表明大尺度模型基本能再現(xiàn)小尺度模型精細(xì)刻畫污染物的運移情況(圖3)，以此可作為非均質(zhì)情況的對照。

圖4為單源情況下，不同滲透系數(shù)場中，PCE運行了2 a(730 d)后，UTCHEM根據(jù)不同網(wǎng)格剖分情況計算出的PCE運移結(jié)果，表2為污染羽的空間矩分析結(jié)果。大尺度模型對于含水層中污染物質(zhì)量(零階矩)的估算準(zhǔn)確，相對誤差絕對值均小于2%，但隨著滲透系數(shù)對數(shù)方差的增大，誤差增大。拉普拉斯-外殼法較算術(shù)平均法更為準(zhǔn)確，且精確度受滲透系數(shù)對數(shù)方差的影響較小。對于質(zhì)心位置(一階矩)和污染羽在空間上的展布范圍(二階矩)，大尺度模型也能較好地刻畫小尺度模型中的情形。拉普拉斯-外殼法與算術(shù)平均法的結(jié)果對比，在垂直水流方向上，兩方法得到的效果接近；沿水流方向，拉普拉斯-外殼法能更好地刻畫PCE的運移情況，各滲透系數(shù)場下，不僅質(zhì)心位置更接近小尺度情況，而且污染物的最遠(yuǎn)運移距離也更接近。由于進(jìn)行含水層污染物治理時，多選擇在污染羽的上游和下游設(shè)井，因此，經(jīng)過拉普拉斯-外殼法尺度提升后得到的大尺度模型能更好地為含水層修復(fù)提供參考。

表1 單一污染源均質(zhì)滲透系數(shù)場中污染羽在空間上的分布特征

圖4 單一污染源非均質(zhì)滲透系數(shù)場中PCE污染Fig.4 PCE saturation distribution in the heterogeneous permeability field under the single source situation

表2 單一污染源非均質(zhì)滲透系數(shù)場中污染羽在空間上的分布特征

3.2 雙重污染源

圖5為均質(zhì)滲透系數(shù)場下，雙重污染源泄漏的PCE運行了2 a(730 d)后，在不同尺度下由UTCHEM模擬的PCE運移結(jié)果，表3為該污染羽的空間矩計算結(jié)果。由表3可知，同單一污染源情況一樣，均質(zhì)情況下，大尺度模型和小尺度模型得到的污染羽的二階矩存在一定誤差，但其零階矩和一階矩誤差很小，表明二源情況下，大尺度模型也能基本再現(xiàn)小尺度模型精細(xì)刻畫污染物的運移情況(圖5)，以此可作為非均質(zhì)情況的對照。

表3 雙重污染源均質(zhì)滲透系數(shù)場中污染羽在空間上的分布特征

圖6為雙重污染源不同滲透系數(shù)場下，PCE運行了2 a(730 d)后，UTCHEM根據(jù)不同的網(wǎng)格剖分情況計算出的PCE的運移結(jié)果，表4為該污染羽的空間矩分析結(jié)果。與單源情況相比，由于非均勻離散中，雙重污染源細(xì)化的網(wǎng)格更多，因此在污染物質(zhì)量估計方面，大尺度模型的準(zhǔn)確度都較單源情況有所提升。雙重污染源得到的分析結(jié)果與單源情況類似，進(jìn)一步印證了拉普拉斯-外殼法的優(yōu)越性。因此，在不同污染情況下，拉普拉斯-外殼法提升后的滲透系數(shù)場都能更好地等效小尺度滲透系數(shù)場情況。

3.3 計算效率

由于算術(shù)平均法和拉普拉斯-外殼法進(jìn)行尺度提升的運算時間均極短，因此將UTCHEM模型的運行時間作為計算成本進(jìn)行對比。本算例采用配置為Intel(R) Core(TM) i5處理器、3.2 GHz、3.4 G內(nèi)存的計算機，具體運行時間見表5。

從表5中可以看出，尺度提升帶來的計算時間減少是很可觀的，并且小尺度模型計算時間越長，尺度提升后節(jié)約的計算時間越多。利用算術(shù)平均進(jìn)行尺度提升后的模型計算時間為原模型的2.47%～4.28%，利用拉普拉斯-外殼法尺度提升后的模型計算時間為原模型的1.85%～3.45%。將二者的運行時間進(jìn)行比較，所有情形下均是拉普拉斯-外殼法尺度提升后的模型計算效率更高，這是由于拉普拉斯-外殼法中外殼區(qū)域不用進(jìn)行計算，因此拉普拉斯-外殼法得到的模型需要計算的網(wǎng)格數(shù)更少，從而具備更高的計算效率。計算效率極大的提高表明，滲透系數(shù)場的尺度提升對運移模型計算成本的降低是非常有效的。

圖6 雙重污染源非均質(zhì)滲透系數(shù)場中PCE污染Fig.6 PCE saturation distribution in the heterogeneous permeability field under the dual sources situation

表4 雙重污染源非均質(zhì)滲透系數(shù)場中污染羽在空間上的分布特征

4 結(jié)論

(1)零階矩的計算結(jié)果表明，算術(shù)平均法和拉普拉斯-外殼法對滲透系數(shù)場進(jìn)行尺度提升后得到的大尺度運移模型都能準(zhǔn)確地估計含水層中污染物質(zhì)量。相較而言，拉普拉斯-外殼法得到的結(jié)果更準(zhǔn)確，且含水層非均質(zhì)性越強，拉普拉斯-外殼法的優(yōu)越性越明顯。

表5 不同尺度提升方法的模型運行時間對比

注：a時間比：尺度提升后模型計算時間與尺度提升前模型計算時間的百分比。

(2)一階矩和二階矩的計算結(jié)果表明，在水流方向上，拉普拉斯-外殼法對質(zhì)心位置的確定和污染羽的運移范圍的計算更準(zhǔn)確；在垂直于水流方向上，兩種計算方法計算效果類似。因此，經(jīng)過拉普拉斯-外殼法尺度提升后得到的運移模型刻畫小尺度模型更準(zhǔn)確，可更好地為后期的含水層污染修復(fù)提供參考。

(3)計算時間的對比結(jié)果表明，算術(shù)平均法和拉普拉斯-外殼法對滲透系數(shù)場進(jìn)行尺度提升都能有效節(jié)省計算成本。