常 興，駱乾坤，鄧亞平，馬 雷，錢家忠

（合肥工業(yè)大學(xué)資源與環(huán)境工程學(xué)院，安徽 合肥 230009）

地下水環(huán)境中存在大量難以消除的污染物，其中包括重非水相液體(dense nonaqueous-phase liquids，DNAPLs)。因?yàn)槠渑c水不相容，在運(yùn)移過程中由于毛管壓力、黏滯力和重力的共同作用下會殘存在運(yùn)移路徑中，長期存在且難以消除。裂隙介質(zhì)是地下水三大賦存介質(zhì)之一[1]。由于裂隙介質(zhì)的各向異性和非均質(zhì)性明顯強(qiáng)于孔隙介質(zhì)，所以裂隙中地下水受到DNAPLs污染后更難修復(fù)。為此，研究裂隙介質(zhì)DNAPLs 的運(yùn)移規(guī)律對此類污染物的修復(fù)尤為重要。

目前DNAPLs運(yùn)移研究大多集中在孔隙介質(zhì)[2-3]，主要是通過數(shù)值模擬[4-7]和室內(nèi)試驗(yàn)[8-9]來研究孔隙介質(zhì)的非均質(zhì)性[10-12]以及地下水流流速[13]等對DNAPLs運(yùn)移的影響。裂隙介質(zhì)中的DNAPLs運(yùn)移研究還處于起步階段[14-15]，主要研究裂隙的相關(guān)參數(shù)以及裂隙的物理性質(zhì)對DNAPLs運(yùn)移及分布的影響。就目前來說，研究裂隙寬度空間的變異性主要圍繞單裂隙介質(zhì)[16-17]，裂隙隙寬的改變影響裂隙中水流的流態(tài)、流速。從裂隙水流態(tài)的研究出發(fā)延伸至單裂隙中裂隙隙寬的變異性對DNAPLs運(yùn)移的影響，主要通過室內(nèi)試驗(yàn)[18]、數(shù)值模擬[19-23]及求解對應(yīng)數(shù)學(xué)模型[24]開展研究。裂隙隙寬變化，會改變DNAPLs 的流量、流速，最終影響DNAPLs 在裂隙介質(zhì)中的運(yùn)移范圍和空間分布。泄漏速率對DNAPLs運(yùn)移的研究集中在孔隙介質(zhì)中[25]，研究結(jié)果表明DNAPLs 泄漏速率不同，泄漏總量不同，導(dǎo)致污染羽的空間展布不同。關(guān)于泄漏位置對DNAPLs 在裂隙中運(yùn)移影響的研究很少[26]。

綜上所述，國內(nèi)外對DNAPLs運(yùn)移規(guī)律的研究大多在單裂隙或多孔介質(zhì)中進(jìn)行，而對隨機(jī)網(wǎng)絡(luò)裂隙中DNAPLs運(yùn)移規(guī)律及影響因素的研究不足，綜合考慮裂隙寬度空間變異性和泄漏條件對DNAPLs 在裂隙中的運(yùn)移影響的研究更少。因此本文結(jié)合真實(shí)網(wǎng)絡(luò)裂隙數(shù)據(jù)，通過蒙特卡羅方法建立裂隙寬度空間變異性不同的隨機(jī)網(wǎng)絡(luò)裂隙，采用PetraSim軟件模擬并分析裂隙寬度空間變異性、泄漏速率和泄漏位置對DNAPLs運(yùn)移的影響，豐富裂隙介質(zhì)中DNAPLs運(yùn)移機(jī)理研究，為裂隙介質(zhì)中DNAPLs 污染修復(fù)提供模型參考。

1 隨機(jī)網(wǎng)絡(luò)裂隙場的生成

為模擬天然裂隙網(wǎng)絡(luò)中裂隙寬度空間變異性對DNAPLs運(yùn)移的影響，本文基于已有統(tǒng)計(jì)資料[27]，使用蒙特卡羅方法(蒙特卡羅方法可以根據(jù)裂隙各幾何參數(shù)所服從的分布函數(shù)，編程生成隨機(jī)網(wǎng)絡(luò)裂隙。此方法能最大程度模擬真實(shí)裂隙，是目前廣泛使用的生成隨機(jī)網(wǎng)絡(luò)裂隙的方法)生成裂隙寬度空間變異性（方差不同，表明生成的裂隙場和裂隙寬度不同，此不同即為空間變異性）不同的3 組網(wǎng)絡(luò)裂隙物理模型，主要參數(shù)見表1。

表1 不同裂隙場中的裂隙幾何參數(shù)Table 1 Parameter setted for different fracture

將生成的隨機(jī)網(wǎng)絡(luò)裂隙坐標(biāo)化。步驟如下：

第一步：按表1 中的裂隙幾何參數(shù)生成隨機(jī)線段，設(shè)定隨機(jī)線段結(jié)果輸出區(qū)域大小為100 cm×100 cm 的正方形區(qū)域，對生成的線段進(jìn)行統(tǒng)計(jì)分析；

第二步：進(jìn)行像素掃描，將上述區(qū)域離散為水平和垂直間距均為1 cm 的小網(wǎng)格共10 000 個(gè)，使用Python語言編程識別生成的線段所覆蓋的網(wǎng)格塊，并對其進(jìn)行標(biāo)記；

第三步：輸出所有標(biāo)記線段首尾網(wǎng)格的坐標(biāo)及寬度信息，生成隨機(jī)網(wǎng)絡(luò)裂隙場?；诓煌严秾挾瓤臻g變異性（σ2）的3 組隨機(jī)網(wǎng)絡(luò)裂隙場如圖1所示。

圖1 基于不同裂隙寬度空間變異性下的隨機(jī)網(wǎng)絡(luò)裂隙場Fig.1 Stochastic network fractures with spatial variability of different widths

2 算例介紹

基于生成的隨機(jī)網(wǎng)絡(luò)裂隙場，建立DNAPLs運(yùn)移數(shù)值模擬模型。假設(shè)模擬區(qū)域?yàn)?00 cm×100 cm 的二維XZ剖面，X軸方向和Z軸方向剖分步長均為1 cm，共離散10 000 個(gè)網(wǎng)格。地下水面位于模型頂部。左右邊界均為定水頭邊界，上下邊界為無通量邊界。模型首先達(dá)到水-氣平衡，初始水相飽和度Sw為1.0。模型區(qū)域內(nèi)水力梯度J設(shè)置為0.01，水流方向?yàn)閺淖笙蛴遥▓D2）。因?yàn)樗穆纫蚁≒CE）是目前地下水中普遍存在的污染物，且已經(jīng)成為國內(nèi)外專家研究地下水生物修復(fù)的熱點(diǎn)問題，所以選取PCE 作為DNAPLs 示例。由于PCE 很難溶解于水相中(溶解度約為150 mg/L)，故本文不考慮PCE 的溶解及吸附，假設(shè)其完全以非水相形式存在。由于本次模擬溫度設(shè)定為20℃，所以不考慮熱物理相關(guān)參數(shù)。PCE 以線源方式（共3 個(gè)網(wǎng)格）從模型頂部3 個(gè)可能位置（圖2）以恒定的速率泄漏。綜合考慮模型尺寸和泄漏速率，設(shè)定總模擬時(shí)間為600 s。模型主要參數(shù)設(shè)置見表2。

圖2 模擬區(qū)域概念模型圖Fig.2 Conceptual model of simulated area

表2 模型參數(shù)值Table 2 Model calculation parameters

3 試驗(yàn)方法

本次采用多相流數(shù)值模擬軟件PetraSim 進(jìn)行模擬。將生成的隨機(jī)網(wǎng)絡(luò)裂隙讀入PetraSim 模擬軟件，進(jìn)行不同裂隙寬度空間變異性、泄漏速率和泄漏位置下PCE 的運(yùn)移模擬。通過空間矩評估PCE 污染羽的空間變化特征，以此分析PCE 的運(yùn)移規(guī)律，其中分別用零階矩、一階矩和二階矩來表征PCE 的總質(zhì)量、質(zhì)心位置和空間展布。由于本例不考慮污染物運(yùn)移過程中的化學(xué)反應(yīng)，所以污染物的總質(zhì)量（零階矩）在運(yùn)移過程中保持不變。

3.1 裂隙寬度空間變異性

選擇3 種不同裂隙寬度方差（σ2=0.25，0.50，0.75）的網(wǎng)絡(luò)裂隙場進(jìn)行模擬，泄漏位置為圖2 中位置2，并對模擬結(jié)果進(jìn)行對比分析。

3.2 泄露速率

在同一網(wǎng)絡(luò)裂隙中，選取3 種泄漏速率在裂隙場中進(jìn)行模擬，泄漏位置為圖2 中泄漏位置2。高、中、低泄露速率分別為1.4×10-4，1.4×10-5，1.4×10-6kg/s[25]。3 種裂隙場中的模擬結(jié)果規(guī)律相似，受文章篇幅所限，僅選擇 σ2=0.25 的裂隙場作為代表，并對模擬結(jié)果進(jìn)行對比分析。

3.3 泄漏位置

在同一網(wǎng)絡(luò)裂隙場中，設(shè)置3 個(gè)泄漏位置（圖2），進(jìn)行不同泄漏位置PCE 運(yùn)移模擬。3 種裂隙場中的模擬結(jié)果規(guī)律相似，受文章篇幅所限，僅選擇 σ2=0.75的裂隙場作為代表，并對模擬結(jié)果進(jìn)行對比分析。

4 結(jié)果

圖3 顯示了不同裂隙寬度空間變異性的網(wǎng)絡(luò)裂隙中PCE 在模擬末時(shí)刻的飽和度分布，裂隙寬度空間變異性的改變使得PCE 空間展布發(fā)生改變，PCE 到達(dá)模型底部的時(shí)間不同，PCE 飽和度不同。

圖3 網(wǎng)絡(luò)裂隙不同空間變異性下末時(shí)刻PCE 的飽和度分布Fig.3 Corresponding saturation distribution of DNAPLs for different spatial variability

圖4 為3 種不同泄漏速率下在模擬末時(shí)刻的PCE 飽和度分布圖。由于泄漏速率不同，相同時(shí)間內(nèi)釋放的PCE 總量不同，空間展布不同，最終PCE 飽和度不同。

圖4 不同泄漏速率條件下末時(shí)刻PCE 的飽和度分布Fig.4 Corresponding saturation distribution of DNAPLs for different leakage rate

圖5 為同一模型不同泄漏位置末時(shí)刻的PCE 飽和度分布圖。在同一模型中，泄漏位置的不同導(dǎo)致相同時(shí)間內(nèi)PCE 污染羽分布范圍不同，飽和度不同。

圖5 不同泄漏位置條件下末時(shí)刻PCE 的飽和度分布Fig.5 Corresponding saturation distribution of DNAPLs for different leakage location

5 分析與討論

5.1 裂隙寬度空間變異性對PCE 運(yùn)移范圍的影響

裂隙寬度空間變異性不同，PCE 運(yùn)移范圍不同。分析原因，由于寬度較大裂隙內(nèi)毛管壓力較小，構(gòu)成優(yōu)勢通道，控制PCE 運(yùn)移路徑，PCE 能夠很快到達(dá)模型底部邊界，所以PCE 多集中在寬度較大的裂隙中。同時(shí)，由于裂隙寬度在垂向和橫向上的分布不同，致使PCE 飽和度空間展布不同。圖6（a）顯示了具有不同裂隙空間變異性下的PCE 在Z軸的質(zhì)心位置。空間變異性大的網(wǎng)絡(luò)裂隙（σ2=0.75），PCE飽和之后通過寬度較大的裂隙通道繼續(xù)運(yùn)移至模型底部其他橫向分布的裂隙中，因而PCE 質(zhì)心位置隨時(shí)間降低減緩。裂隙寬度空間變異性小( σ2=0.50，0.25)時(shí)，泄漏的PCE 到達(dá)模型底部后大多在運(yùn)移通道內(nèi)堆積。另外，由于裂隙的寬度突變和交叉角度的改變使得局部裂隙通道內(nèi)蓄積的PCE 的飽和度增大，減小了PCE 的運(yùn)移速率，因而PCE 質(zhì)心位置在模擬后期呈上升趨勢。圖6（b）顯示了裂隙寬度不同空間變異性的網(wǎng)絡(luò)裂隙中污染羽圍繞質(zhì)心的空間分布，橫軸表示不同方差的裂隙場，縱軸表示不同裂隙場對應(yīng)的橫向，縱向的空間分布。不同裂隙寬度方差條件下污染羽繞質(zhì)心的垂向展布與橫向展布均存在差異。裂隙空間變異性的改變，使得PCE運(yùn)移路徑和空間展布方向不同。裂隙空間變異性小，裂隙寬度離散程度較小，隙寬相對集中，PCE 通過裂隙通道向周圍運(yùn)移，PCE 飽和度的空間變異性較小。隨著裂隙空間變異性的增大，裂隙隙寬離散程度增大，相對較窄和相對較寬的裂隙占比增大，相對較寬的裂隙通道構(gòu)成優(yōu)勢通道，PCE 通過優(yōu)勢通道快速向周圍運(yùn)移，使得PCE 飽和度的空間變異性增大。

圖6 不同空間變異性下污染羽空間分布特征Fig.6 Spatial moments of DNAPLs contaminant plume for different spatial variability

5.2 泄漏速率對PCE 運(yùn)移速率和空間分布的影響

圖7（a）顯示了3 種不同泄漏速率情形下PCE 在Z軸的質(zhì)心位置。從圖中可得，高泄漏速率由于流量大，相同時(shí)間內(nèi)PCE 泄露總量最多，受重力影響最大，因此能最快克服毛管壓力，集中向下運(yùn)移，因此PCE能最快到達(dá)模型底部并在裂隙通道內(nèi)堆積，導(dǎo)致質(zhì)心位置快速降低。相較于高泄漏速率，中泄漏速率相同時(shí)間內(nèi)PCE 泄漏總量減少，導(dǎo)致其重力不易克服毛管壓力，雖然污染羽能到達(dá)模型底部并在模型底部裂隙通道內(nèi)堆積，但是質(zhì)心位置較高泄漏速率降低較慢。低泄漏速率由于流量最小，PCE 在運(yùn)移過程中受毛管壓力影響最大，導(dǎo)致運(yùn)移速率最慢，運(yùn)移深度最小，污染羽大多集中在注入點(diǎn)下方，質(zhì)心位置降低最慢。圖7（b）顯示了3 種不同泄漏速率條件下PCE 污染羽圍繞質(zhì)心的空間分布。由于相同時(shí)間內(nèi)高泄漏速率釋放的總量最多，PCE 在裂隙相對較寬的通道最先飽和，待飽和之后，PCE 會進(jìn)入相對較窄的裂隙通道中，所以空間展布范圍較中、低泄漏速率大。相同時(shí)間內(nèi)，中泄漏速率情形下PCE 與高泄漏速率均運(yùn)移至模型底部并形成高飽和度區(qū)；但由于入滲速率較低，到達(dá)模型底部的時(shí)間較長，所以相同時(shí)間內(nèi)中泄漏速率情形下PCE 空間展布范圍變化較大。低泄漏速率情形下由于流量小，相同時(shí)間內(nèi)釋放總量最小，難以克服毛管壓力，因此PCE 集中在泄漏點(diǎn)下方且空間展布范圍最小。

圖7 不同泄漏速率下空間分布特征Fig.7 Spatial moments of DNAPLs contaminant plume for different leakage rate

5.3 泄漏位置對PCE 質(zhì)心位置和空間分布的影響

泄露位置不同，PCE 污染羽的質(zhì)心位置和空間分布不同。圖8（a）為同一模型不同泄漏位置PCE 在Z軸質(zhì)心位置圖。從圖中可得，由于較寬的裂隙通道內(nèi)毛管壓力較小，利于優(yōu)勢通道的形成，從而控制PCE 的運(yùn)移路徑。PCE 通過優(yōu)勢通道快速運(yùn)移，使得污染羽在垂直方向上的分布范圍明顯增大，所以PCE能夠很快運(yùn)移至模型底部。因此，即使同一模型中，由于泄漏位置不同，污染物運(yùn)移的重力方向的空間變異性也不相同，導(dǎo)致PCE 質(zhì)心位置也不相同。圖8（b）為同一模型不同泄漏位置的PCE 圍繞質(zhì)心的空間分布。PCE 受重力作用向下運(yùn)移，同一模型不同泄漏位置其重力方向網(wǎng)絡(luò)裂隙分布不同，污染物到達(dá)模型底部時(shí)間不同。同樣由于網(wǎng)絡(luò)裂隙空間變異性的不同，PCE 空間分布也不相同。

圖8 不同泄漏位置空間分布特征Fig.8 Spatial moments of DNAPLs contaminant plume for different leakage location

6 結(jié)論

（1）隨機(jī)網(wǎng)絡(luò)裂隙寬度的空間變異性是影響DNAPLs 飽和度分布的一個(gè)重要因素。網(wǎng)絡(luò)裂隙寬度空間變異性的增大，增大了網(wǎng)絡(luò)裂隙中相對較寬的裂隙的比例，相對較寬的裂隙會形成優(yōu)勢通道，DNAPLs 通過這些優(yōu)勢通道快速向周邊運(yùn)移，導(dǎo)致水平方向和垂直方向上DNAPLs 的一階矩和二階矩發(fā)生改變。

（2）隨機(jī)網(wǎng)絡(luò)裂隙中DNAPLs運(yùn)移范圍受泄漏速率影響。隨著泄漏速率的增大，相同時(shí)間內(nèi)，污染物泄漏總量增大，運(yùn)移速率加快，空間分布范圍增大，一階距和二階矩變化增大。

（3）同一模型中，泄漏位置不同，導(dǎo)致DNAPLs運(yùn)移重力方向的裂隙介質(zhì)的空間變異性不同，DNAPLs運(yùn)移路徑和污染羽分布范圍也不同。