某氨氮污染地下水体抽出-处理系统优化模拟研究

王 平，韩占涛，张海领，孔贝贝，张鑫馨

(1.中国地质科学院水文地质环境地质研究所，河北 石家庄 050061；2.中国地质大学(北京)中国地质科学院，北京 100083；3.河北省&中国地质调查局地下水污染机理与修复重点实验室，河北 石家庄 050061；4.河北伟诚环境工程有限公司，河北 石家庄 050000；5.中国地质科学院，北京 100037)

近十几年来我国许多城市老工业区化工企业搬迁或停产[1]，企业搬迁后遗留的场地需要进行环境调查或修复才能转为其它用地[2]。针对污染场地地下水污染修复研究，我国目前多处于研究或中试阶段，相关工程实践不多[3]。目前常用的地下水污染修复技术分为三大类[4-6]：物理修复技术(帷幕法、抽出-处理技术等)[7]、化学修复技术(化学氧化与还原、可渗透反应墙等)[8-9]、生物修复技术(生物通气法、生物注射法、植物修复等)[10-13]。其中抽出-处理技术作为一种常用的地下水污染修复技术，在可溶性地下水污染物修复中应用广泛，其优点是简便易行，可直接移除地下水中污染物并控制污染物的扩散，对含水层破坏性低；缺点是由于污染物在含水层中的吸附/解吸反应以及水文地质条件复杂，导致修复时间长、处理费用昂贵[3,14]，若设计不合理还可能造成经济损失甚至水资源的巨大浪费。解决该问题的关键在于，工程中布设的井群系统能否高效地控制地下水污染水体的流动。这就需要了解污染羽随时间的变化，选择适当数量的抽水井，并设计合理的井位、井间距及抽水量。

利用数值模拟的手段，可以基于场地水文地质条件及特征资料，对污染物在地下水中的迁移进行模拟，设计抽水方案并模拟其修复效果，从中选出最优组合，为抽出-处理技术方案的设计提供依据[3]。早期抽出处理方案优化研究多为理论研究，其研究对象较为规则或理想化[15-16]。随后，采用GMS软件结合加权平均法优化抽出-处理系统的理论研究逐渐应用于野外场地[17]，如宫志强等[18]对某铬污染场地抽出处理修复中的抽水方式进行优化模拟，认为持续性抽出30%污染物后改用间歇性抽水方式，可降低修复成本。近年来，MODFLOW软件内嵌的MGO模块在抽出-处理方案优化中应用愈加广泛[19-20]：姜烈等[21]利用其进行抽水方案优化，认为最佳井位置应位于污染羽中游及下游的中轴线上，且中游抽水速率应比下游大；王燕[22]采用MGO对污染场地的抽出-处理井群进行优化，认为分阶段抽水模式比连续抽水模式更有效率。利用MGO进行抽出系统优化的特点是优化算法(包括遗传算法和模拟退火法)与模型的耦合，其优点是节省优化时间，可自动进行大量计算，从中选出最优方案；缺点是MGO模块只能调用稳定流模型，而自然界中水流的运动和溶质的迁移转化均随时间不断变化。另外，目前抽出-处理系统优化模拟研究中关于模型验证方面的论述较少，一般仅验证水流模型，少有溶质运移模型的进一步验证。而服务于地下水污染修复设计应多关注溶质运移方面的验证，以提高模型准确性，更好地预测污染羽的迁移转化、提高修复效果。

本研究对象为我国北方某化肥厂氨氮污染场地，为提高该场地高浓度氨氮污染地下水体的抽出处理修复效率，降低修复成本，选择广泛应用的GMS软件[23-24]的MODFLOW和MT3DMS模块进行地下水水流和溶质运移模拟，针对井数、抽水天数和总抽水量3个变量反复调整、试算，筛选出3种较好的抽水方案并进一步模拟比选，得出最优的抽水井布设位置及抽水方案，为抽出-处理技术方案设计提供依据。模型将非稳定流数值模拟应用于野外场地的溶质运移模拟，并进行了水流和多期监测井实测溶质浓度的拟合验证，提高模型准确度。

1 研究区概况

研究区位于河北省西北部，属温带大陆性季风气候、半干旱地区，年降水量为330～400 mm，主要集中在7—8月份，水资源严重不足。研究区地貌属浅山丘陵区，地势南高北低，场区位于山前到河流之间的冲洪积扇上，地势相对平坦。地下水类型为第四系松散岩类孔隙水，主要接受大气降水及侧向流入补给，地下水流向为西南向东北，水力梯度为0.000 5，实际流速约为0.013 m/d。调查区域北侧和东侧均为季节性河流，雨季充水、旱季断流。

该化肥厂污染场地总面积约0.24 km2，2016年化肥厂已依法破产并拆除。基于前期的环境调查及风险评估结果，场区地下水中主要污染物为氨氮，其最高浓度达578 mg/L；根据场区环境调查及风险评估报告，氨氮的修复目标值为64.9 mg/L。

结合场地地下水污染程度及其物理化学性质、修复时限要求以及场地施工条件限制等，经筛选得出比较适合本场地污染地下水的修复技术为抽出处理技术。为了设计和优化抽水井的布设及抽水方案，提高地下水抽出效率，降低成本，采用GMS7.1软件，建立了该场地地下水流及溶质运移数值模型，用以模拟不同修复方案下场地污染羽的变化情况，为场地修复方案的设计和优化提供依据。

2 模型建立

2.1 概念模型

构建合理的水文地质概念模型对于数值模拟的准确性十分重要[25]。基于前期钻孔资料及水文地质调查情况，根据场区所处区域地形地貌、水文地质条件和地下水运动特征，以场区为中心向四周扩大圈定了面积约12.48 km2的模拟区，其南北边界基本平行于等水位线，东西边界基本平行于地下水流向(图1)。

模拟区含水层主要为第四系松散岩类孔隙含水层，目标含水层为潜水，埋深8.2～9.1 m，厚约6.5 m(底部为连续的黏土层)，其岩性主要为中细砂、粉砂，渗透性好，故概化为单层均质含水层。补给来源主要为降水入渗和来自山前的侧向补给；排泄主要为侧向排泄，由于埋深8 m以上，不考虑蒸发排泄。地下水径流主要为自西南山前地带向东北部平原区的侧向径流。根据水文地质条件及建模目的，将地下水流模型概化为二维非稳定流。模型上边界为降雨入渗边界，底部为隔水边界。

图1 研究区概况及边界条件示意图

2.2 数学模型

采用GMS7.1对研究区进行水流和溶质运移模拟，选用其中的MODFLOW和MT3DMS模块。

(1)地下水流场数学模型

按照水文地质概念模型，将模拟区地下水概化为空间二维流；地下水流场的输入输出随时间、空间变化，为非稳定流：

(1)

式中：K——含水层渗透系数/(m·d-1)；

H——水位标高/m；

W——源汇项(降雨、蒸发等)/(m3·d-1)；

μ——潜水层给水度；

t——时间/d；

Ω——渗流计算区域；

Kn——边界法线方向的渗透系数/(m·d-1)；

Γ2——流量边界，包括隔水边界(零流量边界)；

n——边界Γ2的外法线方向；

H0(x,y)——已知初始水位分布；

H(x,y,t)——t时刻的水位标高/m。

(2)溶质运移数学模型

考虑对流、弥散、污染物的源汇项及化学反应项，非稳定流溶质运移模拟：

(2)

式中：θ——地层介质的孔隙度；

Ck——k组分的溶解相浓度/(mg·L-1)；

t——时间/d；

xi，j——沿坐标系轴向的距离/m；

Dij——水动力弥散系数张量/(m2·d-1)；

vi——孔隙水平均实际流速/(m·d-1)；

qs——单位体积含水层流量，源为正，汇为负/d-1；

∑Rn——化学反应项/(mg·(L·d)-1)；

C0(x,y,0)——已知的初始浓度分布；

Ω——代表整个模型区域；

Γ1——表示定浓度边界；

C(x,y,t)——沿Γ1所给的浓度/(mg·L-1)。

在MT3DMS运移模型中采用已知浓度边界(Dirichlet条件)作为进入模型区域的溶质的源。

3 数值模型的识别与验证

利用GMS软件建立地下水流数值模型。模拟区面积约12.48 km2，剖分为250×250个网格，并且在拟修复区内进行网格加密。含水层为上层潜水，地下水运动以水平为主，因此模型概化为二维非稳定流。

根据水文地质条件概化及模拟区四周边界的水文地质特征，将四周边界均设为通用水头边界。其中南北边界近似平行于等水位线，且不远处即为山前地段，故根据当地水力梯度将边界处的水头外推得到1 km外的水头设为边界；东西边界由于垂直于等水位线，认为1 km外的水头与模拟区边界处水头相等。各边界上的水力传导系数为：

C=KB/L

式中：K——渗透系数/(m·d-1)；

B——含水层厚度/m；

L——距离，取1 000 m。

研究区源汇项主要为降雨入渗和侧向补给/排泄。根据场区内抽水试验结果及以往经验值，求取初始的水文地质参数并进行模拟。初始流场参考2016年4月水位统测结果，根据该地区水力梯度对其余地区进行外推概化，得到潜水含水层的初始流场。氨氮初始浓度参照2016年3月地下水水质监测结果，对其余区域进行插值，得到潜水含水层中氨氮污染羽分布。

模型模拟期为1975年1月1日至2035年1月1日，共107个应力期，1975年1月至2017年6月期间每年1个应力期，每个应力期2个时间步长；自2017年7月开始，每个月为1个应力期，每个应力期5个时间步长，以期详细刻画流场变化对溶质运移的影响。

由于研究区流场较为稳定，随时间变化不大，首先建立稳定流模型，并利用地下水位统测资料对模型进行识别和调试，将稳定流模型改为非稳定流模型。本项研究较为关注污染羽的变化，利用场区内监测井的多期动态浓度数据对模拟污染羽进行验证与拟合。

图2 2016年4月地下水流场拟合效果图

用稳定流模拟2016年4月的流场，主要识别大气降水入渗系数、边界侧向量、渗透系数等水文地质因素，经过反复调试，地下水流场拟合效果较好(图2)。可见，建立的地下水流动模型已基本达到模型精度要求，符合水文地质条件，基本反映了目前地下水流场的水力特征，可预测地下水流场及污染物浓度。最终识别的水文地质参数为：渗透系数7.18 m/d，给水度0.15，孔隙度0.3，降雨入渗系数0.14。

在水流模型基础上，利用地下水溶质迁移模拟模块进行污染羽拟合验证。

污染物迁移预测从2013年4月1日起至2042年4月30日，共30年。以1个月为一个应力期，总计360个应力期。模拟期间(2013-04-01—2042-04-30)的降雨量根据研究区1993—2003年的月平均降水量数据推算。场区有两个主要污染源：A区和B区，经查阅该化肥厂历史资料及咨询原厂工人，估计其污染时间分别为：1992—2010年和1975—1991年(B区为老旧车间，后淘汰)，污染浓度约为4 000 mg/L。场地内各时期地下水中氨氮计算浓度与实测浓度拟合效果见图3。

图3中监测井旁标注的实测污染物浓度与模拟污染羽对比可以看出，污染物的模拟浓度与实测浓度拟合较好，说明该溶质运移模型可以较好地展现含水层中污染物的迁移转化。最终识别的溶质运移参数为：分子扩散系数0.001 m2/d，纵向弥散度17 m，一阶吸附常数2×10-10m3/mg，速率常数0.000 3 d-1。

4 抽水方案设计与比选

4.1 抽水方案设计原则与方案概况

由于污染羽的不规则形态和流场变化等复杂因素，需布置不规则分布抽出井群方能达到最佳修复效果。因此需统筹考虑污染羽形状及其所在场地的水文地质条件，通过数值模拟对抽出井群进行优化设计，达到经济高效捕获污染羽的目的。本项目中抽水方案设计原则：合理布设，既要有效控制污染羽扩散，又要将其快速消除；优先在高浓度污染区布井；尽量连续稳定抽水(考虑污水处理厂处理污水的连续性及抽水井利用率)，不疏干即可。氨氮的修复目标值为64.9 mg/L，因此模拟中场区内浓度均小于该目标值且在一段时间内不反弹，视为抽水修复已完成。

参考前人经验初步估计井的数量和位置[26]。首先，通过数值模拟计算单井捕获带和最大抽水量[27]。在模型中进行单井抽水模拟，得出单井捕获带宽度，结合污染羽大小估算所需抽水井数量约为5～7口。其次，布设不同数量及位置的抽水井，通过模拟试算，获得不同布井方案的修复效果，通过对比修复效果确定最佳抽水方案(包括抽水井分布和每口井的抽水量变化)。

常用的布井方法包括直线布井法(中轴线布井法)[21]、三角形布井法、正方形布井法、规则多边形布井等[15]。本研究中，污染羽形状近似略宽的椭圆形，较为适合中轴线布井法与三角形布井法相结合，即在中轴线上布设3～4口抽水井，确保高浓度区中心1口井，上游污染源区附近1～2口井。另外，在中轴线两侧各布设1～2口井，与中轴线上的抽水井呈三角形分布，以确保将所有污染水体抽出。在此基础上，结合大量试算得出3种较好的抽水方案。其中方案一5口井(图4a)，其中3口位于污染羽中轴线(含高浓度区中心1口，源区、下游各1口)，中轴线两侧各1口井；方案二7口井(图5a)，其中3口位于污染羽中轴线(含高浓度区中心1口，源区2口)，下游高浓度区中轴线两侧各2口井；方案三6口井(图6a)，其中4口位于污染羽中轴线(含高浓度区中心1口，源区2口，下游1口)，中轴线两侧各1口井。进行抽出-处理修复的起始时间设为2017年8月1日。

4.2 抽水方案模拟与比选

(1)方案一抽水方案模拟及结果

方案一抽水井共5口，其中井1，3，5位于中轴线上，井2、4位于两侧。抽水方案见表1。在抽水第一阶段，5口井均抽水，使得污染羽不再扩散，进而缩小。第二阶段，井1已不在污染羽范围内，故停止抽水，剩余4口井进行抽水，污染羽进一步缩小。第三阶段，井4不在污染羽范围内，停止抽水，剩余3口井加大抽水量，直到污染物浓度均小于64.9 mg/L，完成修复。

模拟结果显示，地下水经过抽出处理在修复第47个月水质达标(图4d)，修复周期为46个月7天，总抽水量为48.97×104m3。可以看出，由于井数略少，导致捕获污染羽的速度较慢。尤其是第24个月起(图4c)，虽然残余污染羽面积不大，但由于该阶段抽水井并未位于污染羽中心，而在其周围，污染物抽出效率变低，无法很快去除局部小范围污染羽残留。

图4 方案一修复过程中污染羽的变化

表1 方案一抽水方案

(2)方案二模拟结果

考虑到方案一井数不足导致修复期太长，方案二将抽水井增加为7口，中轴线上井数仍为3口，井7位于高浓度区，井5和6分别位于源区B区和A区，另外，在中轴线两侧各布设两口抽水井(井1～4)。采用先四周后中间的抽水方式。第一阶段，外围五口井(井1～5)进行连续抽水(60～80 m3/d)；第二阶段，污染羽缩小(图5b)，井1～5均已不在污染羽范围内，为提高效率，只对中间两口抽水井井6，7以不致疏干的最大流速抽水(80 m3/d)。方案二抽水方案见表2。模拟结果显示，地下水经修复在第25个月水质达标(图5d)，修复周期为24个月10天，比方案一缩短22个月，总抽水量为21.22×104m3。

图5 方案二修复过程中污染羽的变化

表2 方案二抽水方案

(3)方案三模拟结果

方案三设计了6口抽水井，其中4口位于污染羽中轴线上：井2位于高浓度区，井4、6分别位于源区A区和B区，井1在下游，确保污染羽不向下游扩散；井3和井5分布在中轴线两侧。抽水共分为四个阶段，根据污染羽的缩小情况及时调整抽水方案，保证在污染羽中心一直有较大流量的抽水井进行抽水，力图快速降低污染物浓度。抽水方案见表3。第一阶段，首先以位于高浓度区中的井2为主要抽水井，进行最大流量的抽水(75 m3/d)，周围井1，3，5，6用较小流量抽水，避免污染羽向四周扩张；第二阶段，高浓度区浓度显著降低后，井2和井4位于相对较高浓度区中(图6b)，对其进行较大流量抽水(60 m3/d和75 m3/d)，其余4口井进行小流量抽水，为避免疏干，距离井2，4距离较近的井3和5抽水流量最小；第三阶段，井4位于此阶段的污染羽中心(图6c)，则对其进行较大流量抽水(75 m3/d)，井2，3，5进行小流量抽水(35 m3/d)，而井1和6停止抽水；在第四阶段，有小范围污染羽残余(图6d)，仅对污染羽内的井4进行大流量抽水(115 m3/d)，以快速完成修复。

模拟结果显示，地下水经过抽出处理修复在第23个月水质达标(图6e)，修复周期仅23个月，总抽水量为17.07×104m3。可以看出，该方案随污染羽形状改变及时调整抽水方案后，修复效率进一步提高，修复周期比方案二进一步缩短41 d，总抽水量也进一步减少。

4.3 抽水方案比选

抽水方案的比选一般考虑以下重要因素：抽水井数、总抽水量及修复完成时间[17]。因为抽出污水最终都要经过处理，井数和总抽水量决定了修复工程的经济成本，而修复完成时间也影响了整个工程的经济与时间成本。三个抽水方案对比结果见表3。

图6 方案三修复过程中污染羽的变化

表3 方案三抽水方案

表4 抽水方案对比

方案三修复时间最短，且井数适中，总抽水量最小，污水处理的负担较小，为最优方案。根据优化过程可以看出，当井数过少时，会存在污染羽死角，延缓修复进程；当井数过多时，井数较密部位的含水层很容易疏干，各井抽水流速均不大，相应的，修复时间也很难进一步缩短。方案三之所以效果最好，一方面是由于污染羽中轴线上的井数量最多(4口)，保证了初始污染羽高浓度区和污染源区均有抽水井将高浓度污染羽尽快抽出，另一方面是由于该方案根据模拟过程中污染羽的形状变化及高浓度区的位置变化，随时调整抽水方案，在23个月的修复周期内细分了多达四个阶段，对各阶段的污染羽中心进行大流量抽水。该优化思路对于地下水污染场地的修复方案设计有一定的参考意义，即在场地修复方案设计中，密切关注污染羽的变化情况并采用溶质运移模拟的手段进行预测，可以为下一步抽水方案的调整提供依据，有助于提高修复效率、缩短修复周期、节约修复成本。

5 结论

(1)地下水流和溶质运移模拟是地下水污染羽的抽出处理效果预测的重要手段，也可为抽水方案的优化提供重要帮助。

(2)抽水井布设应优先在污染羽中轴线上布井，保证高浓度区和污染源区均设有抽水井；另外结合三角形布井法，在中轴线两侧酌情布井。

(3)抽水过程中，抽水方案应随着污染羽外围形状及其高浓度区的变化及时调整，调整的核心原则是保证污染羽中心有抽水井进行较大流量抽水。

(4)由于污染羽在含水层中的变化受多种因素影响，在修复工程实施过程中，应持续监测流场及污染羽变化，并据此对模型进行调整，随时调整抽水方案，以提高修复效率。