高彥斌，張松波，李 韜，沈 超

（1.同濟大學(xué)土木工程學(xué)院，上海200092；2.中建三局集團有限公司工程總承包公司，湖北武漢430064；3.湖北中建三局建筑工程技術(shù)有限責任公司，湖北武漢430070；4.上?？辈煸O(shè)計研究院（集團）有限公司，上海200093）

非水相有機污染物（NAPL）是地下環(huán)境污染防治焦點之一，其按照密度分為兩種：比水重的稱為重質(zhì)非水相有機污染物（DNAPL），比水輕的稱為輕非水相有機污染物（LNAPL）。從20世紀60年代開始，國內(nèi)外對NAPL類水土污染問題開展了大量研究工作，尤其是DNAPL在重力作用下的豎向遷移。目前大部分研究工作集中在DNAPL在砂質(zhì)含水層中的與重力作用有關(guān)的縱向遷移、與地下水徑流有關(guān)的橫向遷移以及黏土質(zhì)隔水層頂部自由相DNAPL池的形成［1-4］，在遷移特征和遷移機理方面積累了豐富的研究成果，普遍認為DNAPL-水兩相流以及對流彌散作用是控制DNAPL在砂質(zhì)含水層中遷移的主要因素。

那么，到達滲透率較低的飽和黏性土隔水頂板處的DNAPL是否會再繼續(xù)向下侵入到較深的位置，是一個需要探討的問題。Essaid等［5］在對NAPL的遷移和歸趨的研究綜述中指出，在20世紀80年代以前，人們普遍認為所有黏土或淤泥沉積物都可以有效防止地下水被污染，但隨著對一些地區(qū)污染場地的調(diào)查，發(fā)現(xiàn)承壓含水層也會有這類污染物的存在。我國沿海地區(qū)地表十幾米深度范圍內(nèi)廣泛分布著飽和黏性土（包括粉質(zhì)黏土和黏土兩大類），舊農(nóng)藥廠和化工廠排放了大量的DNAPL污染物。上海地區(qū)舊農(nóng)藥廠區(qū)多次發(fā)現(xiàn)DNAPL在這類場地中遷移到十余米深度的案例。余梅［6］也曾報道了武漢市某化工廠的氯苯類有機物在飽和黏性土地層中遷移至較深位置。國外有學(xué)者認為黏土層中的裂隙（fractures）或孔洞是DNAPL進入土層的優(yōu)先選擇通道［7］；也有學(xué)者認為由于飽和黏性土中的地下水流速較低，分子擴散作用是DNAPL在飽和黏性土中遷移的主要方式［8］?？傮w上講，目前國內(nèi)外對飽和黏性土場地的DNAPL遷移規(guī)律和機理的認識還非常有限。

由于滲透性低，飽和黏性土中DNAPL的遷移會表現(xiàn)出與砂土含水層中不同的特征：①地下水流速?。◣缀蹯o止），因而與流速有關(guān)的對流彌散作用會大大消弱；②土中孔隙較小，DNAPL難以大量侵入，造成了室內(nèi)試驗研究的一些困難，如難以獲得DNAPL的準確含量；③土的滲透性小，遷移速度慢，室內(nèi)試驗周期長。這些客觀因素制約了相關(guān)試驗研究的進展，也是導(dǎo)致目前國內(nèi)外對飽和黏性土場地的DNAPL遷移規(guī)律和機理的認識非常有限的一個重要原因。本文采用數(shù)值分析法進行這方面的研究，討論DNAPL-水兩相流以及對流彌散作用這兩種物理作用對飽和黏性土場地中DNAPL縱向遷移的影響，揭示其主要規(guī)律及機理。

1 飽和土中DNAPL遷移控制方程

1.1 DNAPL的存在狀態(tài)及轉(zhuǎn)化過程

DNAPL在飽和黏土中的存在狀態(tài)有3種：水中溶解相、吸附于土顆粒的固態(tài)以及可流動的自由相。圖1給出了這3種狀態(tài)（相）之間的轉(zhuǎn)換關(guān)系。自由相的DNAPL可以溶解于水中，溶解相DNAPL可以被黏土顆粒部分吸附。溶解相濃度的最大值為其在水中的溶解度。

圖1 DNAPL在土中存在的三種狀態(tài)Fig.1 Three phases of DNAPL in soil

自由相DNAPL在水中的溶解過程可用下式模擬：

式中：De為DNAPL在土中的分子擴散系數(shù)，m2·s-1；d50為土顆粒平均粒徑，m。大量的砂土試驗的結(jié)果證明，參數(shù)Sh的大小與水的流速有關(guān)［9-11］。

溶解相DNAPL吸附于土顆粒的過程用以下線性平衡吸附過程模擬：

式中：Esn/w為吸附速率，g·cm-3·s-1；ρb為土的干密度，g·cm-3；ρ為溶解相DNAPL的質(zhì)量濃度，kg·m-3；Kd為吸附分配系數(shù)，cm3·g-1，吸附分配系數(shù)Kd的大小主要與土中有機質(zhì)含量有關(guān)。

1.2 DNAPL遷移控制方程

DNAPL污染源在飽和黏性土中的縱向遷移主要包括兩種過程：一種是自由相DNAPL以連續(xù)流體方式的遷移，這個過程受DNAPL-水毛管壓力的控制，稱為兩相流；另外一種是DNAPL溶解于水中后在水中的遷移，也被稱為溶質(zhì)運移，這個過程受對流彌散作用的控制，但在靜水條件下，主要受分子擴散作用的控制。根據(jù)土孔隙中各物質(zhì)的質(zhì)量平衡，水溶液、自由相DNAPL和溶解相DNAPL的控制方程分別為

式中：ε為土的孔隙度；Sw和Snw分別為水和自由相DNAPL的飽和度；vw和vnw分別為水和自由相DNAPL的達西流速，m·s-1；Dw為水動力彌散系數(shù)，m2·s-1；ρw和ρnw分別為水溶液和自由相DNAPL的密度，g·cm-3和分別為溶解速率和吸附速率。方程（4）和（5）構(gòu)成了經(jīng)典的兩相流方程，方程（6）是經(jīng)典的對流彌散溶質(zhì)運移方程。

水動力彌散系數(shù)Dw是個二階張量，包括機械彌散（與水的流速有關(guān)）和分子擴散（與水的流速無關(guān)）兩部分，表達式為

式中：uwx、uwy和uwz為水的真實流速分量，m·s－1；|uw|為水的真實流速大小，m·s－1；αL和αT分別為縱向和橫向彌散度，m，對于黏性土，典型數(shù)值為αL＝0.1～1.0 m，αT＝（0.1～0.3）αL；De為有效分子擴散系數(shù)，m2·s－1，與DNAPL的分子擴散系數(shù)以及土的結(jié)構(gòu)性有關(guān)，即

式中：Daq是DNAPL在水中的分子擴散系數(shù)，m2·s－1），一般在10-10～10-9m2·s－1；ω為彎曲因子，與土的結(jié)構(gòu)性有關(guān)，一般在0.1～0.5之間。

根據(jù)達西定律，真實流速uw（水）、unw（NAPL）以及達西流速vw（水）、vnw（NAPL）之間的關(guān)系可表示為

式中：κ為土體固有滲透率，m2；krw和krnw分別為水和DNAPL的相對滲透率；pw和pnw分別為水和DNAPL的壓力，Pa。

毛管壓力pc（=pnw－pw）、相對滲透率kr和飽和度S三者之間的關(guān)系，即kr-S-p關(guān)系，是兩相流的重要本構(gòu)關(guān)系。根據(jù)Parker等［12］提出的縮放理論，由van Genuchten［13］模型給出的DNAPL-水毛管壓力pc與水的飽和度Sw的關(guān)系如下：

式中：β 為水-氣表面張力 σaw（=72 mN·m-1）與DNAPL-水表面張力σow（mN·m-1）的比值，對于大多數(shù)DANPL，σow=30～45 mN·m-1；m和n為土水特征曲線的參數(shù)，且m=1－1/n。

采用 Mualem［14］提出的方法，由 van Genuchten模型得到的相對滲透率krw（水）、krnw（DNAPL）與水的飽和度Sw的關(guān)系

注意相對滲透系數(shù)只與參數(shù)m（m=1－1/n）有關(guān)，而和α無關(guān)。

1.3 方程求解

基于COMSOL Multiphysics軟件開發(fā)了求解偏微分方程（4）～（6）的模塊，需求解的變量為pw、pnw和ρ。COMSOL Multiphysics軟件是一種自動有限元軟件，采用有限元方法求解偏微分方程組。計算結(jié)果輸出中土中自由相DNAPL的含量用飽和度表示，溶解相DNAPL的含量用濃度表示，二者均可以轉(zhuǎn)化為質(zhì)量分數(shù)表示，即單位質(zhì)量的土（飽和土）中DNAPL的質(zhì)量。

2 分析模型及基本參數(shù)

DANPL在飽和黏性土中一維遷移數(shù)值分析模型如圖2所示。模擬對象為深度為10 m的均質(zhì)飽和黏性土土柱。DNAPL選取氯苯，根據(jù)國家相關(guān)標準［15］規(guī)定，這種污染物在土中的篩選值為68 mg·kg－1（即0.068 mg·g-1）。土柱頂面的水壓pw為0，自由相DANPL壓力pnw為5 kPa（即0.5 m水頭），溶解相DANPL的質(zhì)量濃度ρ0為溶解度S*。土柱底面為一排水界面，水壓pw維持在100 kPa。黏性土中的孔隙水并沒有賦予一個流速，地下水處于靜止狀態(tài)。即使在自由相DNAPL的驅(qū)替下，黏性土中水的流速也非常低，這一點從后面給出的分析結(jié)果中可看出。

圖2 數(shù)值分析模型Fig.2 Model of numerical analysis

計算參數(shù)分為三類：第一類為兩相流方程參數(shù)，第二類為對流彌散方程參數(shù)，第三類為溶解作用參數(shù)。這些參數(shù)的取值如表1所示，所取參數(shù)值均為典型數(shù)值。土的參數(shù)選用的是沿海地區(qū)典型飽和粉質(zhì)黏土的參數(shù)，孔隙度ε取0.5，固有滲透率κ為1×10－15m2（飽和滲透系數(shù)為1×10－6cm·s－1）。土的Sw-Pc模型參數(shù)取值參照Lu等［16］給出的黏土的參數(shù)的變化范圍，α在0.001～0.01 kPa－1之間，n在1.1～2.5之間。對于溶解作用參數(shù)Sh=2，參照了一些學(xué)者在水流緩慢情況下得到的NAPL在砂土中的溶解速率值［9-11］。

表1 參數(shù)基本值Tab.1 Basic values of parameters

3 靜水條件下溶解相DNAPL分子擴散

首先僅模擬溶解相DNAPL自土柱頂面向下遷移的過程，不考慮自由相DNAPL的兩相流遷移，物理上為靜水條件下的分子彌散過程。土柱頂面的邊界條件為：①溶解相DANPL的質(zhì)量濃度ρ0為溶解度S*；②水壓力pw為0；③自由相DANPL的壓力為0。

分析得到的質(zhì)量分數(shù)wad（包括溶解相和吸附相，為其質(zhì)量與土的總質(zhì)量之比）以及溶解相的質(zhì)量濃度比ρ/ρ0隨深度h的變化如圖3所示。在選取的參 數(shù) 值 下（De=2.9×10－10m2·s－1，Kd=3.63cm3·g－1），遷移過程是非常緩慢的，10年的最大遷移深度為0.5 m，平均每年0.05 m。溶解相DNAPL含量的最大值出現(xiàn)在地表，質(zhì)量分數(shù)為1.4 mg·g－1，隨深度的增大而減小。

有效分子擴散系數(shù)De和分配系數(shù)Kd是影響靜水條件下的分子擴散的最重要的兩個參數(shù)，下面分析這兩個參數(shù)的取值對計算結(jié)果的影響。

3.1 有效分子擴散系數(shù)De

有效分子擴散系數(shù)De的大小決定了溶解相DNAPL分子在土中擴散的速度。統(tǒng)計分析結(jié)果表明，黏土中NAPL有效分子彌散系數(shù)De在（1～5）×10－10m2·s－1之間［17-19］。將De分別取 1×10－10m2·s－1和1×10－9m2·s－1，其他參數(shù)按照表1 中的對流彌散方程參數(shù)設(shè)置，計算得到的溶解相DNAPL在土中的遷移過程如圖4所示。有效分子擴散系數(shù)增大10倍，DNAPL遷移速度僅增大到原來的3倍左右，10年的最大遷移深度也僅為1 m左右。

3.2 分配系數(shù)Kd

圖3 靜水條件下溶解相DNAPL的分子彌散Fig.3 Molecular diffusion of dissolved DNAPL under hydrostatic condition

圖4 有效分子擴散系數(shù)參數(shù)分析Fig.4 Parameter analysis of effective molecular diffusion coefficient

黏土中有機質(zhì)含量較高，故研究吸附作用對溶解相DNAPL遷移的影響很有必要。線性吸附分配系數(shù)Kd分別取0 cm3·g－1（不考慮吸附）和10 cm3·g－1，其他溶質(zhì)運移參數(shù)為表1中給出的數(shù)值，得到的計算結(jié)果如圖5所示。從圖5可知，分配系數(shù)Kd越大，阻滯作用越明顯，DNAPL遷移速度越慢。在不考慮吸附的情況下（Kd=0 cm3·g－1），10年的遷移深度增大到2 m左右。可見，飽和黏性土的吸附作用會顯著影響溶解相DNAPL的遷移。

以上分析結(jié)果表明，水流靜止條件下，加上吸附所導(dǎo)致的阻滯，DNAPL在飽和黏性土中的彌散遷移速度非常緩慢，十幾年和幾十年內(nèi)不會造成飽和黏性土的深部污染。

4 連續(xù)注入下兩相流-對流彌散遷移分析

圖5 線性吸附分配系數(shù)參數(shù)分析Fig.5 Parameter analysis of distribution coefficient of linear sorption

下面進行連續(xù)注入情況下的兩相流與對流彌散耦合分析。采用圖2所示的邊界條件，在土柱頂面的自由相DNAPL施加0.5 m的水頭，即5 kPa的壓力（也就是約45 cm氯苯的高度）。這種情況下自由相的DNAPL會在重力作用下逐漸向下遷移。下面主要介紹分析得到的自由相DNAPL的縱向遷移規(guī)律，以及參數(shù)取值對分析結(jié)果的影響。

4.1 自由相DNAPL含量、壓力以及流速

計算得到的不同時間的自由相DNAPL的質(zhì)量分數(shù)wnw（自由相DNAPL的質(zhì)量與土的總質(zhì)量之比）和飽和度Snw隨深度h的變化如圖6所示?？梢钥闯?，10年后的遷移深度達到2.5 m左右，大于圖3給出的不考慮兩相流情況下的遷移深度，平均遷移速度為0.25 m·年－1。自由相DNAPL的飽和度在0.003以內(nèi)，質(zhì)量分數(shù)在1 mg·g－1以內(nèi)，與前面給出的溶解態(tài)DNAPL的質(zhì)量分數(shù)處于同一個數(shù)量級。盡管這個數(shù)量級是非常低的，但是仍然超出污染物氯苯的篩選值（0.068 mg·g－1）。因此，兩相流遷移可能是導(dǎo)致深部飽和黏性土中DNAPL檢出值超標的原因。

圖7給出了DNAPL壓力pnw、水壓力pw以及毛管壓力pc隨深度的變化。可以看出，水壓力pw幾乎始終維持在靜水壓力，而DNAPL壓力pnw則隨著遷移深度的增加逐漸向深部傳遞。毛管壓力pc決定了自由相DNAPL的飽和度。毛管壓力最大值為5 kPa，出現(xiàn)在地表。

圖6 飽和黏性土中自由相DNAPL的含量Fig.6 Content of free DNAPL in soil

圖7 pnw、pw和pc隨深度的變化Fig.7 DNAPL pressure,water pressure，and capillary pressure

圖8 給出了vnw和vw隨深度的變化。在自由相DNAPL遷移深度范圍內(nèi)，水的達西流速隨深度的增加而增大，DNAPL的達西流速隨深度的增加而減小。下面來分析造成這種現(xiàn)象的原因。隨著深度的增大，DNAPL的飽和度減小導(dǎo)致相對滲透率減小，而壓力梯度略有增大，因而達西流速減小。對于孔隙水，隨著深度的增大，水的飽和度增大導(dǎo)致滲透系數(shù)增大，而壓力梯度基本不變，因此水的達西流速增大。在DNAPL侵入深度以下，水的飽和度為100%而水力梯度不變，因此水的速度沿深度基本不變。

圖8中給出的DNAP和水的達西流速為一個數(shù)量級，在0～0.3 mm·年－1之間，這個速度是非常小的。然而由于二者飽和度的巨大差別（自由相DNAPL的飽和度在0.003以內(nèi)，而水的飽和度接近1），因此根據(jù)式（9）得到的DNAPL和水的真實流速也差別較大。圖9給出了DNAPL和水的真實流速uw和unw，DNAP的真實流速達到100～200 mm·年－1，是水的幾百倍，這與圖6中給出的平均每年遷移0.25 m的速度是吻合的。注意DNAPL真實流速隨深度的增大而增大，與達西流速隨深度的變化正好相反。在DNAPL飽和度接近0的深度處，DNAPL真實速度會出現(xiàn)突變的現(xiàn)象，這是因為計算真實速度時，當式（9）分母中的飽和度接近0時得到的結(jié)果會不穩(wěn)定。

圖8 vw和vnw隨深度的變化Fig.8 Darcy velocities of water and DNAPL

圖9 uw和unw隨深度的變化Fig.9 Real velocities of water and DNAPL

4.2 自由相DNAPL的溶解

為了對自由相NDAPL在土中的溶解速度有一個直觀的印象，這里先分析靜止條件下（DNAPL和水的流速均為0）自由相DNAPL在土中的溶解。靜止條件下自由相DNAPL的溶解可以表示為

上述微分方程的解為

取氯苯的 De=3×10－10m2·s－1，黏土 d50≈0.01 mm，孔隙度ε=0.5，水的飽和度Sw=0.9，修伍德準數(shù)Sh=2，得到

令ρ（t）/S=0.999，計算得到t=0.52 s。意味著在t=0.52 s時，水中氯苯的質(zhì)量濃度就達到了溶解度。溶解平衡過程是迅速的，其原因是黏土的平均粒徑d50要遠遠小于砂土，因此按式（2）計算得到的溶解平衡速率 kf較大。Lenczewski等［20］通過實驗發(fā)現(xiàn)，TCE（三氯乙烯）能夠迅速進入富含黏土的細粒土中，并快速溶解在孔隙水中。因此可以認為自由相DNAPL在飽和黏性土中的溶解平衡過程幾乎是瞬時的。

圖13給出了分析得到的溶解相DNAPL在土中的分布。與圖3對比可以看出，自由相DNAPL的溶解使得溶解相DNAPL的分布發(fā)生了顯著的變化。在自由相DNAPL遷移深度范圍內(nèi)，溶解相DNAPL的質(zhì)量濃度均達到溶解度。且由于機械彌散作用，10年時溶解相的最大遷移深度比自由相的（約2.5 m）還要深，達到3.0 m左右。以上分析表明，飽和黏性土中溶解相DNAPL含量受兩相流和溶解作用控制，而不是對流彌散作用。兩相流分析對土中溶解態(tài)DNAPL含量的分析同樣非常重要。

圖10 溶解相DNAPL的分布（Sh=2）Fig.10 Profile of dissolved DNAPL(Sh=2)

4.3 滲透率的影響

土體固有滲透率κ是影響自由相DNAPL遷移的一個最重要的參數(shù)，它決定了自由相DNAP的遷移速率。 將κ從1×10－15m2減小至1×10－16m2和1×10－17m2，也就是飽和滲透系數(shù) ks從 1×10－6cm·s－1減小至 1×10－7cm·s－1和 1×10－8cm·s－1，其余參數(shù)均不變。這種滲透性的變化相當于由滲透性較大的粉質(zhì)黏土變?yōu)闈B透性低的壓實黏土。計算得到的自由相DNAPL的遷移剖面如圖11所示。

圖11 固有滲透率不同取值的分析結(jié)果Fig.11 Paramter analysis of intrinsic pemibility

從圖11可以看出，滲透系數(shù)降低以后，自由相DNAPL的遷移速度和遷移深度均降低。當κ減小100倍時，遷移深度減小至原來的1/10。只有當飽和黏性土的滲透率很低時，自由相DNAPL向下的遷移才會非常緩慢。對于天然飽和粉質(zhì)黏土，其滲透系數(shù)往往會達到1×10－6cm·s－1的數(shù)量級，這種地層中自由相DNAPL的向下遷移是不可忽視的。

4.4 Sw-pc關(guān)系的影響

Sw-pc關(guān)系對兩相流有直接的影響。改變參數(shù)α和n的數(shù)值，參數(shù)α分別取0.002 5（減?。┖?.010 0（增大），參數(shù)n分別取1.2（減?。┖?.0（增大）以了解Sw-pc關(guān)系對自由相DNAPL遷移的影響。計算得到的自由相DNAPL的遷移剖面分別如圖12和圖13所示。

從計算結(jié)果來看，參數(shù)α對DNAPL飽和度的影響較大，對DNAPL侵入深度的影響較小。參數(shù)n對DNAPL遷移速度和遷移深度影響很大，在n=1.2的情況下，連續(xù)注入10年的遷移深度可達到10 m；而在n=2.0的情況下，連續(xù)注入10年的遷移深度只有1 m。因此，遷移速度和深度對參數(shù)n的取值非常敏感。這和公式（11）是一致的，相對滲透率僅與參數(shù)n有關(guān)，而和α無關(guān)。

圖12 參數(shù)α不同取值的計算結(jié)果Fig.12 Parameter analysis of α

圖13 參數(shù)n不同取值的計算結(jié)果Fig.13 Parameter analysis of n

4.5 DNAPL污染源消失后的縱向遷移

前面分析了土柱表面DNAPL污染源在5 kPa壓力下連續(xù)注入10年情況下的DNAPL遷移特性。在連續(xù)污染10年后，將土柱頂部的DNAPL壓力設(shè)為0，代表工廠關(guān)閉等原因?qū)е碌奈廴驹聪?，分析得到了污染源消失?0年DNAPL在土中的遷移過程。

圖14給出了污染源消失后（階段B）自由相DNAPL的飽和度剖面，并與前期連續(xù)注入10年（階段A）的計算結(jié)果對比。頂部污染源消失后，DNAPL完全在重力作用下繼續(xù)向下遷移，最大飽和度出現(xiàn)的深度也逐漸下移。由于DNAPL總量基本不變，因此隨著遷移深度的增加，最大飽和度逐漸減小。在第60年時最大遷移深度達到近6 m，平均遷移速度約0.08 m·年－1，比階段A的平均遷移速度0.25 m·年－1要小一些。

圖14 污染源消失后自由相DNAPL的分布Fig.14 DNAPL distribution after removing pollution source

5 結(jié)論

基于COMSOL Multiphysics軟件平臺，開發(fā)了DNAPL-水兩相流分析以及考慮吸附、溶解作用的溶質(zhì)運移分析模塊，進行了均質(zhì)飽和黏土中DNAPL縱向遷移過程的分析，得到了一些有意義的結(jié)論。

（1）飽和黏性土中溶解態(tài)DNAPL的分子擴散遷移過程非常緩慢，不能在幾十年內(nèi)造成土體深部污染。自由相DNAPL的縱向遷移速度主要受固有滲透率k以及毛管壓力-飽和度關(guān)系的影響，可在一定條件下成為飽和黏性土中DNAPL豎向遷移的主控因素。

（2）毛管壓力-飽和度模型中參數(shù)n以及土的固有滲透率k控制著自由相DNAPL的縱向遷移速度以及含量。盡管在低毛管壓力下自由相DNAPL難以大量侵入飽和黏性土（飽和度只有千分之幾），但仍然會超過污染物的環(huán)境篩選值。對于滲透系數(shù)較大的飽和黏性土地層（如粉質(zhì)黏土層），這種微量的、真實流速可達每年幾十厘米的流動足以造成DNAPL在幾十年內(nèi)遷移到十幾米的深度。

以上基于數(shù)值分析得到的結(jié)論，能夠初步解釋一些舊農(nóng)藥廠廠區(qū)飽和黏性土場地中深部存在DNAPL污染的現(xiàn)象。關(guān)于飽和黏性土中的DNAPL-水兩相流分析，仍然需要進行更為深入和系統(tǒng)的研究工作。