基于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.