污染場(chǎng)地多相抽提低碳修復(fù)布井方案優(yōu)化模擬研究

金嘉路，楊逸文，王 震，張婉瑩，周 龑，張 峰，李 磊，崔長(zhǎng)征

1. 華東理工大學(xué)資源與環(huán)境工程學(xué)院，國(guó)家環(huán)境保護(hù)化工過(guò)程環(huán)境風(fēng)險(xiǎn)評(píng)價(jià)與控制重點(diǎn)實(shí)驗(yàn)室，上海市環(huán)境保護(hù)化學(xué)污染物環(huán)境標(biāo)準(zhǔn)與風(fēng)險(xiǎn)管理重點(diǎn)實(shí)驗(yàn)室，上海 200237

2. 上海市環(huán)境科學(xué)研究院，國(guó)家環(huán)境保護(hù)城市土壤污染控制與修復(fù)工程技術(shù)中心，上海 200233

3. 中石化第五建設(shè)有限公司，廣東 廣州 510145

4. 上海格林曼環(huán)境技術(shù)有限公司，上海 200001

5. 中石化煉化工程（集團(tuán)）股份有限公司洛陽(yáng)技術(shù)研發(fā)中心，河南 洛陽(yáng) 471003

多相抽提(multi-phase extraction，MPE)是一種通過(guò)“抽提井”真空提取污染區(qū)域地下的氣體、地下水和非水相液體(non-aqueous phase liquids，NAPL)，將污染物分離和處理的原位場(chǎng)地修復(fù)技術(shù)，作為氣相抽提(soil vapor extraction，SVE)的改進(jìn)技術(shù)已被廣泛應(yīng)用于實(shí)踐，可以與空氣注入與原位熱脫附等其他技術(shù)聯(lián)合應(yīng)用強(qiáng)化修復(fù)效果[1-3]. MPE 修復(fù)對(duì)于中高滲透性的土壤具有較好的修復(fù)效果，在我國(guó)，張峰[4]在華東地區(qū)某工廠甲苯污染砂質(zhì)粉土場(chǎng)地進(jìn)行了MPE修復(fù)中試試驗(yàn)，在25 d 內(nèi)總共去除甲苯污染物約125 kg，占總?cè)コ齆APL 污染物的85%，部分監(jiān)測(cè)點(diǎn)修復(fù)后未檢出NAPL 殘余. 張晶等[5]在某有機(jī)復(fù)合污染砂質(zhì)粉土場(chǎng)地應(yīng)用MPE 和原位化學(xué)氧化技術(shù)進(jìn)行了聯(lián)合修復(fù)，45 d 的修復(fù)周期內(nèi)收集了約100 L 的LNAPL，污染物去除率大于90%.

修復(fù)井間距對(duì)場(chǎng)地修復(fù)具有重要意義，NAPL 污染場(chǎng)地修復(fù)常存在拖尾與反彈現(xiàn)象[6]，Bass 等[7]對(duì)44 個(gè)實(shí)際場(chǎng)地的修復(fù)案例的調(diào)查結(jié)果顯示，當(dāng)單井所需修復(fù)的半徑小于5 m 時(shí)，往往能無(wú)反彈修復(fù)，而所需覆蓋半徑超過(guò)7 m 時(shí)，修復(fù)往往失敗. 現(xiàn)有的研究顯示，土壤滲透率等水文地質(zhì)因素對(duì)抽提井影響半徑和修復(fù)效率有著極大的影響[8-10]，然而，目前MPE現(xiàn)場(chǎng)實(shí)踐通常依據(jù)工程經(jīng)驗(yàn)均勻布設(shè)抽提井，沒有針對(duì)場(chǎng)地滲透率等水文地質(zhì)條件及污染羽分布設(shè)計(jì)布井方案，常采用密集的梅花樁布點(diǎn)法布設(shè)抽提井，導(dǎo)致MPE 修復(fù)效率低、成本高. 因此，如何根據(jù)場(chǎng)地情況優(yōu)化MPE 布井方案和工藝，降低修復(fù)成本，提升抽提效率，降低碳排放量備受關(guān)注[11-12].

由于實(shí)際場(chǎng)地的復(fù)雜性和不確定性，單獨(dú)探尋各滲透率等地質(zhì)參數(shù)下的抽提井影響半徑?jīng)]有實(shí)際意義[11,13]，近年來(lái)，越來(lái)越多的研究者通過(guò)TMVOC、GMS、MT3DMS 等商業(yè)軟件及模型，利用數(shù)值模擬研究土壤地下水修復(fù)過(guò)程，優(yōu)化布井方案和操作工藝[14-17]. 這一過(guò)程中，部分研究者通過(guò)模擬軟件刻畫實(shí)際污染場(chǎng)地模型，經(jīng)人工調(diào)參的方式，修改模擬參數(shù)，比選最佳修復(fù)工藝和布井方案[18-19]. 但人工調(diào)參試錯(cuò)法機(jī)械重復(fù)性強(qiáng)，往往只能涉及50~100 次的工作量，為減少人工試參的過(guò)程，也有研究者基于模擬軟件的模擬結(jié)果，建立簡(jiǎn)化替代模型，利用MATLAB等數(shù)據(jù)處理軟件結(jié)合優(yōu)化算法計(jì)算最優(yōu)方案[20]，但替代模型的精確度會(huì)影響優(yōu)化結(jié)果的正確性. 為此，部分研究者直接利用外部程序調(diào)用地下水模擬軟件，結(jié)合優(yōu)化算法計(jì)算最優(yōu)布井方案和操作工藝，其中基于MT3DMS 開發(fā)的模塊化地下水優(yōu)化軟件(modular groundwater optimizer，MGO)被廣泛應(yīng)用在地下水抽出問題中[21]，該方法在保證模擬準(zhǔn)確性的同時(shí)也免除了人工調(diào)參過(guò)程. 耿國(guó)婷等[22]利用MGO 模擬三氯乙烯污染含水層的抽水修復(fù)過(guò)程，利用遺傳算法優(yōu)化抽提井布井位置及抽水量等工藝參數(shù)，減少了18%的運(yùn)行成本，提升了經(jīng)濟(jì)效益. 但MGO 基于的模擬軟件MT3DMS 主要適用于地下水抽出修復(fù)，其修復(fù)區(qū)域主要為地下水含水層，而MPE 過(guò)程涉及氣、液、NAPL三相變化，修復(fù)區(qū)域包含包氣帶和含水層，現(xiàn)有研究普遍采用針對(duì)多相運(yùn)移問題開發(fā)的TMVOC 軟件模擬MPE 修復(fù)過(guò)程[9,23-24].

TMVOC 是由美國(guó)伯克利實(shí)驗(yàn)室開發(fā)的基于質(zhì)能守恒、亨利定律等物理化學(xué)理論規(guī)律模擬多相流體在飽和-非飽和帶運(yùn)移的模擬軟件，能夠有效模擬抽提修復(fù)過(guò)程，反應(yīng)修復(fù)過(guò)程中污染物濃度等參數(shù)變化. 孫超等[25]和王穎[9]分別通過(guò)TMVOC 模擬了場(chǎng)地中熱強(qiáng)化SVE 修復(fù)過(guò)程和苯系物在一維土柱的垂向遷移過(guò)程，TMVOC 模擬結(jié)果與試驗(yàn)數(shù)據(jù)的擬合優(yōu)度(R)分別為0.995 和0.89. 在對(duì)場(chǎng)地模擬還原的基礎(chǔ)上，研究者基于TMVOC 模型探究了抽提修復(fù)過(guò)程中的主要影響因素與最優(yōu)參數(shù)，王穎[9]通過(guò)TMVOC模擬對(duì)15 個(gè)場(chǎng)地的SVE 修復(fù)效果的影響因素進(jìn)行了灰色關(guān)聯(lián)度分析，結(jié)果表明土壤滲透率的關(guān)聯(lián)度最高，為0.841. 田蕾等[23]通過(guò)TMVOC 對(duì)苯和乙烯污染的低滲透場(chǎng)地模型的MPE 修復(fù)過(guò)程進(jìn)行了模擬，探究抽提井開篩位置的最優(yōu)設(shè)計(jì)參數(shù)，結(jié)果顯示，抽提井開篩最優(yōu)位置應(yīng)處于含水層深部. 近年來(lái)，已有研究者開始嘗試將TMVOC 與優(yōu)化算法相結(jié)合，如Sookhak 等[11]利用TMVOC-MP 與遺傳算法相結(jié)合，優(yōu)化了MPE 修復(fù)場(chǎng)地中抽提井水平布井方案，指出均勻布井方案并非最佳修復(fù)方案. 現(xiàn)有的優(yōu)化研究尚缺乏對(duì)污染場(chǎng)地內(nèi)每口抽提井開篩位置的優(yōu)化探討，實(shí)際修復(fù)過(guò)程中，抽提井開篩位置位于非飽和區(qū)時(shí)修復(fù)以氣相污染物抽取為主，而位于飽和區(qū)時(shí)修復(fù)以地下水和自由相NAPL 抽取為主，統(tǒng)一的抽提井的開篩位置難以達(dá)到最優(yōu)抽提效果[26]，針對(duì)開篩位置的優(yōu)化修復(fù)方案具有必要性.

因此，本研究利用MATLAB 實(shí)現(xiàn)遺傳算法和TMVOC 軟件調(diào)用，構(gòu)建了基于場(chǎng)地污染羽分布和土壤滲透率、有機(jī)碳含量等水文地質(zhì)條件優(yōu)化MPE 修復(fù)方案的MATLAB-TMVOC 聯(lián)合優(yōu)化程序，在優(yōu)化布井水平位置的同時(shí)優(yōu)化各抽提井的開篩位置，并探究了不同數(shù)量井群的最優(yōu)布設(shè)及開篩位置，總結(jié)了不同土壤滲透率下的最優(yōu)布井規(guī)律，以期為實(shí)際現(xiàn)場(chǎng)抽提井布設(shè)提供指導(dǎo)，實(shí)現(xiàn)高效和低碳修復(fù).

1 研究方法

1.1 MPE 模擬與優(yōu)化軟件

本研究所涉及的軟件包括數(shù)據(jù)處理軟件MATLAB(R2021a，MathWorks，Inc.，美國(guó))、地下水模擬軟件TMVOC(TOUGH2-TMVOC， Lawrence Berkeley National Laboratory，美國(guó))和地下水模擬及可視化軟件PetraSim(Version 5.1，RockWare，Inc.，美國(guó)).

1.2 MPE 修復(fù)模型設(shè)置

根據(jù)現(xiàn)有場(chǎng)地尺度模擬研究[9,11,18,23]，本研究以苯作為唯一污染物，建立了二維切片概念模型，模型尺寸及網(wǎng)格劃分如圖1 所示. 模型水平方向上長(zhǎng)度為100 m，豎直方向上深度為15 m，厚度為1 m. 此外，模型上部設(shè)置高度為0.001 m 的網(wǎng)格作為大氣邊界，左右兩側(cè)設(shè)置水平長(zhǎng)度為0.001 m 的網(wǎng)格作為定壓水頭邊界. 模型網(wǎng)格均勻劃分，在水平方向上以2 m 為單位長(zhǎng)度劃分網(wǎng)格，以最左側(cè)作為相對(duì)水平距離等于0 m的原點(diǎn)；在豎直方向上以1 m 為單位長(zhǎng)度劃分網(wǎng)格，將地下水埋深位置作為相對(duì)深度等于0 m 的原點(diǎn). 模型頂部相對(duì)深度為5 m，表示其位于地下水位線以上5 m；模型底部相對(duì)深度為——10 m，表示其位于地下水位線以下10 m. 模型網(wǎng)格總數(shù)為832 個(gè)，按照由左至右、由下往上的方式逐行編號(hào). 在場(chǎng)地內(nèi)進(jìn)行為期1 年的釋放以模擬污染物的泄漏，并進(jìn)行1 年的污染物遷移擴(kuò)散模擬以設(shè)置污染物初始分布，并模擬了為期90 d的MPE 修復(fù).

圖1 模型尺寸與網(wǎng)格劃分Fig.1 Model size and grid division

模型巖性設(shè)置中采用修正后的Stone 模型[27]描述場(chǎng)地相對(duì)滲透率，其方程形式如式(1)~(3)所示：

式中，krg、krw、krn分別表示氣相、液相和NAPL 相的相對(duì)滲透率，krn′為修正后的近似NAPL 相相對(duì)滲透率，Sg、Sw、Sn分別為氣相、液相和NAPL 相飽和度，Sgr、Swr、Snr分別為氣相、液相和NAPL 相殘余飽和度，A為擬合參數(shù). 其中TMVOC 所需輸入?yún)?shù)為Sgr、Swr、Snr以及A，本研究依照經(jīng)驗(yàn)值和TMVOC 手冊(cè)[23,28]分別取0.01、0.2、0.05 和3.

采用Parker 模型[29]描述三相系統(tǒng)的毛細(xì)壓力，具體公式如式(4)~(7)所示：

式中：Sw、Sm、Sn分別為液相、殘余相和NAPL 相飽和度；和分別為由式(4)和式(5)計(jì)算得到毛細(xì)管壓力函數(shù)的斜率；Pcgn和Pcgw分別為氣相-NAPL 相間毛細(xì)壓力以及氣相-水相間毛細(xì)壓力，Pa；g為重力加速度，m/s2；αnw和αgn分別為氣相-NAPL 相之間進(jìn)氣壓力的倒數(shù)，以及NAPL 相和水相之間進(jìn)氣壓力的倒數(shù)，1/Pa；B為Parker 模型的擬合參數(shù)；B′為計(jì)算過(guò)程的中間變量，數(shù)值上B′=1——1/B. 其中TMVOC 所需輸入?yún)?shù)為Sm、B、αnw和αgn，本研究依照經(jīng)驗(yàn)值和TMVOC 手冊(cè)[23,28]分別取0.00、1.84、11.0 和10.0.

模型的水文地質(zhì)參數(shù)、污染物設(shè)置參數(shù)與修復(fù)工藝參數(shù)取值范圍參考了現(xiàn)有的TMVOC 研究[9,18,23,30-32]和TMVOC 手冊(cè)[28]中案例所選取的參數(shù)范圍以及工程經(jīng)驗(yàn). 對(duì)于巖性參數(shù)，現(xiàn)有文獻(xiàn)報(bào)道[30-32]顯示，土壤滲透率參數(shù)取值介于10——14~10——10m2之間，垂直滲透率為水平滲透率的0.25 倍；孔隙度隨土壤滲透率的減小而增大，一般而言，黏土孔隙度為0.4~0.6，砂土孔隙度為0.3~0.4；巖層密度普遍采用2 600~2 700 kg/m3；常見的土壤其有機(jī)碳含量為0%~5%[33]. 其他參數(shù)中，抽提井真空度普遍取1×104~3×104Pa[9]，污染物釋放速率普遍為10——6~10——5kg/s，降雨入滲設(shè)置可按照TMVOC手冊(cè)[28].

本研究在上述范圍內(nèi)依照工程經(jīng)驗(yàn)設(shè)定了模型參數(shù)，具體數(shù)值如表1 所示. 此外，為探究土壤滲透率對(duì)抽提修復(fù)和最優(yōu)抽提井布設(shè)方案的影響，在此模型的基礎(chǔ)上設(shè)定了8 種污染羽分布相同、土壤滲透率不同的場(chǎng)景模型，其水平滲透率參數(shù)取值分別為2×10——13、4×10——13、6×10——13、8×10——13、2×10——12、4×10——12、6×10——12和8×10——12m2，其余參數(shù)同表1.

表1 模型中水文地質(zhì)參數(shù)取值與修復(fù)工藝參數(shù)Table 1 Hydrogeological and remediation process parameters in the model

1.3 MPE 布井及開篩位置優(yōu)化算法參數(shù)設(shè)置及約束條件

MPE 布井及開篩位置優(yōu)化算法采用遺傳算法，目標(biāo)函數(shù)為抽提修復(fù)90 d 后的污染物去除率，分別以錦標(biāo)賽選擇法、單點(diǎn)變異和交換法實(shí)現(xiàn)種群的選擇、變異和交換，群數(shù)大小設(shè)置為150，代溝為0.8，變異概率為0.01，迭代25 次[34].

TMVOC 所遵循質(zhì)能平衡方程的形式如式(8)[28,35]所示：

式中：Va為流動(dòng)單元體的體積，m3；Γa為表面積，m2；Mk為組分k在單位土壤介質(zhì)中的質(zhì)量，kg；Fk為進(jìn)入流動(dòng)單元體中組分k的總通量，kg；qk為組分k在單元體的源匯項(xiàng)，kg；a為流動(dòng)區(qū)單元體表面的外法向單位矢量.

TMVOC 以GENER 模塊描述特殊源匯項(xiàng)，其中同處一列的DELV 屬性網(wǎng)格表示抽提井，因此設(shè)置抽提井位置時(shí)應(yīng)當(dāng)避開模型兩側(cè)及上方邊界，故布井位置和開篩深度的約束條件可見式(9)~(12)：式中：ωc,i,j為第i口抽提井的開篩位置底部上方第j個(gè)網(wǎng)格在模型中的列位置；ωr,i,j為該網(wǎng)格在模型中的行位置；X和Z分別為模型總列數(shù)與總行數(shù)；L為設(shè)定的抽提井開篩長(zhǎng)度對(duì)應(yīng)的網(wǎng)格行數(shù)，本研究開篩長(zhǎng)度為5 m，取5 行；Y為降雨入滲屬性所占網(wǎng)格行數(shù)，數(shù)值取1.

1.4 經(jīng)濟(jì)性評(píng)估與碳排放量核算

本研究對(duì)成本計(jì)算和碳排放量核算方法進(jìn)行了適當(dāng)簡(jiǎn)化[13]，以定性分析各修復(fù)方案的成本和碳排放量. 以抽提井?dāng)?shù)量和深度評(píng)價(jià)施工成本，以抽提井運(yùn)行過(guò)程中的電能耗為唯一因素核算運(yùn)行成本及碳排放量[36]，計(jì)算方法如式(13)~(15)[13]所示，為評(píng)定最優(yōu)修復(fù)方案，依據(jù)經(jīng)驗(yàn)簡(jiǎn)單設(shè)定了式中變量的參考值以量化修復(fù)成本和碳排放量，在實(shí)際修復(fù)過(guò)程中，應(yīng)當(dāng)結(jié)合具體情況設(shè)定其數(shù)值.

式中：N為抽提井?dāng)?shù)量；θbu為建井施工成本，元；Di為第i口抽提井的鉆井深度，m；kd為鉆井深度費(fèi)用系數(shù)，元/m，取經(jīng)驗(yàn)值250 元/m；θbu′為其他施工成本，元，取1 000 元；θru為運(yùn)行成本，元；Ei為第i口抽提井的每日能耗(kW·h/d)，與抽提井真空度和流量相關(guān)，為簡(jiǎn)化研究，假設(shè)每口抽提井每日能耗相同，統(tǒng)一取60 kW·h/d；θru′為其他每日運(yùn)行成本，取1 000 元；T為達(dá)到修復(fù)目標(biāo)值的時(shí)間，d；kE為能耗價(jià)格系數(shù)，即單位電價(jià)，以1 元/(kW·h)計(jì)；C為碳排放量，kg；kC為單位能耗碳排放量系數(shù)，kg/(kW·h)，取0.960 kg/(kW·h)[36].

2 結(jié)果與討論

2.1 MPE 布井方案及開篩位置優(yōu)化算法構(gòu)建

利用MATLAB 搭建了基于TMVOC 與遺傳算法的聯(lián)合優(yōu)化程序，利用程序代替了人工調(diào)參試錯(cuò)過(guò)程，輸入文件改寫和讀取輸出文件所用時(shí)間僅為0.8 s，相較人工設(shè)置參數(shù)可忽略不計(jì)；遺傳算法的運(yùn)用相較人工優(yōu)化可以在有限的次數(shù)中確定最優(yōu)布井位置. 其主要步驟包括建立輸入文件與讀取模型數(shù)據(jù)、參數(shù)集編碼、修改輸入文件、運(yùn)行TMVOC 并讀取模擬結(jié)果、優(yōu)化結(jié)果計(jì)算.

2.1.1 建立輸入文件與讀取

使用MATLAB-TMVOC 優(yōu)化程序前，先利用PetraSim 快速建立場(chǎng)地水文地質(zhì)和污染物分布模型，基于該軟件獲得初始模型的TMVOC 輸入文件，MATLAB-TMVOC 優(yōu)化程序可以通過(guò)對(duì)該輸入文件的改寫和模擬實(shí)現(xiàn)優(yōu)化過(guò)程.

根據(jù)TMVOC 用戶手冊(cè)，抽提井網(wǎng)格的參數(shù)被編寫在輸入文件的“GENER”(源和匯)模塊中，該模塊中每一行為一個(gè)網(wǎng)格所具有的源匯屬性，共6 列，其中第1 列內(nèi)容為對(duì)應(yīng)的網(wǎng)格編號(hào)，同一列連續(xù)的網(wǎng)格被視作一口抽提井；第4 列為屬性，其中DELV 為抽提井，COM1 常用作記錄降雨入滲，COM2 則被設(shè)置為空氣注入井；第5 列、第6 列分別為流量(曝氣)和壓力(抽提).

利用程序逐行讀取文件后判定該行內(nèi)容中是否存在“COM1、COM2、DELV”這三類TMVOC 中“GENER”模塊的“源與匯”屬性，最終記錄下“GENER”模塊的所有參數(shù)，計(jì)算模型中井?dāng)?shù)量以及網(wǎng)格數(shù)量.程序巧妙地利用了“COM1”入滲網(wǎng)格占據(jù)模型中一整行且上方僅存一行大氣邊界網(wǎng)格的原理，計(jì)算出了模型的總網(wǎng)格數(shù)、列數(shù)、行數(shù). 總網(wǎng)格數(shù)計(jì)算公式：

式中，X為模型總列數(shù)， ?MAX與 ?MIN分 別 表 示 具 有“COM1”屬性的網(wǎng)格的最大網(wǎng)格編號(hào)和最小網(wǎng)格編號(hào).

因此總網(wǎng)格數(shù)為?MAX與α之和，行數(shù)為總網(wǎng)格數(shù)與列數(shù)之商.

2.1.2 參數(shù)集編碼

井網(wǎng)格的參數(shù)集確定及編碼對(duì)遺傳算法的運(yùn)行至關(guān)重要. 依照TMVOC 將所處同一列連續(xù)的“DELV”屬性網(wǎng)格判定為一口抽提井的原則，井網(wǎng)格參數(shù)集可以轉(zhuǎn)化為每口抽提井最底部的井網(wǎng)格列編碼和行編碼.

顯然對(duì)于井深度可以通過(guò)單井范圍給定參數(shù)集，而對(duì)于井水平位置，可以采用“nchoosek”函數(shù)窮舉所有列組合，然后通過(guò)選擇組合的編號(hào)確定每口抽提井所在列.

2.1.3 修改輸入文件、運(yùn)行TMVOC 并讀取模擬結(jié)果

在獲得初始二進(jìn)制編碼轉(zhuǎn)換的輸入?yún)?shù)集后，修改輸入文件，建立新的工藝參數(shù)下的修復(fù)模型，通過(guò)“system”函數(shù)調(diào)用TMVOC 便可以求解該模型，生成輸出文件，讀取輸出文件中“TOTAL VOCs”模塊內(nèi)容的污染物殘余量評(píng)價(jià)該參數(shù)的適應(yīng)度.

2.1.4 優(yōu)化結(jié)果計(jì)算

基于遺傳算法對(duì)于不同的布井方案進(jìn)行模擬以及優(yōu)化迭代，最終確定最優(yōu)布井位置. 利用MATLABTMVOC 優(yōu)化布井方案及開篩位置可以有效地減少模擬時(shí)間和運(yùn)算次數(shù)，優(yōu)化程序可以自動(dòng)化修改輸入文件和讀取模擬結(jié)果，所用時(shí)間可忽略不計(jì)，進(jìn)行單次模擬的時(shí)長(zhǎng)僅為TMVOC 模擬運(yùn)行時(shí)長(zhǎng)，約為70 s，而對(duì)于相同的模型，通過(guò)人工設(shè)置參數(shù)進(jìn)行模擬的過(guò)程時(shí)長(zhǎng)為120~180 s，自動(dòng)化運(yùn)行可以節(jié)省75%以上的單次模擬時(shí)間. 整個(gè)優(yōu)化過(guò)程中，人工調(diào)參優(yōu)化需要模擬嘗試的次數(shù)具有不確定性，而遺傳算法通過(guò)隨機(jī)搜索可以在有限的次數(shù)中尋找到最優(yōu)布井位置，快速高效的處理速度更有利于抽提井布設(shè)方案與開篩位置的優(yōu)化.

在進(jìn)行抽提井布設(shè)方案優(yōu)化探究時(shí)，布設(shè)方案眾多，本模型水平方向可選網(wǎng)格位置為50 種，最小布井間距為單元格水平距離2 m，豎直方向可選網(wǎng)格位置為9 種. 當(dāng)布設(shè)5 口抽提井時(shí)，僅在水平方向有2 118 760〔C(50,5)〕種組合可能，但對(duì)于這種簡(jiǎn)單組合問題，遺傳算法所需種群大小遠(yuǎn)小于參數(shù)集大小[11]，當(dāng)遺傳算法種群數(shù)為150 時(shí)，25 次迭代的運(yùn)算過(guò)程中即可收斂獲得最優(yōu)解，總模擬次數(shù)為3 150 次，相較窮舉所有可能性，計(jì)算量簡(jiǎn)化率超過(guò)90%.

2.2 不同MPE 抽提井?dāng)?shù)量方案的抽提井布點(diǎn)及開篩位置優(yōu)化

通過(guò)PetraSim 模擬了污染物遷移釋放及抽提修復(fù)過(guò)程并繪制了污染羽分布圖，各井群數(shù)量下布井方案及開篩位置優(yōu)化前后污染羽分布情況及去除率時(shí)間曲線見圖2，在此基礎(chǔ)上通過(guò)MATLAB-TMVOC聯(lián)合程序優(yōu)化了抽提井水平布置和開篩位置，經(jīng)優(yōu)化后的抽提方案可以于65 d 內(nèi)達(dá)到90%的污染物去除率.

圖2 優(yōu)化前后修復(fù)60 d 后殘余污染物分布及90 d 內(nèi)污染物去除率變化Fig.2 Distribution of residual pollutants after 60 d and changes in removal rate within 90 d of remediation before and after optimization

均勻布井修復(fù)前后污染與分布和去除率情況如圖2(a)(b)所示. 修復(fù)前的初始污染羽水平長(zhǎng)度約為58 m，豎直長(zhǎng)度約為6 m，超過(guò)50%的污染物處在模型相對(duì)水平位置17~34 m 處，即污染源附近形成了一個(gè)污染核心區(qū)域，該區(qū)域內(nèi)地下水中的污染物濃度超過(guò)其他區(qū)域2 倍. 本研究設(shè)置了不同數(shù)量的抽提井群，分別以24、16、12、10、8 m 的間隔均勻布設(shè)抽提井以保證污染物的去除，抽提井開篩位置根據(jù)核心區(qū)域污染物分布情況設(shè)置，篩管底部的相對(duì)深度設(shè)置為——4 m，篩管頂部的相對(duì)深度設(shè)置為1 m，具體布井點(diǎn)位如表2 所示. 在90 d 的修復(fù)周期內(nèi)，隨著布井?dāng)?shù)量由3 口增至7 口，污染物去除率由58.5%增至89.2%，修復(fù)效率增加了52.4%. 均勻布設(shè)間距為8 m 的7 口抽提井可以在93 d 達(dá)到90%的去除率.

表2 優(yōu)化前后的抽提井位置布設(shè)Table 2 Layout of extraction well positions before and after optimization

利用MATLAB-TMVOC 聯(lián)合優(yōu)化程序優(yōu)化了不同數(shù)量井群的布井方案和每口井的開篩位置，其中開篩長(zhǎng)度統(tǒng)一為5 m，優(yōu)化后修復(fù)前后殘余污染物的分布和去除率情況如圖2(c)(d)所示，在污染物去除率上，90 d 的周期內(nèi)，7 口抽提井MPE 修復(fù)方案的污染物去除率可達(dá)97.0%，相比優(yōu)化前提高了7.9%，優(yōu)化前7 口抽提井需要93 d 才能達(dá)到90%的污染物去除率，而優(yōu)化后7 口抽提井在63 d 內(nèi)即可達(dá)到目標(biāo)修復(fù)值，縮短了30 d 的修復(fù)時(shí)間，優(yōu)化方案提升了32.9%的修復(fù)速率.

優(yōu)化后布井方案及開篩位置如表2 所示，在布井規(guī)律上，相較均勻布井，優(yōu)化后的抽提井排布是根據(jù)污染物濃度及抽提井?dāng)?shù)量設(shè)計(jì)的非等間距布點(diǎn)，隨著污染羽的分布，在水流方向上由密至疏，抽提井的數(shù)量在4 口以上時(shí)，位于核心區(qū)域的兩口抽提井間距縮小至2~4 m，而污染羽尾部低濃度區(qū)抽提井間距為10~20 m，這與王穎[9]對(duì)于多個(gè)低滲透場(chǎng)地SVE 模擬修復(fù)結(jié)果相近，當(dāng)抽提真空度為10 kPa 時(shí)，單井的影響半徑為3.5~16 m. Sookhak 等[11]利用遺傳算法與TMVOC 的并行版本TMVOC-MP 相結(jié)合研究了多相抽提修復(fù)某低滲透性污染場(chǎng)地的平面布井模型，得到了相同的結(jié)論，均勻布井并非是最優(yōu)布井方案. 當(dāng)抽提井開篩長(zhǎng)度固定為5 m 時(shí)，最優(yōu)布井方案的開篩位置在水平方向上深度逐漸增加，在污染羽核心區(qū)域，即污染羽前段，篩管頂部相對(duì)深度為1 m，略低于污染羽，修復(fù)以多相抽提為主，隨著污染羽核心區(qū)域水平距離的增加，篩管頂部相對(duì)深度為0 m，以液相抽提為主，這是因?yàn)槲廴疚镏饕ㄟ^(guò)地下水流動(dòng)而發(fā)生橫向遷移，下游區(qū)域主要為液相LNAPL[37-38]. 田蕾等[23]探討了低滲透場(chǎng)地中抽提井開篩位置深度對(duì)于污染物去除率的影響，得到了相同的結(jié)論，當(dāng)抽提井開篩位置最底端降至NAPL 附近時(shí)，修復(fù)效率最高，達(dá)到90.4%，而僅對(duì)包氣帶進(jìn)行抽提時(shí)，去除率僅為29.1%. de Souza 等[17]通過(guò)STOMP 探究了MPE 修復(fù)過(guò)程中抽提井底部的最佳位置，結(jié)果顯示，當(dāng)開篩位置底端高程從14 m 提升至17 m 的過(guò)程中，最佳修復(fù)位置位于15.8 m 處，為氣相-NAPL 與液相-NAPL 交界面之間. 這些研究結(jié)果都顯示，MPE 抽提井開篩位置最底端應(yīng)當(dāng)?shù)陀诘叵滤唬瑢?duì)于NAPL 和地下水直接的抽提有利于污染物的去除. 4 種優(yōu)化布井方案都選擇相對(duì)水平距離25 m 處(污染物釋放點(diǎn))布設(shè)篩管底部深度為地水位線下3~4 m 的抽提井，可見此處是污染物MPE 修復(fù)的關(guān)鍵布井點(diǎn)，在實(shí)際場(chǎng)地布設(shè)抽提井時(shí)，應(yīng)當(dāng)重點(diǎn)關(guān)注污染源，布設(shè)抽提井群中的主井.

2.3 不同MPE 抽提井?dāng)?shù)量修復(fù)方案的修復(fù)成本和碳排放量核算

當(dāng)污染物目標(biāo)去除率為90%時(shí)，不同抽提井?dāng)?shù)量下優(yōu)化前后MPE 修復(fù)成本及碳排放量核算結(jié)果如表3 所示，優(yōu)化前，6 口抽提井修復(fù)運(yùn)行成本、總成本及總碳排放量(能耗)最低，分別為12.9×104元、14.8×104元與32.9 t，而5 口抽提井的總碳排放量最高，能耗最高. 優(yōu)化后6 口抽提井運(yùn)行成本最低，為8.4×104元；5 口井的總成本與碳排放量最低，分別為10.2×104元與19.6 t，盲目地縮小布井間距、增加抽提井?dāng)?shù)量會(huì)導(dǎo)致修復(fù)效率低，運(yùn)行成本高. 相比優(yōu)化前，優(yōu)化后的方案可以節(jié)約29.7%~34.7%的運(yùn)行成本、22.7%~29.3%的總成本以及21.4%~40.4%的碳排放量，MATLAB-TMVOC 聯(lián)合優(yōu)化后的抽提方案能夠有效提高抽提效率，具有綠色、低碳和高效的特點(diǎn).

2.4 不同滲透率場(chǎng)地中6 口MPE 抽提井修復(fù)方案的布井位置優(yōu)化

2.2 節(jié)分析結(jié)果顯示90 d 內(nèi)6 口抽提井MPE 修復(fù)方案對(duì)土壤滲透率為4×10——12m2的模擬場(chǎng)地的污染物去除率可達(dá)90%，因此本節(jié)基于1.2 節(jié)構(gòu)建的污染場(chǎng)地模型和污染羽分布，模擬了不同滲透率場(chǎng)地下6 口抽提井的MPE 修復(fù)過(guò)程，并利用MATLAB-TMVOC優(yōu)化程序探究最優(yōu)抽提井布井方案規(guī)律.

優(yōu)化前后不同土壤滲透率下6 口抽提井布井方案的布井點(diǎn)位水平位置和開篩位置如表4 所示. 優(yōu)化前的布井方案為均勻布設(shè)，具體點(diǎn)位同2.2 節(jié)中6 口抽提井均布方案，根據(jù)污染羽分布情況，以10 m 為間隔設(shè)置了均勻布井方案，篩管底部深度位于地下水位線下4 m；優(yōu)化后的布井方案為非均勻的布井方案，與2.2 節(jié)相似，位于核心區(qū)域的抽提井間距較短，為2~4 m，而污染羽尾部低濃度區(qū)抽提井間距較長(zhǎng)，為10~20 m. 優(yōu)化布井方案隨土壤滲透率的降低而發(fā)生改變，在土壤滲透率由8×10——12m2降至2×10——13m2的過(guò)程中，6 口抽提井的布井位置逐漸向污染物核心區(qū)域移動(dòng)，6 口井首尾距離由71 m 縮短至44 m，平均井間距由11.8 m 降至7.3 m，降低了38%. 4 口靠近污染羽核心區(qū)域(水平位置15~50 m)的抽提井的平均間距由10.5 m 降至5.5 m，減少了47.6%. 優(yōu)化后的布井平均間距與土壤滲透率呈正相關(guān)，當(dāng)土壤滲透率處在10——12~10——13m2范圍內(nèi)，土壤滲透率每減少10——13m2，優(yōu)化后的抽提井間距會(huì)減小約0.13 m，即當(dāng)土壤滲透率下降一個(gè)數(shù)量級(jí)時(shí)，抽提井間距應(yīng)減少17.8%~48.3%.

表4 優(yōu)化前后不同滲透率下的抽提井位置布設(shè)Table 4 Layout of extraction well positions for different permeability rates before and after optimization

優(yōu)化前后不同土壤滲透率下6 口抽提井的布井方案與修復(fù)60 d 后污染羽分布如圖3(a)(b)所示，污染羽分布變化顯示，污染物修復(fù)效果隨土壤滲透率的降低而減小，以土壤滲透率為8×10——12m2的模擬場(chǎng)地的均勻布井修復(fù)結(jié)果為例，經(jīng)60 d 的MPE 修復(fù)后污染羽水平相對(duì)位置由17~80 m 范圍內(nèi)縮小至20~30 m范圍內(nèi)，相對(duì)深度位于包氣帶的污染羽發(fā)生消失；而當(dāng)土壤滲透率為8×10——12m2時(shí)，經(jīng)均勻布井修復(fù)60 d后的污染羽大小幾乎不變，包氣帶仍存在大量污染物. 對(duì)比結(jié)果表明，優(yōu)化后的方案較優(yōu)化前可以縮小修復(fù)后的污染羽范圍，以土壤滲透率為8×10——12m2的場(chǎng)地為例，相同比例尺下，優(yōu)化前仍有少部分可見污染羽，而優(yōu)化后幾乎未見污染羽分布.

圖3 優(yōu)化前后不同滲透率下的布井位置及污染物去除率Fig.3 Well location and pollutant removal rates before and after optimization for different permeability conditions

優(yōu)化前后不同土壤滲透率下修復(fù)過(guò)程中苯去除率隨時(shí)間的變化曲線如圖3(c)(d)所示，隨著土壤滲透率由8×10——12m2降至2×10——13m2，均勻布井方案的污染物去除率由97.1%降至10.2%，而優(yōu)化布井方案的污染物去除率由99.2%降至21.5%，優(yōu)化后的布井方案可以提升污染物修復(fù)效果，各土壤滲透率下MPE修復(fù)污染物去除率大于均勻布井方案，土壤滲透率處于2×10——13m2~8×10——12m2范圍時(shí)，優(yōu)化后的修復(fù)方案可使污染物去除率增加2.1%~16.2%. 王穎[9]利用TMVOC模擬了15 個(gè)土壤滲透率介于10——14~10——12m2之間的場(chǎng)地模型中的苯系物泄漏、遷移以及SVE 修復(fù)過(guò)程，15 個(gè)模型除土壤滲透率外水文地質(zhì)參數(shù)相近，苯泄漏速率均設(shè)置為2×10——5kg/s，泄漏與擴(kuò)散時(shí)間均為1 年，在布井間距為7.5~10 m、抽提真空度為1×104Pa 的等間距布井方案下經(jīng)SVE 修復(fù)2 年，污染物去除率為9%~99%；其中包氣帶土壤滲透率為1×10——12m2與9×10——13m2的場(chǎng)地，經(jīng)SVE 修復(fù)后污染物去除率分別為60%與50%. 本研究中場(chǎng)地土壤滲透率為2×10——12m2與8×10——13m2的場(chǎng)景經(jīng)優(yōu)化后布井方案修復(fù)后，180 d 內(nèi)的污染物去除率分別達(dá)到90%與60%. 對(duì)比顯示，在土壤滲透率、污染物羽分布和抽提真空度相似的條件下，基于MATLAB-TMVOC 聯(lián)合優(yōu)化程序設(shè)計(jì)的方案可以顯著提升污染物的修復(fù)效率.

總體而言，當(dāng)土壤滲透率大于2×10——12m2時(shí)，最小布井距離為4 m，90 d 內(nèi)污染物去除率高于90%，當(dāng)土壤滲透率處于10——13m2數(shù)量級(jí)時(shí)，污染物去除率低于60%，最小布井間距達(dá)模型布井間距下限(2 m)，為達(dá)到90%的去除率，應(yīng)當(dāng)增加抽提井?dāng)?shù)量或修復(fù)時(shí)間. 這是因?yàn)橥寥罎B透率的降低會(huì)使流體流動(dòng)能力減弱，使得污染物向氣相、液相的遷移能力降低，從而使得抽提效率降低，抽提井有效影響半徑減小[39].

2.5 不確定性分析

本研究以數(shù)值模擬的方式探究MPE 修復(fù)抽提井布設(shè)的優(yōu)化方案，其結(jié)果具有不確定性，主要體現(xiàn)在兩個(gè)方面：①在優(yōu)化效果上，雖然2.2 節(jié)與2.4 節(jié)的模擬結(jié)果顯示，優(yōu)化后的方案相較優(yōu)化前的均勻布井方案，能夠提升2.1%~32.9%的修復(fù)效率，但是單一遺傳算法的應(yīng)用存在產(chǎn)生局部最優(yōu)解的可能，優(yōu)化后方案可能并非全局最優(yōu)布井方案，需要結(jié)合諸如改進(jìn)遺傳-模擬退火算法等聯(lián)合算法的應(yīng)用以降低遺傳算法發(fā)生“早熟”的可能[40]. ②研究模型僅為理想化的二維切片，在模型構(gòu)建和經(jīng)濟(jì)及碳排放評(píng)估過(guò)程中選取的部分參數(shù)存在理想化和簡(jiǎn)化，對(duì)計(jì)算結(jié)果存在影響. 薛成杰等[36]的研究結(jié)果顯示，使用氣相抽提等土壤原位修復(fù)技術(shù)修復(fù)1 m3土壤的成本為500~1 000元，修復(fù)1 t 土壤的碳排放量為40~400 kg，以此參數(shù)計(jì)算本場(chǎng)地的修復(fù)成本與碳排放量，可知2.2 節(jié)場(chǎng)地的修復(fù)成本為2×105~4×105元，碳排放量為23~230 t，該結(jié)果均存在一定程度的低估，主要原因在于，為簡(jiǎn)化分析，本研究?jī)H以抽提井運(yùn)行過(guò)程中的電能耗為唯一因素核算運(yùn)行成本及碳排放量. 但實(shí)際上多相抽提工藝中諸如管理等諸多環(huán)節(jié)都涉及碳減排潛能，特別是對(duì)于抽提尾氣處理環(huán)節(jié)，需要對(duì)收集抽提尾氣的活性炭進(jìn)行無(wú)害化處理，該過(guò)程也是MPE 修復(fù)中重要的碳排放源[36].

3 結(jié)論

a) 本研究建立的MATLAB-TMVOC 優(yōu)化程序，能夠有效聯(lián)合TMVOC 與優(yōu)化算法，自動(dòng)化實(shí)現(xiàn)MPE 修復(fù)模擬，免除數(shù)值模擬繁重的人工調(diào)參過(guò)程，縮短單次模擬的時(shí)間，在有限的計(jì)算模擬次數(shù)中高效求解最優(yōu)布井方案.

b) 多相抽提井布設(shè)間距應(yīng)當(dāng)隨著土壤滲透率的降低而縮小，當(dāng)土壤滲透率處在10——13~10——12m2范圍時(shí)，抽提井間距與滲透率成正比. 使用MPE 技術(shù)修復(fù)土壤滲透率低于10——13m2的黏性地層時(shí)，場(chǎng)地污染羽核心位置的布井間距應(yīng)當(dāng)小于4 m.

c) MPE 布井方案和抽提井開篩位置經(jīng)優(yōu)化后，可以減少抽提井?dāng)?shù)量并縮短修復(fù)時(shí)間，節(jié)約29.7%~34.7%的運(yùn)行成本，降低22.7%~29.3%的總成本，減少21.4%~40.4%的碳排放量，實(shí)現(xiàn)污染場(chǎng)地修復(fù)的低碳減排.