扶曉琴,施小清,蔣建國(guó),吳吉春

寒區(qū)多孔介質(zhì)中LNAPL遷移的多場(chǎng)耦合模擬

扶曉琴,施小清*,蔣建國(guó),吳吉春

(南京大學(xué)地球科學(xué)與工程學(xué)院,表生地球化學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210023)

建立了包含水、氣、NAPL、固(冰)四相的溫度場(chǎng)、流場(chǎng)及化學(xué)場(chǎng)多場(chǎng)耦合模型,適用于定量模擬LNAPL泄漏過程中發(fā)生凍結(jié)后的遷移分布.通過與非凍條件下LNAPL遷移模型對(duì)比,模擬結(jié)果顯示,凍結(jié)條件下因冰的形成LNAPL會(huì)發(fā)生與非凍條件下不同的遷移行為,冰的形成阻礙LNAPL下滲,并產(chǎn)生凍結(jié)勢(shì)使LNAPL向地表方向遷移.相比非凍結(jié)條件,凍結(jié)條件下滲達(dá)到穩(wěn)定時(shí)LNAPL相的質(zhì)量通量偏差為-51.1%.凍結(jié)條件LNAPL溶解相在水平方向的遷移速度大于非凍結(jié)條件,凍結(jié)條件下如果不考慮凍結(jié),下滲過程中LNAPL溶解量與揮發(fā)量最大偏差分別為2.4%、-24.9%.而凍結(jié)條件對(duì)LNAPL下滲達(dá)到穩(wěn)定時(shí)的溶解量、揮發(fā)量影響很小,僅揮發(fā)量產(chǎn)生1.8%的偏差.但會(huì)減少LNAPL的吸附量,其偏差為-6.6%.

輕非水相流體(LNAPL)；多相多場(chǎng)耦合；凍結(jié)勢(shì)；溶質(zhì)運(yùn)移

隨著石油化工業(yè)的發(fā)展,非水相有機(jī)污染物在其煉制、貯存及運(yùn)輸過程中泄漏至環(huán)境,對(duì)人類健康造成嚴(yán)重威脅.建立凍結(jié)條件下飽和非飽和帶的LNAPL遷移模型有利于預(yù)測(cè)其在寒區(qū)泄漏后于土壤-地下水中的遷移分布,為L(zhǎng)NAPL污染防控提供技術(shù)支撐[1].

NAPL根據(jù)密度大小,可分為輕非水相(Light NAPL)和重非水相(Dense NAPL),其控制方程基于多相流理論構(gòu)建.數(shù)值模擬作為一種預(yù)測(cè)NAPL在土壤—地下水中遷移分布的重要手段,在NAPL污染防治及油氣開采中都得到了廣泛應(yīng)用.目前關(guān)于三相流在飽和-非飽和帶中的模擬已經(jīng)相對(duì)成熟.有關(guān)在飽和帶中不同降雨速率、介質(zhì)非均質(zhì)性及NAPL泄漏速率、地下水流速下的NAPL兩相流模型已比較完善[2-10].在非飽和帶,相關(guān)研究也對(duì)NAPL遷移過程進(jìn)行了數(shù)值模擬[11-14].而針對(duì)NAPL從非飽和帶向飽和帶遷移的過程,Lenhard等通過建立自由相、殘余相和截留相NAPL與流體水頭的函數(shù)關(guān)系式建立了遷移模型,預(yù)測(cè)了其分布[15-16].

迄今,寒區(qū)凍結(jié)條件下LNAPL的遷移研究仍相對(duì)較少,研究表明LNAPL在凍結(jié)條件下泄漏后在土壤包氣帶及地下水中的遷移因?yàn)楸в绊憰?huì)與非凍結(jié)條件下的遷移圖像有所不同,冰晶有效阻止石油垂直遷移,加強(qiáng)石油沿優(yōu)先流路徑下滲及在地表淺層水平流動(dòng),但即使完全冰飽和的土也表現(xiàn)出一定的冰滲透性[17-19].且在凍結(jié)條件下,土壤包氣帶伴隨凍結(jié)滯水的形成,飽和水汽壓作用減小,而毛管作用得以加強(qiáng)[21],地下水中神經(jīng)節(jié)狀Ganglia在凍結(jié)勢(shì)和基質(zhì)勢(shì)作用下會(huì)向凍結(jié)峰積聚,且有部分LNAPL匯聚于表層[20,22].國(guó)內(nèi)也有少數(shù)學(xué)者通過室內(nèi)凍結(jié)箱實(shí)驗(yàn)開展LNAPL于凍結(jié)條件下的遷移規(guī)律研究[22-23].然而由于凍結(jié)固相(冰)的存在,需將溫度場(chǎng)與LNAPL遷移所涉及的滲流場(chǎng)、化學(xué)場(chǎng)耦合,同時(shí)多相多場(chǎng)耦合問題求解往往不易收斂,有關(guān)凍結(jié)條件下氣相、水相、NAPL相三相運(yùn)移模擬尚未見相關(guān)報(bào)道.

COMSOL Multiphysics軟件常用于多場(chǎng)多相問題的求解.而凍結(jié)條件下LNAPL的遷移,涉及到水、氣、NAPL三相在多孔介質(zhì)中運(yùn)移,是溫度場(chǎng)、流場(chǎng)及化學(xué)場(chǎng)多場(chǎng)耦合的問題.已有的COMSOL模塊設(shè)置中缺少凍結(jié)模塊,無法直接用于建立凍結(jié)條件下LNAPL遷移模型.本研究通過在COMSOL中增設(shè)凍結(jié)過程控制方程對(duì)凍結(jié)條件下LNAPL入滲過程進(jìn)行模擬,建立固相-氣相-水相-NAPL相四相共存,溫度場(chǎng)、流場(chǎng)、化學(xué)場(chǎng)三場(chǎng)耦合模型,模擬凍結(jié)條件下LNAPL在飽和-非飽和帶的遷移行為及分布形態(tài),旨在克服LNAPL遷移模擬模型中未考慮冰相而無法適用于寒區(qū)的缺陷.

1 凍結(jié)條件下多孔介質(zhì)中LNAPL遷移的數(shù)學(xué)模型

1.1 數(shù)學(xué)模型的建立

凍結(jié)條件下多孔介質(zhì)LNAPL遷移因有凍結(jié)固相(冰)存在,需要在以往三相流運(yùn)移方程的水相運(yùn)移質(zhì)量守恒方程中加入冰相,并通過未凍水和凍結(jié)冰的經(jīng)驗(yàn)關(guān)系式實(shí)現(xiàn)溫度場(chǎng)和流場(chǎng)的耦合.

凍結(jié)條件下LNAPL的遷移包括兩個(gè)過程,一是LNAPL自由相在水-氣-LNAPL-冰四相共存時(shí)的毛細(xì)壓力作用下遷移,二是溶于水的LNAPL及其揮發(fā)在空氣中的LNAPL的溶質(zhì)遷移過程,以對(duì)流-彌散作用為主.該過程涉及水-冰相變、水-氣-NAPL三相流遷移、多孔介質(zhì)溫度傳導(dǎo)及溶質(zhì)運(yùn)移模型耦合,其控制方程如下:

1.1.1 水-氣-NAPL-固(冰)四相控制方程 據(jù)各相的質(zhì)量守恒定律[16,24]得到:

自由相LNAPL在水中的溶解過程用式(4)表達(dá)[26]:

自由相LNAPL揮發(fā)于空氣可用式(5)表達(dá)[26]:

溶于水的LNAPL揮發(fā)于空氣可用式(6)表達(dá)[26]:

溶于水中的LNAPL吸附于土壤骨架可用式(7)表示[27]:

式中:oc為土體中有機(jī)碳的質(zhì)量分?jǐn)?shù);oc為有機(jī)碳分配系數(shù).

而對(duì)于毛管壓力CNW、CGN以及相對(duì)滲透率rW、rN、rG,需構(gòu)建--關(guān)系式.本模型-關(guān)系將van Genuchten方程拓展到三相流[28],得到:

而對(duì)于飽和滲透率,它隨凍結(jié)過程孔隙度的變化而變化(將冰視為土壤固體骨架來實(shí)現(xiàn)水-氣-NAPL-固四相向水-氣-NAPL三相的轉(zhuǎn)換),采用Kozeny-Carman模型[29],表達(dá)式如下:

式中:為顆粒直徑.

1.1.2 多孔介質(zhì)傳熱 多孔介質(zhì)傳熱方程即熱量守恒方程如下:

式中:為流體密度;p為比定壓熱容;為達(dá)西流速;為溫度;eff為熱導(dǎo)率;f為融化潛熱[24].

為了實(shí)現(xiàn)溫度場(chǎng)與三相流場(chǎng)的耦合,一方面需要引入凍結(jié)條件下未凍水和凍結(jié)冰的經(jīng)驗(yàn)比例關(guān)系式[30](式16),另一方面,通過將冰視為土壤固體骨架來實(shí)現(xiàn)水-氣-NAPL-固四相向水-氣-NAPL三相的轉(zhuǎn)換,即凍結(jié)過程中固體骨架的孔隙度將隨冰的形成而不斷變化,且傳熱方程中的比熱容和導(dǎo)熱系數(shù)是水-氣-NAPL-固四相體積占比的加權(quán)函數(shù).

凍結(jié)場(chǎng)和三相流場(chǎng)的耦合關(guān)鍵即未凍水與凍結(jié)冰的經(jīng)驗(yàn)比例式如下:

式中:()為比例系數(shù);f為起始凍結(jié)溫度(絕對(duì)溫度);系數(shù)1.1為水與冰的密度之比.

1.1.3 溶解相-揮發(fā)相LNAPL質(zhì)量平衡方程 LNAPL溶解相的溶質(zhì)運(yùn)移方程為:

LNAPL揮發(fā)相的溶質(zhì)運(yùn)移方程為:

式中:為L(zhǎng)NAPL的物質(zhì)的量質(zhì)量分?jǐn)?shù);W、G為水動(dòng)力彌散系數(shù).

非凍結(jié)條件下LNAPL遷移的控制方程除水相運(yùn)移方程(式1)不包含含冰項(xiàng)及無需進(jìn)行多孔介質(zhì)溫度傳導(dǎo)模塊設(shè)置外,其余控制方程與凍結(jié)條件下LNAPL遷移控制方程一致.

1.2 模型求解

采用COMSOL Multiphysics的多孔介質(zhì)多相流模塊、多孔介質(zhì)溫度傳導(dǎo)模塊、多孔介質(zhì)物質(zhì)傳遞模塊及達(dá)西定律模塊輸入并求解凍結(jié)條件下LNAPL遷移控制方程.COMSOL Multiphysics是一款強(qiáng)大的多物理場(chǎng)耦合軟件,常用于多場(chǎng)多相問題的求解.而凍結(jié)條件下LNAPL的遷移,涉及到水、氣、NAPL三相在多孔介質(zhì)中運(yùn)移,是溫度場(chǎng)、流場(chǎng)及化學(xué)場(chǎng)高度耦合的問題.COMSOL中水相方程并未包含冰相,故在模擬凍結(jié)時(shí)需要在水相質(zhì)量守恒方程中加入冰相隨時(shí)間變化(式1).而氣相、LNAPL相則采用其內(nèi)置方程(式2～3).COMSOL內(nèi)置的多孔介質(zhì)相傳遞模塊可用來模擬三相流運(yùn)移,但其未考慮四相質(zhì)量轉(zhuǎn)換過程,需要在模塊中增設(shè)(式4～8).毛管壓力以及相對(duì)滲透率也需要在COMSOL中自行構(gòu)建其與相對(duì)飽和度間的關(guān)系式(式9～13).多孔介質(zhì)溫度傳導(dǎo)模塊可用來模擬介質(zhì)傳熱,但其熱傳導(dǎo)方程缺少冰的融化潛熱項(xiàng),需要在源項(xiàng)增設(shè)(式15).而且其缺少凍結(jié)方程設(shè)置模塊,需將凍結(jié)條件下未凍結(jié)水和凍結(jié)冰的經(jīng)驗(yàn)比例關(guān)系式(式16)在COMSOL中進(jìn)行設(shè)置.

對(duì)于非凍條件下LNAPL遷移模型,可采用COMSOL的多孔介質(zhì)多相流模塊、多孔介質(zhì)稀物質(zhì)傳遞模塊及達(dá)西定律模塊輸入并求解LNAPL遷移控制方程.同樣將四相質(zhì)量轉(zhuǎn)換過程在多孔介質(zhì)相傳遞模塊中進(jìn)行增設(shè).而三相流運(yùn)移控制方程均采用COMSOL內(nèi)置方程.

COMSOL通過有限單元法離散上述模型,以求解凍結(jié)、非凍結(jié)條件下LNAPL遷移的非穩(wěn)定問題.凍結(jié)、非凍結(jié)條件下LNAPL遷移模擬輸出結(jié)果中自由相LNAPL的含量用飽和度表示,溶解相及揮發(fā)相LNAPL含量用質(zhì)量濃度表示.

2 算例驗(yàn)證

2.1 算例概述及初始邊界條件

算例改自文獻(xiàn)[26]所建立的數(shù)值模擬算例,文獻(xiàn)[26]中的算例研究了LNAPL定點(diǎn)泄漏后于土壤多孔介質(zhì)中的遷移行為.本文采用與該算例相同的水文地質(zhì)條件和污染物泄漏條件建立凍結(jié)條件下LNAPL遷移模型.為驗(yàn)證本次所構(gòu)建凍結(jié)條件LNAPL遷移模型的可靠性,首先利用所建立的模型開展非凍結(jié)條件下模擬,并與文獻(xiàn)[26]中結(jié)果進(jìn)行對(duì)比驗(yàn)證.進(jìn)而在非凍條件下LNAPL遷移模型中加入冰相、溫度場(chǎng)和流場(chǎng)的耦合方程及凍結(jié)相關(guān)邊界條件進(jìn)行了凍結(jié)條件下LNAPL在多孔介質(zhì)中的二維遷移模擬.

建立的凍結(jié)條件下LNAPL在非飽和帶中二維遷移數(shù)值算例見圖1.研究區(qū)為二維剖面,長(zhǎng)1.5m,厚1.5m,下部0.5m為地下水飽和區(qū),取地下水面為基準(zhǔn)面.區(qū)域網(wǎng)格剖分采用COMSOL的物理場(chǎng)細(xì)化網(wǎng)格剖分.

至于初始條件,對(duì)于非飽和帶與飽和帶分別進(jìn)行如下設(shè)置:

模型的邊界條件如下設(shè)置:

(1)流場(chǎng)

(2)化學(xué)場(chǎng)

(3)溫度場(chǎng)

初始邊界條件如圖1.對(duì)于非凍結(jié)條件下,無需設(shè)置溫度場(chǎng)邊界及初始條件,流場(chǎng)和溶質(zhì)場(chǎng)邊界條件及初始條件與凍結(jié)條件下保持一致.

2.2 基本參數(shù)

模型中采用的參數(shù)值如表1,包括土壤、水、LNAPL、冰的特征參數(shù)以及三相質(zhì)量轉(zhuǎn)換的相關(guān)參數(shù).LNAPL在此以甲苯為模擬對(duì)象.

表1 模型相關(guān)參數(shù)

2.3 情景設(shè)置

設(shè)置凍結(jié)條件與非凍結(jié)條件LNAPL下滲2種情景進(jìn)行模擬,通過模擬結(jié)果對(duì)比來探究?jī)鼋Y(jié)條件對(duì)LNAPL遷移影響.非凍結(jié)條件下除溫度場(chǎng)外,其它條件與凍結(jié)條件下保持一致.

3 結(jié)果與討論

3.1 溫度場(chǎng)

凍結(jié)過程中溫度隨凍結(jié)時(shí)間增長(zhǎng)在地表淺層由0℃逐漸降低至-14℃(圖2),冰飽和度在地表淺層逐漸增加,凍結(jié)12h時(shí),地表淺層冰飽和度達(dá)到0.8(圖3).因?yàn)闇囟戎辉谏喜勘影l(fā)生變化,下部恒為5℃,故取土壤表層0.1m厚的溫度場(chǎng)模擬結(jié)果進(jìn)行分析.可知凍結(jié)場(chǎng)在很短時(shí)間內(nèi)達(dá)到穩(wěn)定,溫度在土壤表層驟降,在土壤表層會(huì)形成薄冰層.這再現(xiàn)了低溫條件下多孔介質(zhì)LNAPL、水遷移實(shí)驗(yàn)所得實(shí)驗(yàn)結(jié)果:低溫環(huán)境中LNAPL相遷移時(shí)高溫區(qū)域面積較大,等溫線在土壤表層較為密集,溫度梯度變化較快[23].說明LNAPL泄漏會(huì)對(duì)土壤的對(duì)流傳熱產(chǎn)生影響,對(duì)土壤起保溫作用.LNAPL流經(jīng)薄冰層時(shí),順著薄冰層做水平方向流動(dòng).而當(dāng)流體做水平方向流動(dòng)時(shí),在垂直于流動(dòng)方向的熱量傳遞,主要以熱傳導(dǎo)方式進(jìn)行.LNAPL流體的導(dǎo)熱系數(shù)較小,傳熱熱阻主要集中在該層中,溫差也集中在該層[31].且達(dá)西數(shù)較小時(shí),也會(huì)明顯削弱對(duì)流傳熱[32].

在近地表凍層處取坐標(biāo)為(0.75m,0.98m)的點(diǎn)A(圖1)為代表性點(diǎn)進(jìn)行進(jìn)一步分析,該點(diǎn)處空隙里的水組分和冰含量隨溫度變化的情況如圖4(a)所示.隨溫度降低至凍結(jié)溫度0.3℃時(shí),冰飽和度由0逐漸增加,而水飽和度隨溫度降低逐漸減小.溫度達(dá)到-10℃時(shí),水飽和度僅為0.17,而冰飽和度達(dá)到0.80.而該點(diǎn)處的水相、冰相隨時(shí)間的變化及由水冰相變?cè)斐傻臐摕嶙兓鐖D4(b)所示.水飽和度在未凍結(jié)前期隨LNAPL下滲水壓增大而增加,后趨于穩(wěn)定(圖4(c)).而隨溫度下降至0攝氏度時(shí)(圖4(a)),即在第2h時(shí)冰飽和度由0逐漸增加,同時(shí)水飽和度減小,水相轉(zhuǎn)化為冰相.而潛熱伴隨這一過程釋放,達(dá)到峰值后下降,隨凍結(jié)速度變慢而降低,當(dāng)?shù)竭_(dá)12h時(shí),凍結(jié)達(dá)到穩(wěn)定,潛熱趨于零(圖4(b)).

圖2 溫度場(chǎng)隨時(shí)間的演變過程

圖3 含冰率隨時(shí)間的演變過程

3.2 流場(chǎng)

3.2.1 凍結(jié)條件下LNAPL遷移特征 如圖5(a)所示,凍結(jié)前5h LNAPL在重力作用下向下遷移,在含冰層下端有突變界面,LNAPL含量在泄漏點(diǎn)的凍層處聚集.對(duì)凍結(jié)層點(diǎn)A處(圖1)的孔隙率和滲透率進(jìn)行分析,可知當(dāng)溫度隨時(shí)間增長(zhǎng)降至凍結(jié)溫度時(shí),由于水相轉(zhuǎn)換為冰,介質(zhì)孔隙度變小,而滲透率隨孔隙度的減小而降低(圖6),即凍層對(duì)LNAPL的下滲有明顯阻擋作用.而隨著凍結(jié)作用增強(qiáng),LNAPL不再繼續(xù)向下遷移.凍結(jié)5h后,其在凍結(jié)勢(shì)的作用下克服重力勢(shì)作用向上遷移[21],遷移程度較小,隨后達(dá)到穩(wěn)定.而對(duì)于包氣帶中的水分,其在凍結(jié)前1h內(nèi)向上運(yùn)移(圖5(b)),由圖5(c)壓力分布及流線圖,地表土層處水壓逐漸升高,流線向上,可知趨使水分向上運(yùn)移的是隨LNAPL下滲產(chǎn)生的壓力變化.隨后其飽和度隨LNAPL污染暈的移動(dòng)而變化,含冰層處水分轉(zhuǎn)化為冰,水飽和度相對(duì)較低(圖5(b)).包氣帶中的氣體則隨LNAPL的下滲逐漸向邊緣擴(kuò)散,直至整個(gè)區(qū)域被LNAPL和水分填充.且其在含冰層下部出現(xiàn)模擬高值,發(fā)生聚集(圖5(d)).

如圖7所示,由凍結(jié)5h和凍結(jié)12h的量分布剖面圖可知,LNAPL向下泄漏,其量由泄漏表面至深部逐漸減少.而且LNAPL從凍結(jié)5h至凍結(jié)12h量分布曲線上移,剖面上LNAPL量均有不同程度減小,在凍結(jié)勢(shì)作用下向地表方向遷移,暴露在表層.以凍結(jié)層濃度減小最明顯,最高可達(dá)5mg/L.而在凍結(jié)層下部量減少較小,這是因?yàn)橄虏縇NAPL向上遷移,會(huì)受到凍層阻滯,形成累積,這與李興柏等[22]對(duì)凍結(jié)條件下的石油遷移試驗(yàn)結(jié)果一致.故下部遷移的LNAPL累積量抵消了部分LNAPL的向上遷移量.由模擬結(jié)果計(jì)算得凍結(jié)條件下前5h LNAPL向下遷移量為21.24mL,而凍結(jié)5h～凍結(jié)12h LNAPL向上遷移至土壤表層的量為1.04mL.

圖6 滲透率及孔隙度隨時(shí)間變化

圖7 LNAPL運(yùn)移量的垂向分布

3.2.2 非凍結(jié)條件下LNAPL遷移特征 如圖8所示,LNAPL隨時(shí)間推移逐漸向下滲透,前5h變化程度大,隨后達(dá)到區(qū)域飽和穩(wěn)定狀態(tài).非凍結(jié)條件下LNAPL運(yùn)移5和12h相比的量分布差異以及向下泄漏量分布剖面變化如圖7所示,可知LNAPL在非凍結(jié)條件下在重力勢(shì)作用下向下滲,其量同樣由泄漏表面至深部逐漸減少.在第5h～第12h,上部土壤水帶和中間帶LNAPL濃度略有減少,LNAPL向下遷移,在下部地下水飽和帶以及毛細(xì)水帶累積,而且泄露源對(duì)上部土壤水帶和中間帶的補(bǔ)給量小于其向地下水飽和帶及毛細(xì)水帶的下滲量.由模擬結(jié)果計(jì)算得非凍結(jié)條件下前5h LNAPL向下遷移量為38.59mL,至第12h達(dá)到41.30mL.

非凍結(jié)條件與凍結(jié)條件下LNAPL下滲模擬結(jié)果相比,因其不存在凍結(jié)層的阻擋作用,下滲速度快,下滲量遠(yuǎn)大于凍結(jié)時(shí)下滲量.12h后地下水面處體積分?jǐn)?shù)可達(dá)0.040(圖8),而凍結(jié)條件下地下水面土壤處LNAPL體積分?jǐn)?shù)僅達(dá)0.016(圖5(a)).非凍結(jié)條件LNAPL向地下水面遷移,不存在向上遷移現(xiàn)象,而凍結(jié)條件下LNAPL在凍結(jié)5h后在凍結(jié)勢(shì)作用下發(fā)生向上遷移(圖7).凍結(jié)條件下LNAPL相的下滲量相比非凍結(jié)條件下的量產(chǎn)生-51.1%偏差.

圖8 常溫條件下LNAPL飽和度變化

3.3 化學(xué)場(chǎng)

3.3.1 凍結(jié)條件LNAPL溶解揮發(fā)相遷移特征 如圖9(a)所示,隨LNAPL下滲,其在水中的溶解度逐漸增加,隨后趨于穩(wěn)定狀態(tài).凍結(jié)1h時(shí),LNAPL溶解相分布圖像呈橢圓型,在地下水中濃度達(dá)116.5mg/L.遷移1h時(shí),LNAPL揮發(fā)相在包氣帶中濃度可達(dá)118.7mg/L.凍結(jié)穩(wěn)定后,在LNAPL以0.138cm/min的泄漏速度下,其在包氣帶中濃度可達(dá)186.7mg/L (LNAPL在水中的溶解度為515mg/L).而凍結(jié)條件下LNAPL揮發(fā)相隨氣相逐漸向兩邊擴(kuò)散(圖9(b)),在土壤表層濃度可達(dá)101.0mg/L(LNAPL在空氣中的溶解度為139mg/L).

模擬可得凍結(jié)條件下不同狀態(tài)LNAPL含量關(guān)系圖,如圖10所示.可知凍結(jié)條件下凍結(jié)1h時(shí)總?cè)芙饬繛?.304mL,揮發(fā)量為0.195mL.凍結(jié)穩(wěn)定時(shí)LNAPL總下滲量為23.74mL,溶于水的量為0.48mL,氣相LNAPL含量為0.165mL,而LNAPL吸附量為2.71mL(圖10(a)).凍結(jié)條件下12h時(shí)可自由移動(dòng)LNAPL量(溶解相LNAPL與自由相LNAPL的和)為20.68mL,其含量在凍結(jié)前5h逐漸增加,從凍結(jié)5h后的21.72mL逐漸減少,減少量為凍結(jié)勢(shì)作用下向土壤表層遷移量.而吸附量隨凍結(jié)時(shí)間增長(zhǎng)增加到2.71mL后保持穩(wěn)定(圖10(b)).

3.3.2 非凍結(jié)條件LNAPL溶解揮發(fā)相遷移特征 如圖11(a)所示,水中的溶解相隨LNAPL下滲逐漸增加,隨后穩(wěn)定.遷移1h時(shí),LNAPL溶解相分布圖像呈楔形,其在垂直方向擴(kuò)散速度大于水平方向,在地下水飽和區(qū)濃度可達(dá)117.2mg/L.遷移1h時(shí),LNAPL揮發(fā)相在包氣帶中濃度可達(dá)143.2mg/L.下滲穩(wěn)定后,在LNAPL以0.138cm/min的泄漏速度下,LNAPL溶解相在包氣帶中濃度可達(dá)185.8mg/L(LNAPL在水中的溶解度為515mg/L).而非凍結(jié)條件下LNAPL揮發(fā)相隨氣相逐漸向兩邊擴(kuò)散(圖11(b)),在土壤表層濃度可達(dá)101.9mg/L(LNAPL在空氣中的溶解度為139mg/L).

圖9 凍結(jié)條件下LNAPL溶解相和揮發(fā)相濃度分布

圖10 凍結(jié)條件不同狀態(tài)LNAPL含量變化

非凍結(jié)下不同狀態(tài)LNAPL含量關(guān)系圖如圖12所示.溶解相LNAPL、揮發(fā)相LNAPL及吸附相LNAPL隨凍結(jié)時(shí)間增長(zhǎng)逐漸增多,達(dá)到穩(wěn)定時(shí)含量吸附相LNAPL遠(yuǎn)大于溶解相,而揮發(fā)相最少(圖12(a)).可知非凍結(jié)條件下下滲1h時(shí)總?cè)芙饬繛?.196mL,揮發(fā)量為0.251mL.穩(wěn)定后LNAPL總下滲量為41.30mL,溶于水的量為0.48mL,氣相LNAPL含量為0.168mL,LNAPL吸附量為2.90mL(圖12(a)).吸附量占泄漏總量7%,吸附量大于溶解量,這也與Kim和Corapcioglu 的LNAPL泄漏實(shí)驗(yàn)結(jié)果吻合[33].非凍結(jié)條件LNAPL從泄漏點(diǎn)開始不斷下滲,至第7h達(dá)到38.43mL穩(wěn)定狀態(tài),吸附量隨凍結(jié)時(shí)間增長(zhǎng)增加到2.90mL后保持穩(wěn)定(圖12(b)).而可移動(dòng)的LNAPL相(包括溶解相和自由相)略有減少后趨于穩(wěn)定.

圖11 常溫條件下LNAPL溶解相和揮發(fā)相濃度分布

圖12 常溫條件不同狀態(tài)LNAPL含量變化

對(duì)比LNAPL于凍結(jié)、非凍結(jié)條件下1h時(shí)溶解相、揮發(fā)相濃度分布可知,LNAPL溶解相在凍結(jié)條件下在水平方向上擴(kuò)散速度大于非凍結(jié)條件,地下水飽和區(qū)溶解相濃度差為0.7mg/L,土壤非飽和帶中LNAPL揮發(fā)相最大濃度之差為24.6mg/L.在下滲過程中LNAPL溶解量與揮發(fā)量?jī)鼋Y(jié)條件相比非凍結(jié)條件會(huì)產(chǎn)生2.4%、-24.9%的最大偏差.穩(wěn)定時(shí)LNAPL在水中及空氣中的溶解擴(kuò)散量濃度差別在1mg/L內(nèi),其受冰凍層影響小.又由不同狀態(tài)LNAPL含量變化圖可知兩種情景下LNAPL溶于水的量相等,而氣相LNAPL含量差僅為0.003mL,偏差為-1.8%.即含冰層對(duì)LNAPL的溶解量和揮發(fā)量造成的影響很小,這是因?yàn)閮鼋Y(jié)條件下溶解量和揮發(fā)量主要通過分子擴(kuò)散作用通過冰晶孔隙,含冰層的存在對(duì)分子擴(kuò)散作用影響小.而對(duì)比吸附量,非凍結(jié)條件LNAPL吸附量高于凍結(jié)條件LNAPL吸附量,偏差為-6.6%.這是因?yàn)橐环矫鍸NAPL在非凍結(jié)條件下下滲量大于凍結(jié)條件,另一方面土壤表面被冰包裹會(huì)減少土壤的比表面積從而降低其吸附性能.

3.4 模型驗(yàn)證

所建模型未加入凍結(jié)模塊時(shí),與文獻(xiàn)[26]建立的水文地質(zhì)條件及污染物泄漏條件設(shè)置一致.故可用文獻(xiàn)[26]的模擬結(jié)果與本模型模擬結(jié)果進(jìn)行對(duì)比以驗(yàn)證所建模型可靠性.由圖13可以看出,兩者模擬結(jié)果一致,驗(yàn)證了所建模型在非凍結(jié)條件下的正確性.

圖13 模型與Kacem案例不同相態(tài)LNAPL含量變化對(duì)比

所建模型在凍結(jié)條件下的可靠性,由于相關(guān)研究較少,難以找到相關(guān)實(shí)驗(yàn)數(shù)據(jù)直接驗(yàn)證.后續(xù)可開展凍結(jié)條件下多相流遷移實(shí)驗(yàn)來進(jìn)一步驗(yàn)證模型.

4 結(jié)論

4.1 本文在COMSOL Multiphysics的多孔介質(zhì)相傳遞、多孔介質(zhì)溫度傳導(dǎo)、多孔介質(zhì)物質(zhì)傳遞模塊的基礎(chǔ)上,通過增設(shè)相關(guān)模塊構(gòu)建了包含水、氣、LNAPL、冰四相的溫度場(chǎng)、流場(chǎng)及化學(xué)場(chǎng)三場(chǎng)耦合模型.可用于模擬及定量刻化LNAPL泄漏過程中發(fā)生凍結(jié)后LNAPL的遷移分布,為寒區(qū)LNAPL泄漏后的修復(fù)提供參考.

4.2 LNAPL在地表淺層的泄漏會(huì)影響土壤與環(huán)境的熱量傳導(dǎo)過程.在靠近薄冰層處總有一薄層流體順著薄冰層做水平流動(dòng),而冰層介質(zhì)的達(dá)西數(shù)較小,會(huì)明顯削弱土壤的對(duì)流傳熱,形成熱阻層,使得溫度驟變界面在薄冰層產(chǎn)生.

4.3 凍結(jié)條件下LNAPL遷移及水分運(yùn)移受重力勢(shì)、基質(zhì)勢(shì)及凍結(jié)勢(shì)的控制,LNAPL先在重力勢(shì)下向下遷移,隨凍結(jié)時(shí)間增長(zhǎng),凍結(jié)勢(shì)和重力勢(shì)平衡后在凍結(jié)勢(shì)作用下發(fā)生一定程度的向上遷移.凍結(jié)條件下LNAPL相的下滲量相比非凍結(jié)條件下的量產(chǎn)生-51.1%的偏差.

4.4 凍結(jié)條件LNAPL溶解相在水平方向遷移速度大于非凍結(jié)條件,下滲過程中凍結(jié)條件LNAPL溶解量與揮發(fā)量相比非凍結(jié)條件會(huì)產(chǎn)生2.4%、-24.9%的最大偏差,下滲達(dá)到穩(wěn)定時(shí)僅揮發(fā)量產(chǎn)生1.8%的偏差.溶解量和揮發(fā)量主要通過分子擴(kuò)散作用通過冰晶孔隙,含冰層的存在對(duì)分子擴(kuò)散作用影響小.而含冰層的存在對(duì)吸附量存在影響,凍結(jié)條件下LNAPL吸附量略低于非凍結(jié)條件下LNAPL吸附量,偏差為-6.6%.

[1] 李國(guó)玉,馬 巍,李興柏,等.多年凍土區(qū)石油污染物遷移過程研究回顧與展望 [J]. 冰川凍土, 2011,33(4):947-952.

Li G Y, Ma W, Li X B, et al. Migration of the petroleum pollutants in permafrost regions: Review and prospect [J]. Journal of Glaciology and Geocryology, 2011,33(4):947-952.

[2] 施小清,吳吉春,劉德朋,等.飽和介質(zhì)中重非水相液體運(yùn)移的數(shù)值模擬及敏感性分析 [J]. 南京大學(xué)學(xué)報(bào)(自然科學(xué)版), 2011,47(3):299- 307.

Shi X Q, Wu J C, Liu D P, et al. Numerical simulation of transportation of dense nonaqueous phase liquids in the subsurface environment [J]. Journal of Nanjing University (Natural Sciences), 2011,47(3):299-307.

[3] Yang Z B, Zandin H, Niemi A, et al. The role of geological heterogeneity and variability in water infiltration on non-aqueous phase liquid migration [J]. Environmental Earth Sciences, 2013,68(7): 2085-2097.

[4] 王 涵,盧文喜,李久輝,等.地下水DNAPLs污染多相流的隨機(jī)模擬及其不確定性分析 [J]. 中國(guó)環(huán)境科學(xué), 2018,38(7):2572-2579.

Wang H, Lu W X, Li J H, et al. Stochastic simulation and uncertainty analysis of multi-phase flow of groundwater polluted by DNAPLs [J]. China Environmental Science, 2018,38(7):2572-2579.

[5] 鄭 菲,高燕維,施小清,等.地下水流速及介質(zhì)非均質(zhì)性對(duì)重非水相流體運(yùn)移的影響 [J]. 水利學(xué)報(bào), 2015,46(8):925-933.

Zheng F, Gao Y W, Shi X Q, et al. Influence of groundwater flow velocity and geological heterogeneity on DNAPL migration in saturated porous media [J]. Journal of Hydraulic Engineering, 2015,46 (8):925-933.

[6] Van Genuchten M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils [J]. Soil Science Society of America Journal, 1980,44(5):892-898.

[7] Brooks R H, Corey A T. Properties of porous media affecting fluid flow [J]. Journal of the Irrigation and Drainage Division, 1966,92(2): 61–90.

[8] Falta R W, Pruess K, Javandel I, et al. Numerical modeling of steam injection for the removal of nonaqueous phase liquids from the subsurface: 1. Numerical formulation [J]. Water Resource Research, 1992,28(2):433–449.

[9] Adenekan A E, Patzek T W, Pruess K. Modeling of multi-phase transport of multicomponent organic contaminants and heat in the subsurface: Numerical model formulation [J]. Water Resource Research, 1993,29(11):3727–3740.

[10] White M D, Oostrom M, Lenhard R J. Modeling fluid flow and transport in variably saturated porous media with the STOMP simulator. 1. Nonvolatile three-phase model description [J]. Advances in Water Resources, 1995,18(6):353–364.

[11] 阮冬梅,卞建民,王 倩,等.低滲透介質(zhì)中輕非水相流體遷移轉(zhuǎn)化規(guī)律 [J]. 中國(guó)環(huán)境科學(xué), 2021,41(4):1815-1823.

Ruan D M, Bian J M, Wang Q, et al. The migration and transformation of light non-aqueous fluid in silty clay [J]. China Environmental Science, 2021,41(4):1815-1823.

[12] Hofstee C, Lenhard R J, et al. Flow behavior and residual saturation formation of liquid carbon tetrachloride in unsaturated heterogeneous porous media [J]. Journal of Contaminant Hydrology, 2003,64(1):93-112.

[13] Lo I M, Hu L M, Meegoda J N. Centrifuge modeling of light nonaqueous phase liquids transport in unsaturated soils [J]. Journal of Geotechnical and Geoenvironmental Engineering, 2004,130(5):535-539.

[14] Kechavarzi C, Soga K, Lllangasekare T H. Two-dimensional laboratory simulation of LNAPL infiltration and redistribution in the vadose zone [J]. Journal of Contaminant Hydrology, 2005,76(3/4): 211–233.

[15] Lenhard R J, Rayner J L, Davis G B. A practical tool for estimating subsurface LNAPL distributions and transmissivity using current and historical fluid levels in ground water wells: Effects of entrapped and residual LNAPL [J]. Journal of Contaminant Hydrology, 2017,205: 1–11.

[16] Kacem M, Esrael D, Boeije C S, et al. Multiphase Flow Model for NAPL Infiltration in Both the Unsaturated and Saturated Zones [J]. Journal of Environmental Engineering, 2019,145(11):04019072.

[17] Biggar K W. Site investigations of fuel spill migration into permafrost [J]. Journal of Cold Regions Engineering, 1998,12(2):84–104.

[18] Chuvilin E M, Naletova N S, Miklyaeva E C, et al. Factors affecting spreadability and transportation of oil in regions of frozen ground [J]. Polar Record, 2001,37(202):229-238.

[19] Barnes D L, Wolfe S M. Influence of Ice on the Infiltration of Petroleum into Frozen Coarse-grained Soil [J]. Petroleum Science and Technology, 2008,26(7/8):856-867.

[20] Singh K, Niven R K. Non-aqueous Phase Liquid Spills in Freezing and Thawing Soils: Critical Analysis of Pore-Scale Processes [J]. Critical Reviews in Environmental Science and Technology, 2013, 43(6):551-597.

[21] 徐樹林,那平山,武俊英.關(guān)于凍結(jié)滯水的探討 [J]. 內(nèi)蒙古林學(xué)院報(bào), 1994,(2):46-51.

Xu S L, Na P S, Wu J Y. An approach to frozen perched groundwater [J]. Journal of inner Mongolia forestry college, 1994,(2):46-51.

[22] 李興柏,李國(guó)玉.溫度梯度對(duì)多年凍土區(qū)石油遷移影響的研究 [J]. 甘肅科學(xué)學(xué)報(bào), 2013,25(1):73-76.

Li X B, Li G Y. Experimental study on the influence of temperature gradient on the migrating process of oil contaminants in permafrost regions [J]. Journal of Gansu Sciences, 2013,25(1):73-76.

[23] 吳國(guó)忠,趙文浩,呂 妍,等.多孔介質(zhì)內(nèi)油、水介質(zhì)遷移對(duì)比低溫實(shí)驗(yàn)[J]. 當(dāng)代化工, 2016,45(6):1108-1111.

Wu G Z, Zhao W H, Lv Y, et al. Comparison of the migration of oil and water in porous medium with cryogenic experiment [J]. Contemporary Chemical Industry, 2016,45(6):1108-1111.

[24] Zhang M L, Wen Z, Xue K, et al. A coupled model for liquid water, water vapor and heat transport of saturated–unsaturated soil in cold regions: model formulation and verification [J]. Environmental Earth Sciences, 2016,75(8):1-19.

[25] 胡黎明,邢巍巍,吳照群.多孔介質(zhì)中非水相流體運(yùn)移的數(shù)值模擬 [J]. 巖土力學(xué), 2007,28(5):951-955.

Hu N M, Xin W W, Wu Z Q. Numerical simulation of non-aqueous phase liquids migration in porous media [J]. Rock and Soil Mechanics, 2007,28(5):951-955.

[26] Kacem M, Benadda B. Mathematical model for multiphase extraction simulation [J]. Journal of Environmental Engineering, 2018,144(6): 04018040.

[27] 高彥斌,張松波,李 韜,等.飽和黏性土中重質(zhì)非水相有機(jī)污染物縱向遷移數(shù)值模擬[J]. 同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版), 2020,48(1):24- 32.

Gao Y B, Zhang S B, Li T, et al. Numerical analysis of vertical migration of dense nonaqueous-phase liquids in saturated clay [J]. Journal of Tongji University (Natural Science), 2020,48(1):24-32.

[28] Parker J C, Lenhard R J, Kuppusamy T. A parametric model for constitutive properties governing multiphase flow in porous media [J]. Water Resource Research, 1987,23(4):618-624.

[29] Peng X, Yu B. Developing a new form of permeability and Kozeny–Carman constant for homogeneous porous media by means of fractal geometry [J]. Advances in Water Resources, 2008,31(1):74-81.

[30] 白青波,李 旭,田亞護(hù),等.凍土水熱耦合方程及數(shù)值模擬研究 [J]. 巖土工程學(xué)報(bào), 2015,37(S2):131-136.

Bai Q B, Li X, Tian Y H, et al. Equations and numerical simulation for coupled water and heat transfer in frozen soil [J]. Chinese Journal of Geotechnical Engineering, 2015,37(S2):131-136.

[31] 劉 偉,范愛武,黃曉明.多孔介質(zhì)傳熱傳質(zhì)理論與應(yīng)用 [M]. 北京:科學(xué)出版社, 2006:28-30.

Liu W, Fang A W, Huang X M. Theory and Application of Heat and Mass Transfer in Porous Media [M]. Beijing: Science Press, 2006:28- 30.

[32] Millington R J. Gas Diffusion in Porous Media [J]. Science, 1959,130(3367):100-102.

[33] Kim J, Corapcioglu M Y. Modeling dissolution and volatilization of LNAPL sources migrating on the groundwater table [J]. Journal of Contaminant Hydrology, 2003,65(1):137–158.

Multi-field coupling simulation of LNAPL migration in cold regions.

FUXiao-qin,SHI Xiao-qing*, JIANGJian-guo, WU Ji-chun

(Key Laboratory of Surficial Geochemistry of Ministry of Education, School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023, China)., 2022,42(3)：1346～1358

A multi-field coupling model of temperature, flow and chemical fields including water, gas, NAPL, and solid (ice) phases was established to quantitatively simulate the migration and distribution of LNAPL under freezing condition during LNAPL leakage. By comparing the LNAPL flow model under non-freezing conditions, the simulated results found that LNAPL would undergo different migration behavior due to the formed ice under freezing conditions. The formed ice prevented LNAPL from infiltrating and generated freezing potential to make LNAPL migrate to the surface. The difference of LNAPL mass flux was -51.1% when the LNAPL infiltration reached stable under freezing/non-freezing conditions. Under freezing conditions, the speed of LNAPL plume in the horizontal direction was greater than that under non-freezing conditions. The maximum predicted error of LNAPL dissolved and volatilized components during the infiltration process was 2.4% and -24.9%, respectively, if freezing condition was not considered. However, the freezing condition had little effect on the dissolved and volatilized components of LNAPL when the LNAPL infiltration reached stable, and only the volatilized component had an error of 1.8%. But the adsorbed component was reduced, and the error was -6.6%.

light non-aqueous fluid (LNAPL)；multi-field multi-phase coupling；freezing potential；solute transport

X131

A

1000-6923(2022)03-1346-13

扶曉琴(1999-),女,湖南婁底人,南京大學(xué)碩士研究生,主要從事地下水動(dòng)力及溶質(zhì)運(yùn)移方面研究.

2021-08-06

國(guó)家自然科學(xué)基金資助重點(diǎn)項(xiàng)目(41730856);國(guó)家自然科學(xué)基金資助項(xiàng)目(41977157)

*責(zé)任作者, 教授, shixq@nju.edu.cn