水位波動下包氣帶透鏡體影響LNAPL遷移的數(shù)值模擬研究

潘明浩 ，時 健 ，左 銳 ，趙 曉 ，劉嘉蔚 ，薛鎮(zhèn)坤 ，王金生 ，胡立堂

（1.北京師范大學水科學研究院，北京 100875；2.地下水污染控制與修復教育部工程研究中心，北京 100875）

輕質(zhì)非水相液體（Light Non-Aqueous Phase Liquid，LNAPL）是土壤及地下水環(huán)境中有機污染問題的重要來源，往往具有持續(xù)危害性[1-2]。土壤包氣帶是LNAPL下滲遷移過程中的重要通道[3]，其含水量、氧化還原條件等隨土壤深度變化，使得LNAPL在包氣帶中的多相流遷移行為變得較為復雜[4-5]。一般認為，LNAPL進入包氣帶后首先以垂向遷移為主，到達地下水位上方的毛細帶后將在毛細力作用下驅(qū)替孔隙水，發(fā)生側(cè)向遷移擴展[6-8]。LNAPL在包氣帶中的遷移行為受黏度、密度等流體性質(zhì)的控制[9-10]以及介質(zhì)巖性、地下水位波動等外部因素的影響[7,11]，其中非均質(zhì)性（如透鏡體、非均質(zhì)夾層、裂隙等）[12-13]以及水動力條件[14-15]在很大程度上將影響LNAPL的遷移規(guī)律和相態(tài)分布特征。

國內(nèi)外學者廣泛采用物理模型試驗或數(shù)值模擬方法探究包氣帶內(nèi)的LNAPL遷移過程。非均質(zhì)條件下的介質(zhì)含水率差異以及毛細力作用均被認為是控制LNAPL遷移規(guī)律的重要因素。Schroth[16]在通過砂箱實驗探究LNAPL在非均質(zhì)包氣帶中的遷移過程時發(fā)現(xiàn)含水量對LNAPL分布存在顯著影響。束善治等[7]、陶佳輝等[17]則通過數(shù)值模擬探究不同地層結(jié)構(gòu)中LNAPL的遷移過程，發(fā)現(xiàn)介質(zhì)含水率差異決定了非均質(zhì)結(jié)構(gòu)對LNAPL運動的阻礙或促進作用。此外，Simantiraki[18]與Wipfler[19]分別運用圖像分析與數(shù)值模擬手段分析LNAPL遷移過程，二者均認為毛細力是污染物遷移過程的重要控制因素，LNAPL垂向運動受阻可能與界面形成的毛細屏障有關。水位波動是土壤及地下水污染問題中的另一關鍵因素，對于LNAPL污染物而言，水位波動將顯著影響氣-水-LNAPL三相之間的相互驅(qū)替，進而改變污染物遷移及相態(tài)分布過程[20]。羅凌云[11]通過實驗探究水位波動條件下柴油進入包氣帶的遷移行為時，發(fā)現(xiàn)LNAPL在水位上升時將會優(yōu)先驅(qū)替多孔介質(zhì)中的氣體，在水位下降時則將優(yōu)先驅(qū)替水相；楊明星[15]結(jié)合實驗及數(shù)值模擬手段對石油組分在水位波動帶內(nèi)分異演化機理進行研究，發(fā)現(xiàn)污染物在水位波動條件下受“界面效應”影響，趨向于滯留在較低滲透性的介質(zhì)中，水位波動將加劇LNAPL吸附、溶解等過程并使其對水的驅(qū)替過程減弱。因此可以說，地下水位波動將顯著影響LNAPL的遷移過程[21]。

為探究土壤包氣帶內(nèi)存在透鏡體結(jié)構(gòu)這一復雜環(huán)境條件下LNAPL的遷移過程及其相態(tài)分布，本文基于北京市通州區(qū)張家灣包氣帶試驗基地的場地狀況及其地下水位波動背景，采用TOUGH2（transport of unsaturated groundwater and heat）程序的TMVOC模塊進行算例研究，構(gòu)建二維非均質(zhì)剖面中LNAPL入滲遷移的多相流模型，以透鏡體介質(zhì)巖性、地下水位波動作為共同影響因素，模擬非均質(zhì)包氣帶內(nèi)LNAPL的遷移規(guī)律、分布特征及揮發(fā)過程，以期為復雜條件下的石油污染場地修復提供理論參考。

1 材料與方法

1.1 模型背景

TOUGH2程序可用于模擬多孔介質(zhì)中的地下水流、非等溫熱量運移等多相流過程[22]，TMVOC模塊是該程序后續(xù)版本中用于模擬NAPLs在飽和-非飽和介質(zhì)中遷移、VOCs（揮發(fā)性有機污染物）溶解與揮發(fā)等過程的模塊[23]，采用積分有限差分法進行模擬區(qū)域的空間離散[24]。

通州區(qū)為北京城市副中心，近年來該地區(qū)地下水位總體呈上升趨勢，但伴有水位的持續(xù)波動。在水位波動背景下，其土壤地下水中有機物污染問題愈發(fā)受到關注。通州區(qū)張家灣包氣帶試驗基地位于永定河與潮白河沖洪積扇中下部的交匯地帶，屬于第四系沉積地層，其上部地層巖性主要包括細砂、粉細砂、中砂等，場地介質(zhì)非均質(zhì)性較強且局部存在厚度為幾厘米至幾十厘米不等的微觀非均質(zhì)透鏡體結(jié)構(gòu)。在此基礎上，本文設計包氣帶內(nèi)存在透鏡體的非均質(zhì)地層結(jié)構(gòu)，假定LNAPL污染物由近地表泄漏，采用TMOVC模塊建立XZ方向的二維剖面數(shù)值模型。如圖1所示，模型尺寸為90 cm×10 cm×110 cm，其中X方向剖分為18列網(wǎng)格，Z方向剖分為22層網(wǎng)格，網(wǎng)格水平及垂直間距均為5 cm。透鏡體區(qū)域坐標范圍為X=30～60 cm；Z=-70～-75 cm，LNAPL泄漏源坐標范圍為X=40～50 cm；Z=0～-5 cm。此外，模型共設置6個觀察點用于觀察透鏡體周圍的LNAPL運移狀況。

圖1 概念模型示意圖Fig.1 Schematic diagram of the conceptual model

1.2 模型參數(shù)

本文采用中砂介質(zhì)作為包氣帶內(nèi)的主要介質(zhì)，分別采用細砂及粗砂作為透鏡體介質(zhì)。采用Stone模型作為相對滲透率函數(shù)，其中Swr為水相殘余飽和度；Snr為NAPL相殘余飽和度；Sgr為氣相殘余飽和度；n為擬合參數(shù)。采用van Genuchten模型作為毛細壓力函數(shù)，其中m=1-1/n；Slr為液相殘余飽和度；Sls為飽和液相的飽和度；P0為進氣值；Pmax為最大進氣壓力值。介質(zhì)相關各參數(shù)如表1所示。

表1 介質(zhì)相關參數(shù)[17,25-26]Table 1 Medium-relevant parameters

苯系物（BTEX）是較為常見的LNAPL污染物[27]，本文采用甲苯作為算例中的LNAPL污染物，考慮甲苯的揮發(fā)及溶解等相態(tài)傳質(zhì)過程，相關參數(shù)如表2所示。模型在等溫模式下運行，溫度恒定為25 ℃。

表2 甲苯相關參數(shù)[24]Table 2 Toluene-relevant parameters

1.3 邊界條件及初始條件

模型頂部單元格設定為大氣邊界（Dirichlet邊界），垂向范圍為Z=0～5 cm，壓力值恒定為大氣壓1.013×105Pa，氣體飽和度為0.99。模型底部一層單元格為定壓Dirichlet邊界，其深度范圍Z=-100～-105 cm，初始壓力值按照地下水位-100 cm計算并設定為1.017 9×105Pa。模型構(gòu)建完成后首先計算得到初始水分分布，圖2所示分別為中砂介質(zhì)包氣帶內(nèi)存在細砂、粗砂透鏡體時模型中的初始水分分布情況。

圖2 各透鏡體模型中的含水率變化規(guī)律Fig.2 Variation pattern of water content in each model

2 結(jié)果與討論

2.1 不同巖性透鏡體條件下自由相遷移分布

本節(jié)設計包氣帶內(nèi)存在細砂、粗砂透鏡體兩種模型算例，分別以LNAPL飽和度SNAPL、LNAPL在水相中的質(zhì)量分數(shù)x(aq)作為表征自由相、溶解相分布的定量指標，分析水位恒定時LNAPL在包氣帶內(nèi)的遷移分布。LNAPL釋放時間為0～180 min，釋放速率為4.483 mL/min。圖3及圖4分別代表包氣帶內(nèi)存在細砂、粗砂透鏡體時LNAPL的遷移過程。

圖3 細砂透鏡體中不同時刻LNAPL飽和度分布Fig.3 Distribution of LNAPL saturation at different moments in the fine-sand lens model

圖4 粗砂透鏡體中不同時刻LNAPL飽和度分布Fig.4 Distribution of LNAPL saturation at different moments in the coarse-sand lens model

對于細砂透鏡體模型，LNAPL遷移到達透鏡體上方時存在明顯的蓄積和橫向擴展現(xiàn)象，180 min前透鏡體兩側(cè)SNAPL幾乎為0，且透鏡體上方（觀察點2）及下方（觀察點5）SNAPL相對較大，表明部分LNAPL已優(yōu)先以“穿過透鏡體”形式向下遷移。對于粗砂透鏡體模型，LNAPL更多以垂向遷移為主，表現(xiàn)為150 min時LNAPL分布深度超過細砂模型約17%，105～180 min時LNAPL在透鏡體下方發(fā)生橫向擴展。污染物釋放階段的模擬結(jié)果表明，粗砂透鏡體成為LNAPL垂向遷移的“優(yōu)勢通道”，細砂透鏡體則表現(xiàn)出對LNAPL垂向運動的阻礙，蓄積于透鏡體上方的LNAPL，在其橫向擴展超過透鏡體寬度之前以“穿過透鏡體”為主。上述不同的遷移過程與介質(zhì)參數(shù)及含水率差異存在較大聯(lián)系：由于水-氣-NAPL三相體系中，NAPL驅(qū)替水相、氣相發(fā)生運移的難易程度與NAPL/水界面的進氣壓力值密切相關[7,28]，故粗砂介質(zhì)相對較小的進氣壓力值Pmax（表1）決定了LNAPL更容易驅(qū)替粗砂介質(zhì)中的水相、氣相，因此LNAPL到達粗砂透鏡體附近時具有更強的遷移能力。其次，細砂介質(zhì)較強的持水性顯著改變包氣帶內(nèi)水分分布，水分占據(jù)孔隙空間使細砂透鏡體對LNAPL垂向遷移表現(xiàn)出一定阻礙作用；而粗砂透鏡體則使-70～-90 cm深度范圍內(nèi)的含水率變化更加顯著，促使LNAPL在透鏡體及其下方發(fā)生顯著的橫向擴展。

污染物停止釋放（180 min）后，粗砂透鏡體模型中LNAPL遷移過程與釋放階段一致，而對于細砂透鏡體模型，180～200 min內(nèi)透鏡體兩側(cè)4個觀察點中SNAPL均呈現(xiàn)“迅速上升-逐漸降低”規(guī)律（圖5），表明LNAPL橫向擴展超過細砂透鏡體寬度后將轉(zhuǎn)而主要以“繞流”形式運移，釋放停止約5.8 h后透鏡體上方污染物橫向擴展范圍逐漸減小，“繞流”減弱。

圖5 細砂透鏡體周圍各觀察點SNAPL對比（0～360 min）Fig.5 Comparison of SNAPL at each observation point around the fine-sand lens (0～360 min)

就污染物分布狀況而言，180 min后兩模型-85 cm深度以下的LNAPL橫向擴展情況趨于一致，且兩透鏡體模型在該深度范圍內(nèi)的含水率分布規(guī)律一致（圖2），故可推斷包氣帶內(nèi)的含水率分布規(guī)律是控制LNAPL分布形態(tài)的主要因素。透鏡體巖性差異亦將影響兩模型中的LNAPL分布：由于細砂介質(zhì)孔隙度較大且滲透性較差，10～24 h部分LNAPL在透鏡體內(nèi)部產(chǎn)生滯留現(xiàn)象；而粗砂介質(zhì)較好的滲透性則決定了其對污染物的滯留能力較差，大部分LNAPL集中分布于透鏡體下方。模擬至24 h時，粗砂透鏡體模型中污染物分布深度范圍相比細砂透鏡體模型小28.5%，包氣帶內(nèi)LNAPL分布面積（約1 950 cm2）相比細砂透鏡體模型（約2 300 cm2）小15.2%。

2.2 水位波動條件下自由相遷移分布

各種原因造成的地下水補排不平衡將導致水位短時波動（如地表水-地下水相互轉(zhuǎn)化、降水入滲等[29]），使包氣帶內(nèi)LNAPL遷移分布狀況復雜化[15]。本節(jié)以細砂、粗砂透鏡體模型作為水位波動算例基礎，保持LNAPL瞬時釋放（4.483 mL/min，釋放180 min），通過下層18個網(wǎng)格作為控制水位波動的源匯項，在250 min時控制水位由-100 cm升至-80 cm，并于12 h時恢復至-100 cm。圖6(a)、(b)分別代表水位波動條件下細砂、粗砂透鏡體模型中LNAPL遷移分布。

圖6 地下水位波動情景中不同時刻LNAPL飽和度分布Fig.6 Distribution of LNAPL saturation at different moments in the water table fluctuation scenario

對于細砂透鏡體模型，水位上升階段污染物在浮力作用下向上抬升，至12 h時LNAPL主要分布于透鏡體上方并呈連續(xù)池狀（Pool）[30]，此時透鏡體上方（觀察點2）SNAPL相比水位恒定時上升0.15。同時，水位升高使LNAPL遷移過程發(fā)生顯著變化，圖7所示為100 min至24 h內(nèi)透鏡體上方（觀察點2）、下方（觀察點5）以及透鏡體內(nèi)部SNAPL變化狀況，對比可知，水位升高時細砂透鏡體內(nèi)部SNAPL相比水位恒定情景減小0.13～0.20，透鏡體兩側(cè)SNAPL平均值相比水位恒定時增大257%～386%，表明水位升高時LNAPL將主要以“繞流”形式在包氣帶內(nèi)運動，其原因為水位升高時細砂透鏡體含水率增大至35.8%，更加顯著的透鏡體高含水率特征導致其對LNAPL遷移的阻礙作用增強。

圖7 水位恒定與波動條件下細砂透鏡體附近SNAPL對比Fig.7 Comparison of SNAPL near the fine-sand lens under the stationary and fluctuating water table conditions

水位降低時（12～24 h）包氣帶內(nèi)含水率分布恢復為水位波動前狀況，已大范圍橫向擴展的LNAPL在重力作用下向下遷移，故18 h后透鏡體內(nèi)部及其周圍SNAPL變?yōu)榕c水位恒定情景一致（圖7）。經(jīng)水位波動，-60～-80 cm深度內(nèi)LNAPL分布范圍大幅增加，導致包氣帶內(nèi)LNAPL分布面積相比水位恒定情景中增大約51%。此外，24 h時部分LNAPL蓄積于地下水位以上呈透鏡體狀，其橫向分布范圍、SNAPL最大值相比水位恒定時分別增大12.5%、5.6%，表明水位開始降低時部分LNAPL仍以“繞流”形式運動，由此造成更多污染物向下部區(qū)域遷移。

對于粗砂透鏡體模型，水位升高時大部分LNAPL分布于透鏡體內(nèi)部及上方，結(jié)合圖8可知，250 min后粗砂透鏡體內(nèi)部及其上方SNAPL均出現(xiàn)上升，透鏡體內(nèi)部SNAPL相比水位恒定時平均增大0.19，其原因為：水位升高使透鏡體進一步成為LNAPL遷移分布過程的相對“優(yōu)勢空間”。此外，水位升高時透鏡體兩側(cè)SNAPL相比水位恒定時平均增大0.13～0.16，表明水位抬升使LNAPL轉(zhuǎn)為在透鏡體區(qū)域發(fā)生橫向遷移，透鏡體所處的-70～-75 cm深度內(nèi)污染物橫向擴展范圍相對增大約64.3 cm。

圖8 水位恒定與波動條件下粗砂透鏡體附近SNAPL對比Fig.8 Comparison of SNAPL near the coarse-sand lensunder the stationary and fluctuating water table conditions

水位降低過程中，粗砂透鏡體內(nèi)部及其上方聚集的LNAPL向透鏡體下方遷移，至24 h時同樣在地下水位附近蓄積形成LNAPL污染暈。經(jīng)歷水位波動后，該模型包氣帶內(nèi)污染物分布面積相較水位恒定時增大約63%，對比細砂透鏡體模型可知，水位波動對粗砂透鏡體模型中污染物分布的影響更加顯著。

2.3 溶解相分布特征

以LNAPL在水相中的質(zhì)量分數(shù)x(aq)作為表征LNAPL溶解相分布特征的定量指標，圖9（a）（b）分別代表水位恒定時，兩種透鏡體模型中溶解相LNAPL的分布變化情況。LNAPL釋放階段兩模型中溶解相分布差異主要體現(xiàn)在透鏡體附近，細砂透鏡體模型觀察點1、3中x(aq)平均值均為6.98×10-5，相比粗砂透鏡體模型高81%。其次，由于粗砂透鏡體的“優(yōu)勢通道”作用以及粗砂較好的滲透性，至180 min時粗砂透鏡體下方的溶解相分布面積相比細砂透鏡體大46%左右。釋放結(jié)束后，隨自由相的遷移，兩模型內(nèi)溶解相均出現(xiàn)較大擴展，至24 h時細砂透鏡體模型包氣帶內(nèi)溶解相分布面積（約2 600 cm2）相比粗砂透鏡體模型（約2 150 cm2）大20.9%。上述結(jié)果表明，兩模型包氣帶內(nèi)溶解相的分布規(guī)律均與自由相一致，但其分布范圍的差異相比自由相表現(xiàn)得更為顯著。

圖9 細砂、粗砂透鏡體中不同時刻LNAPL在水相中的質(zhì)量分數(shù)x(aq)Fig.9 Distribution of x(aq) at different moments in the fine sand-lens model

地下水位波動條件下，溶解相分布特征亦與自由相LNAPL密切相關，水位升高時兩模型中自由相均集中分布于-70～-85 cm，由此使該深度成為水位波動條件下溶解相分布最廣泛的區(qū)域。水位波動使溶解相分布狀況趨于一致：12 h時兩模型在-60 cm以下的溶解相分布面積僅相差約100 cm2，至24 h時兩模型的溶解相分布差異僅體現(xiàn)在包氣帶上部。

2.4 LNAPL揮發(fā)過程

LNAPL通過質(zhì)量交換進入氣相，將以揮發(fā)方式向上穿過包氣帶到達地表[31]，本節(jié)通過模型頂部單元格計算得到的平均氣相VOC通量數(shù)據(jù)表征LNAPL揮發(fā)，圖10及圖11所示分別為水位恒定、水位波動條件下兩模型中揮發(fā)通量的變化情況。

圖10 不同透鏡體模型中LNAPL揮發(fā)通量對比Fig.10 Comparison of LNAPL volatile flux in each lens structure model

圖11 不同水位情景中LNAPL揮發(fā)通量對比Fig.11 Comparison of LNAPL volatile flux in different water table scenarios

水位恒定時，由于污染物橫向擴展阻礙了氣體-LNAPL質(zhì)量交換，因此兩模型中揮發(fā)通量在此階段均呈降低趨勢（圖10）。此階段兩模型揮發(fā)過程存在一定區(qū)別：結(jié)合圖3、圖4對比揮發(fā)通量可知，由于粗砂透鏡體對LNAPL垂向運動呈“優(yōu)勢通道”作用，透鏡體下方污染物分布范圍更大，故粗砂透鏡體模型中溶解過程更顯著，其揮發(fā)通量平均值（1.40×10-6mol/(s·m2)）相比細砂模型（1.53×10-6mol/(s·m2)）小8.5%。LNAPL釋放停止后，兩模型中揮發(fā)過程均呈現(xiàn)出“先減弱-后增強”趨勢，并在23 h左右達到峰值（1.3×10-6mol/(s·m2)）。結(jié)合氣體質(zhì)點方向分析，釋放停止后包氣帶上部的LNAPL對水、氣驅(qū)替減緩，部分氣體快速向原LNAPL下滲路徑上運動，填補多孔介質(zhì)中LNAPL減少所產(chǎn)生的“空缺”，且由于180 min后包氣帶下部的LNAPL分布范圍擴大，故揮發(fā)通量上升。

地下水位波動時，由于水位抬升打破飽和帶內(nèi)原有的三相平衡狀態(tài)，且水位升高時LNAPL將優(yōu)先驅(qū)替多孔介質(zhì)內(nèi)的氣體[11]，故水位開始抬升時細砂、粗砂透鏡體模型中揮發(fā)通量分別迅速上升至9.64×10-6mol/(s·m2)、1.09×10-5mol/(s·m2)，此時分別于細砂透鏡體兩側(cè)、粗砂透鏡體上部觀察到揮發(fā)通量快速增大，部分LNAPL迅速以氣相形態(tài)向上逸散。揮發(fā)通量快速到達峰值后，由于兩模型中LNAPL均在透鏡體區(qū)域橫向擴展，阻礙了氣體-LNAPL質(zhì)量交換過程，故此后模型中揮發(fā)通量呈逐漸降低趨勢。兩透鏡體模型在水位上升階段（250 min至12 h）的平均揮發(fā)通量相比水位恒定時分別增大124%、126%，表明水位上升將大幅增強LNAPL的揮發(fā)過程。水位降低時兩模型的揮發(fā)通量均出現(xiàn)小幅度波動，結(jié)合氣體質(zhì)點運動方向分析，隨著水位降低，模型下部介質(zhì)孔隙內(nèi)LNAPL及水分均向下遷移，模型上部揮發(fā)相隨之轉(zhuǎn)為向下部運動，故模型上部揮發(fā)通量均有所降低。

3 結(jié)論

（1）地下水位恒定時，介質(zhì)含水率差異導致粗砂透鏡體表現(xiàn)出對LNAPL垂向遷移的“優(yōu)勢通道”作用，細砂透鏡體則起不完全阻礙作用，使LNAPL呈現(xiàn)“蓄積穿透-橫向擴展-繞流”遷移過程。細砂、粗砂透鏡體模型中LNAPL主要分布于透鏡體內(nèi)部、下方，兩模型包氣帶內(nèi)污染物分布面積相差約15.2%。

（2）地下水位波動將使粗砂透鏡體進一步成為LNAPL遷移分布的“優(yōu)勢空間”，而其導致的細砂透鏡體含水率變化則使“繞流”顯著增強。水位波動使-60～-80 cm深度內(nèi)LNAPL分布范圍大幅增加，細砂、粗砂模型包氣帶內(nèi)LNAPL分布面積相比水位恒定情景分別增大51%、63%。

（3）溶解相分布規(guī)律與自由相遷移密切相關，水位恒定時兩模型中溶解相分布差異主要體現(xiàn)在透鏡體附近，至24 h時溶解相分布面積相差20.9%，地下水位波動則使該差距顯著減小至約100 cm2。

（4）氣體-LNAPL接觸條件、LNAPL分布狀況決定LNAPL揮發(fā)強度變化，兩模型中揮發(fā)強度變化均呈“先減弱-后增強”趨勢，前180 min粗砂透鏡體模型平均揮發(fā)量相比細砂透鏡體模型小8.5%。水位波動使包氣帶內(nèi)三相平衡狀態(tài)受擾動，主要表現(xiàn)在水位抬升階段LNAPL揮發(fā)增強，該階段兩模型中平均揮發(fā)量相比水位恒定時分別增大124%、126%。