基于TMVOC的有機污染場地熱強化氣相抽提修復效果模擬

孫 超，徐 超，孫宗全，沈佳倫，龍 杰，李曉東*，馬福俊，谷慶寶

1. 浙江工業大學環境學院，浙江 杭州 310014

2. 中國環境科學研究院，環境基準與風險評估國家重點實驗室，北京 100012

進入21世紀以來，城市化進程不斷加快，城區建設與污染嚴重的老舊企業之間的矛盾愈發不可調和，伴隨這些工業企業拆除搬遷而來的是大量的遺留污染場地. 在石油化工、金屬冶煉、機械制造、燃氣生產等多種行業企業中，有機化合物作為原料、中間產物或產品在其生產、運輸和儲藏過程中通過“跑冒滴漏”等方式進入土壤造成土壤污染[1-2]. 有機污染土壤具有種類繁多、危害性大等特點，嚴重威脅人類身體健康和生態環境[3-4]. 因此，如何精準有效地對遺留有機污染場地進行修復迫在眉睫[5-6].

土壤氣相抽提技術(Soil Vapor Extraction，SVE)作為一種常見的土壤原位修復技術，具有成本低廉、操作簡單、效率高、對周邊環境影響小等優勢，能夠高效去除土壤中的揮發性有機物[7-9]. 然而在常溫情況下，SVE技術易受到土壤理化性質、污染物特性等因素的局限，通過與原位熱修復技術進行耦合，可以擴大其適用范圍并提高修復效率[10-15]. 其中，原位熱傳導(Thermal Conduction Heating, TCH)技術可通過土壤顆粒及孔隙流體的導熱作用使土壤升溫以實現有機污染物的高效去除[16]，因其修復周期短、二次污染可控、對土壤質地以及污染物性質的適用性強等優點，被普遍認為是一種極具前景的污染場地修復技術[17-19].

目前，污染場地修復工藝的設計與運行大多參照經驗公式與有限的場地實踐經驗. 然而，工藝修復效果受土壤質地[20]、水文地質特性[21]、井位布設[22]、通氣速率[23]等多種因素共同影響，且熱處理過程的熱能損失難以避免，對TCH技術而言，在擬修復區域頂部和底部的熱量逸散尤其明顯[11,24]. 因此，在國家“雙碳”目標和“可持續發展”的政策背景下，減少能源浪費，對污染場地進行系統化、精準化修復，必將是未來污染場地修復行業的趨勢所在. 數值模擬技術作為一種科學研究與工程規劃的重要手段，能夠有效揭示研究區溫度場的時空分布規律[25-27]，深入闡明揮發性有機物(VOCs)在土壤中的遷移轉化規律和SVE治理機制[28-30]. 例如，Wang等[31]通過COMSOL軟件對土壤加熱過程中的傳熱性能進行了研究，評估了不同工況下的溫度分布和加熱周期，模擬結果可為現場試驗提供參考；Huang等[32]建立了一種土壤氣相抽提井運行過程中污染物濃度空間分布及時間演變的模型，可用于設計SVE的修復策略. TOUGH 2是一款能夠精準有效地模擬多相流和熱量運移過程的數值模擬軟件[33-34]，TMVOC作為其子模塊，側重應用于烴類或揮發性有機物在土壤及地下水中的環境污染問題，不僅可以模擬“自然”條件下多維非均相系統中三相(氣、水、NAPL)的非等溫流動，還可以模擬多孔介質的熱量傳遞及污染物的運移、去除過程[34-37].

目前相關領域的模擬研究多集中于小試規模，未充分考慮實際污染場地的復雜性. 筆者以天津某試劑廠中試場地為研究對象，采用TMVOC模型模擬TCH強化氣相抽提過程中場地的溫度變化及目標污染物氯苯的去除效果，結合中試數據進一步闡明場地升溫過程及氯苯的去除規律，并驗證TMVOC模型的可信度，為工程化應用打下堅實的基礎.

1 研究區概況

研究區為天津某化學試劑廠的中試場地，原廠主要生產有機通用試劑、指示劑和基準試劑等，現已關停搬遷. 調查資料顯示，研究區存在較嚴重的有機污染，污染面積278.9 m2，地下水埋深15.0 m，土壤污染深度為10.0 m，主要污染物氯苯的最高濃度達30.00 mg/kg[20]. 土壤采樣點分布見圖1，不同深度處氯苯初始濃度見表1. 地質勘查結果表明，土壤污染區域自上而下分為人工填土、粉質黏土、砂質粉土夾粉質黏土、粉質黏土夾粉土薄層團塊5個工程地質層，各類土壤理化性質見表2.

表1 地下不同深度各采樣點氯苯的濃度Table 1 Concentrations of chlorobenzene at different depths and sampling points in soil

表2 研究區土壤理化性質Table 2 Physical and chemical properties of soil in study area

圖1 研究區采樣點分布Fig.1 Distribution of sampling points in study area

研究區采用“TCH+SVE”組合修復技術，設置3 m(TCH-A區)和4 m (TCH-B區)兩個不同間距的加熱井群進行分區修復. 系統主要由電力控制單元、加熱單元、溫度監測單元、抽提單元和廢水廢氣處理單元5個部分組成. 加熱井采用正三角形布點法，抽提井分布于三角形質心位置. 工藝設備及井位布設如圖2所示.

圖2 TCH原位熱修復系統和井位布設示意Fig.2 In situ thermal remediation system of TCH and layout of wells location

研究區布設加熱井45個，測溫井19個，抽提井19個. 在TCH-A、TCH-B兩個分區內分別針對一加熱井，在其水平方向不同距離處設置一系列測溫井以分析單井有效加熱范圍內的溫度分布規律(見圖3). TCH-A區設T1～T4四個測溫井，分別距熱源中心0.5、1.0、1.5、1.73 m (冷 點). TCH-B區 設T8～T12五個測溫井，分別距熱源中心0.5、1.0、1.5、2.0、2.3 m (冷點). 為監測兩個分區區域內平均溫度，在TCH-A區另設3個冷點測溫井，編號為T5～T7，TCH-B區另設兩個測溫井，編號為T13～T14. 系統運行7周后，TCH-A區平均溫度為100～110 ℃，TCH-B區平均溫度為80～90 ℃. 維持該溫度一周后停止加熱并采樣檢測，土壤中氯苯僅有少量殘余，且濃度低于修復目標值.

圖3 研究區測溫井布設示意Fig.3 Layout of temperature measuring wells in study area

2 數值模型概化

2.1 計算方程

2.1.1多相流控制方程

TMVOC模型假設由水、揮發性有機物和不可壓縮氣體組成多相流系統，基于組成部分相對豐度的不均衡及不同的熱力學條件，系統在進行相間傳質時流體組分可能會出現7種不同的相態組合，如圖4所示.

圖4 TMVOC模型中流體的相態變化Fig.4 The phase change of fluid in TMVOC model

各組分之間的箭頭代表因熱力學條件改變而導致Newton-Raphson迭代過程中的相間傳質路線[38-39].假設同種組分不同相態下的化學勢相同，則有：

模擬過程中設定理想混合組分中的活度系數和逸度系數為1，根據摩爾分數的比例采用平衡常數表示各相態間的分配系數：

TMVOC模型中質量和能量平衡的總控制方程：

式中：Vn為流動單元體體積，m3；Mk為組分k在單位土壤介質中的質量，kg；Γn為表面積，m2；Fk為進入到流動單元體組分k的總通量；n為流動區單元體表面的外法向單位矢量；qk為組分k在單元體的源匯項.

2.1.2熱傳導控制方程

土壤一般被認為是一種含濕多孔介質，在TCH過程中的傳熱可以通過熱量的傳導、對流、相變、輻射以及水分的擴散，壓力驅動等方式進行. TMVOC模型假設條件：土壤外部絕熱；土壤內部質地均勻；土壤熱性能穩定；忽略污染物的傳熱；土壤中水的密度、導熱性及比熱容恒定不變.

根據熱力學第一定律，熱傳導控制方程利用局部容積平均法可表示為

式中：ρeff為有效密度，kg/m3；ceff為有效比熱容，J/(kg·K)；T為溫度，℃；t為時間，s；λ為導熱系數，W/(m·K)；?為蒸發速率，kg/(m3·s)；ΔHvap為水蒸氣的熱量變化，J；θS、θL、θG分別為土壤孔隙率、液態水體積分數和水蒸氣體積分數；ρS、ρL、ρG分別為土壤、液態水和水蒸氣的密度，kg/m3；CS、CL、CG分別為土壤、液態水和水蒸氣的比熱容，J/(kg·K).

濕度控制方程為

式中，DL為表觀液體擴散系數.

依據ANTOINE方程以解釋土壤中氣體的飽和蒸汽壓、平衡蒸汽壓與蒸發速率的關系：

式中：α為比例常數；θL*表示殘余飽和度；kvap為蒸發速率常數，s-1；P*為飽和蒸汽壓，Pa；PG為平衡蒸汽壓，Pa；A、B、C均取常數[40].

2.2 模型構建與網格剖分

將研究區概化為面積278.9 m2、高度11.2 m的三維模型，不考慮0～2 m的人工清挖土壤及地下水.加熱棒為半徑0.2 m的圓柱，加熱深度為地下2～10 m，抽提井采用TMVOC模型內置模塊，抽提深度為地下2.0～10.5 m. 參數設置、工藝布設方式與場地設置一致，研究區概念模型如圖5所示.

圖5 幾何概念模型示意Fig.5 Schematic diagram of geometric conceptual model

模型依據笛卡爾X、Y、Z正交坐標系，X、Y方向上網格剖分采用不規則多邊形填充，網格大小決定模型的細化程度，網格越小，模型的細化程度越高，可分析的網格數據越多，模擬結果越精準. 當網格最大面積設置低于0.5 m2時，網格數量超出模型求解最大范圍導致模型不收斂，故設置網格最大面積為0.5 m2.最小改良角度越大，加熱井周邊網格的細化程度越高，根據模型的實際需求確定最小改良角度為15°.Z方向上共剖分為10層，從上至下依次為1層0.001 m的大氣邊界，1層0.2 m的保溫層，5層1.5 m的污染層，1層0.5 m和1層4.0 m的無污染層，以及1層1.0 m的底部固定邊界. 根據場地土壤初始溫度實測數據，模型初始溫度設置為15 ℃. 研究區四周設有止水帷幕，故將模型四周設置為無通量邊界.

2.3 條件初始化與初始濃度賦值

模型構建完成后，設置土壤參數(見表2)、污染物參數、模型參數和網格中的初始條件，運行程序1年使模型自動計算至重力和毛細壓力平衡以模擬場地自然初始狀態.

研究選取氯苯為污染物進行模擬，初始污染濃度采用“點泄漏”的方式進行賦值，即在包氣帶采樣點處及其附近的網格設置源匯項，以一定的速率注入NAPL相氯苯，模擬氯苯的實際泄漏過程，一定時間后去除泄漏點并使之自由遷移擴散1年，最終得到研究區氯苯的模擬初始濃度狀態. 采樣點位置見圖1，由于各采樣點在不同深度的污染濃度不一致且無分布規律(見表1)，故對污染區域各層網格進行單獨賦值. 為防止賦值過程中的垂向遷移擴散，在進行單層網格濃度賦值時禁用其他層. 為降低賦值方式帶來的污染物初始遷移擴散速率過高的影響，賦值時先將氯苯在空氣中的擴散速率增大為3×10-5m2/s，按一定的速率泄漏1年后去除泄漏點，再將擴散速率改為正常值(7.6×10-6m2/s)，自由遷移擴散1年使之達到目標污染濃度.

2.4 模擬過程

以2.3節污染物賦值完成的模擬結果作為模型的初始運行條件，加熱功率設為520 W/m，抽提速率5×10-4m3/s，抽提壓強70 kPa，加熱井與抽提井同時運行8周. 模擬相關參數設置均與場地中試參數一致，不考慮生物降解作用.

3 模擬結果分析與驗證

3.1 溫升模擬效果

模型在TCH-A區和TCH-B區冷點處分別設置了4個和3個測溫井(冷點為加熱井圍成的正三角形的質心)，根據模型輸出的測溫井數據可得到冷點在不同深度處的溫度，冷點垂向上的溫度取平均值作為其垂向平均溫度，區域內各冷點垂向平均溫度取平均值作為區域整體平均溫度. 根據冷點測溫井模擬數據可知，經過7周熱強化氣相抽提系統的持續運行，加熱井間距為3 m的TCH-A區土壤整體平均溫度達到99.4 ℃，加熱井間距為4 m的TCH-B區為84.8 ℃.結果表明，加熱井間距直接影響場地的升溫效果，間距越小，升溫效果越好. 垂直方向上，地表以下2.75、4.25、5.75、7.25、8.75 m各深度的升溫效果不一致.2.75 m處TCH-A和TCH-B兩個分區的平均溫度分別為94.3和74.0 ℃；4.25 m處分別為97.9和89.1 ℃；5.75 m處分別為106.2和89.2 ℃；7.25 m處分別為99.9和91.8 ℃；8.75 m處分別為98.8和79.8 ℃.

模型區域溫度分布以及垂向溫升差異如圖6所示. 由圖6(a)可知，頂、底兩層土壤溫度較中間三層土壤更低，該現象歸因于該兩層土壤在加熱過程中能量向外界發生了擴散. 頂層土壤接近地表，即使設置保溫層仍有較大熱量散失，底層土壤以下區域也有一定的升溫，說明加熱過程中底層土壤熱量向深部地層發生了擴散. 模擬加熱7周后，研究區溫度達到目標溫度. 隨后維持該溫度1周，TCH-A區的平均溫度降至95.0 ℃，TCH-B區升至88℃ ，其主要原因為土壤介質之間產生了熱量傳遞，溫度高的TCH-A區將一部分熱量傳遞給了溫度更低的TCH-B區〔見圖6(b)〕.

圖6 TCH系統運行8周后溫度場分布Fig.6 Temperature field distribution of TCH system after 8 weeks of operation

為進一步研究加熱間距對溫度場升溫過程的影響，獲取溫度監測井的模擬數據進行分析，兩個分區距單井熱源中心不同位置處的土壤溫升曲線見圖7.土壤溫升過程包括加熱、潛熱和過熱3個階段[26].TCH-A區T1、T2點位加熱7周后的溫度均超過100.0 ℃，因其受到的熱量供給相對較高，在加熱階段達到100.0 ℃以后進入潛熱階段，持續加熱使土壤孔隙水蒸發，待蒸發完全后進入過熱階段，溫度繼續上升. 由于T2所受熱量供給低于T1，故其潛熱階段相對較長. T3和T4點位溫度分別為100.7 和100.5 ℃，均處于潛熱階段. TCH-B區5個點位的模擬溫度均未超過100.0 ℃，其中T8、T9點位溫度接近100.0 ℃，即將進入潛熱階段，T10～T12點位均處于加熱階段. 綜上可知，加熱井間距對場地溫升效果影響顯著，加熱井間距越小，升溫速率越快. 因此，合理布設加熱井間距可在達到修復效果的前提下降低能量損耗，節約成本.

圖7 距單井熱源不同位置處溫升曲線Fig.7 Temperature rise curve at different positions from single well heat source

3.2 VOCs 修復模擬效果

模型初始污染狀態、運行7周及維持目標溫度1周后的污染狀態見圖8. 結果顯示，S01～S04采樣點在土壤不同深度處氯苯的檢出濃度最高為0.1 mg/kg，表明土壤中大部分氯苯被去除. 地下2.0～3.5、3.5～5.0、5.0～6.5、6.5～8.0、8.0～9.5 m處氯苯的平均去除率分別為97.9%、97.0%、96.3%、95.4%、95.1%. 由圖8可知，模型頂、底兩層土壤中的氯苯全部去除，中間三層土壤中的氯苯大部分被去除，TCH-B區的S04采樣點附近還有少量氯苯殘留. 出現該現象的主要原因是：①頂、底兩層土壤氯苯的初始濃度較低，去除速率快；②TCH-B區的平均溫度低于TCH-A區，污染物的揮發速率相對較慢. 繼續加熱停留1周后，污染物基本去除完全，達到修復目標. 由此表明，適當延長目標溫度維持時間可提高污染物的修復效率.

圖8 氯苯濃度隨時間的變化規律Fig.8 Variation of chlorobenzene concentration with time

通過收集模型輸出數據可知，模型初始濃度賦值后氯苯的總質量為13.27 kg，經過8周的熱強化氣相抽提修復后僅剩0.36 kg，去除率達97.3%(見圖9). 通過氯苯總質量變化計算氯苯的去除速率，發現其去除速率呈先增后降趨勢，第1周的去除速率最慢，為0.07 kg/d，第4周的去除速率最快，為0.46 kg/d. 出現該現象的原因主要是：①初始階段溫度場升溫較慢，低溫環境下氯苯分子的熱運動受到限制，導致修復效率較低；②隨著加熱時間延長，土壤溫度逐漸升高進而提升修復效率；③在修復后期由于氯苯濃度的降低，土壤顆?？紫吨休^易揮發位置處的氯苯已基本去除，殘留于土壤顆粒結構內部及液相中的氯苯只能通過緩慢擴散進入氣孔通道被揮發去除[41-42].

圖9 熱強化SVE修復期間氯苯總質量的變化趨勢Fig.9 Trend of total chlorobenzene mass during thermal enhanced SVE remediation

3.3 擬合優度評估

選取TCH-A區和TCH-B區加熱過程中各時間點溫度場平均溫度的試驗值與模擬值、修復后氯苯濃度的試驗值與模擬值分別進行擬合優度評估，統計數據如表3、4所示. 選取R作為擬合優度的評估指標[35,43]，該值越接近1代表模型的擬合程度越高.

表3 TMVOC模型溫度擬合數據Table 3 Temperature fitting data of TMVOC model

表4 TMVOC模型污染濃度擬合數據Table 4 Concentration fitting data of TMVOC model

R的計算公式：

式中：cali為兩個分區溫度的試驗值為試驗值的平均值，g/m3；obsi為兩個分區溫度的模擬值，為模擬值的平均值，g/m3；N表示一組擬合數據的第N個值.

根據計算可得，TCH-A區和TCH-B區溫升數據的擬合優度分別為0.995和0.989，修復后氯苯濃度數據的擬合優度為0.914，表明模型計算的模擬值與試驗值相關性良好，模型可信度較高. 模擬結果與試驗結果產生偏差的原因主要有以下幾點：①中試試驗過程中受水文和地質等復雜環境條件影響，而模擬程序運行時僅受計算機中央處理單元計算精度的影響；②土壤的相對滲透系數和與孔隙壓力等相關參數難以獲取，模型采用經驗參數代替；③TMVOC模型不能直接進行污染物濃度賦值，通過“點泄漏”方式進行初始濃度賦值會帶來一定的誤差. 值得注意的是，TMVOC模型并無內置加熱模塊，該模型所設計的加熱元件實際上為實心管，通過管體加熱向周邊土壤傳熱，與實際工藝元件由“芯-管壁-土壤”的熱量傳導過程不一致，故在加熱前期溫度模擬數據與試驗數據會存在較大偏差，當模擬加熱元件與試驗加熱元件的管壁受熱傳熱穩定后，二者的差距逐漸變小. 整體而言，該模型能夠較好地模擬有機污染場地溫度場的時空分布特征和污染物的動態去除規律.

4 結論

a) 經過7周熱強化氣相抽提系統的持續運行，加熱井間距為3 m的TCH-A區土壤平均溫度達到99.4℃，加熱井間距為4 m的TCH-B區為84.8 ℃. TMVOC模型可很好地模擬中試規模的溫度場分布特征及溫升規律，擬合優度值高于0.989. 通過TMVOC模型計算，可以有效確定加熱井的合理布設間距以及目標溫度的最佳維持時間.

b) 氯苯的去除速率呈先增后降趨勢，第1周的去除速率最慢，為0.07 kg/d，第4周的去除速率最快，為0.46 kg/d. 熱強化氣相抽提8周后氯苯總剩余量為0.36 kg，去除率達97.3%. TMVOC模型可較好地模擬污染場地中氯苯的動態去除過程，擬合優度值達0.914.