地下水典型非水相液體污染運移模擬的尺度提升研究

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

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

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

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

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

1 研究方法

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

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

(1)大尺度網(wǎng)格(block)中的水頭應等于小尺度網(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)格中的流量應等于小尺度網(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代碼進一步計算地下水中NAPLs污染羽隨時間變化的空間分布。

2 算例應用

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泄漏，并穿過弱透水層進入含水層。污染源分別有以下兩種情形：

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

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

對雙重污染源兩個污染點位置分別進行網(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)小尺度模型精細刻畫污染物的運移情況(圖3)，以此可作為非均質(zhì)情況的對照。

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

表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)小尺度模型精細刻畫污染物的運移情況(圖5)，以此可作為非均質(zhì)情況的對照。

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

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

3.3 計算效率

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

從表5中可以看出，尺度提升帶來的計算時間減少是很可觀的，并且小尺度模型計算時間越長，尺度提升后節(jié)約的計算時間越多。利用算術(shù)平均進行尺度提升后的模型計算時間為原模型的2.47%～4.28%，利用拉普拉斯-外殼法尺度提升后的模型計算時間為原模型的1.85%～3.45%。將二者的運行時間進行比較，所有情形下均是拉普拉斯-外殼法尺度提升后的模型計算效率更高，這是由于拉普拉斯-外殼法中外殼區(qū)域不用進行計算，因此拉普拉斯-外殼法得到的模型需要計算的網(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ù)場進行尺度提升后得到的大尺度運移模型都能準確地估計含水層中污染物質(zhì)量。相較而言，拉普拉斯-外殼法得到的結(jié)果更準確，且含水層非均質(zhì)性越強，拉普拉斯-外殼法的優(yōu)越性越明顯。

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

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

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

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