基于COMSOL Multiphysics 软件的地下水污染物运移模拟研究

基于COMSOL Multiphysics 软件的地下水污染物运移模拟研究

张 洋

(天津华北地质勘查局核工业二四七大队,天津 301800)

[摘 要]  地下水环境质量是衡量生态文明建设水平的重要指标,应用科学有效的计算方法模拟地下水污染物渗漏的影响范围及程度,对于地下水污染的防治有着极其重要的意义。COMSOL Multiphysics作为一款多物理场建模与仿真软件,能够较好地结合现实水文地质条件,有效地模拟自然界一些水文地质过程发生和发展过程,对进行以预防和修复治理为目的的地下水污染物数值模拟预测具有十分重要的应用价值。本文根据研究区的水文地质特征,利用COMSOL Multiphysics 软件构建了地下水数值模型,对研究区地下水污染物运移规律进行分析,对污染物运移轨迹进行模拟与预测。预测结果表明,氯仿、石油烃污染物在释放1 000 a时的影响范围分别是6.3×10 6  m 2 和5.7×10 6  m 2 。由于地下水流速相对较慢,污染物扩散速度极其缓慢,研究区内污染物需要在百年时间量级上才可能形成较大规模扩散,短期内对地下水环境的污染风险较小。模拟结果为研究区的地下水环境评价和防范措施提供了科学依据。

[关键词]  COMSOL Multiphysics;地下水;污染物运移;模拟预测

0 引言

地下水既是不可替代的水资源,又是重要的环境要素和生态支撑条件。地下水污染不仅会对生态环境造成严重影响,还会制约工农业发展,对广大人民群众生命安全造成极大危害。地下水污染防治是贯彻习近平生态文明思想,落实《中华人民共和国水污染防治法》《水污染防治行动计划》的重要举措。

近年来,地下水环境质量越来越受到重视,对一些可能会造成地下水污染的企业,应用科学有效的计算方法模拟污染物泄露的影响范围及程度,对于地下水污染的防治有着极其重要的意义。随着地下水科学与计算机技术的快速发展,地下水进行水流和水质运移采用数值模拟进行预测研究也得到了广泛的应用。目前地下水系统数值模拟方法主要有有限单元法 (FEM)、有限差分法 (FDM)、有限分析法 (FAM)和边界元法 (BEM)等 [1] 。此次利用 COMSOL Multiphysics 软件来模拟研究区地下水污染物运移情况,为研究区地下水污染物的防治提供科学依据。

1 COMSOL Multiphysics 软件

COMSOL Multiphysics是瑞典康模数尔公司于1998年发布的一款多物理场建模与仿真软件,涵盖力学、电磁场、流体、传热、化工、MEMS、声学等专业领域。软件使用全局隐式方法同时求解Richards方程的三维形式,可同时求解包气带-饱和带渗流、甚至裂隙流过程。同时,由于其本身基于多场耦合特性,求解所得到的渗流场可直接与三维反应性物质运移模拟(对流弥散方程)耦合,可考虑固、液、气相中自由、饱和、变饱和流动对应的溶质运移,最终求得地下污染物在非饱和带-饱和带中的迁移过程,真实意义上实现多组分、多过程耦合的数值模拟。同时,COMSOL地下水模块可考虑多孔介质中的热传导和对流、相变、线弹性等过程,同时可与任意物理场(化学反应、动力学、电磁场)耦合 [2] 。

2 研究区水文地质条件

2.1 地层岩性

研究区为天津市某工业园区,地貌特征属典型冲海积平原型,地势比较平坦。研究区埋深约30 m深度范围内,地基土按成因年代可分为以下6层,按力学性质可进一步划分为13个亚层 [3] 。自上而下依次为人工填土层(Qml)2个亚层、新近冲积层(Q 4 3N al)、全新统中组海相沉积层(Q 4 2 m)6个亚层、全新统下组沼泽相沉积层(Q 4 1 h)、全新统下组陆相冲积层(Q 4 1 al)2个亚层、上更新统第五组陆相冲积层(Q 3 e al)。

2.2 地下水赋存条件

包气带:主要指地下水位以上的人工填土层(Qml)杂填土(地层编号① 1 )、素填土(地层编号① 2 )及新近冲积层(Q 4 3N al)粉质黏土(地层编号 ③ 1 ),包气带厚度一般约为0.50～1.70 m。

潜水含水层:主要由地下水位以下的人工填土层(Qml)素填土(地层编号① 2 )、新近冲积层(Q 4 3N al)粉质黏土(地层编号③ 1 )、全新统中组海相沉积层(Q 4 2 m)淤泥质黏土(地层编号⑥ 1-1 )、粉质黏土(地层编号⑥ 1-2 )、粉土(地层编号⑥ 1-3 )、粉质黏土(地层编号⑥ 1-4) 、淤泥质黏土(地层编号⑥ 2 )、粉质黏土(地层编号⑥4)组成,底板埋深一般为15.20～18.00 m,厚度约为13.90～17.10 m。

潜水相对隔水层:主要由全新统下组沼泽相沉积层(Q 4 1 h)粉质黏土(地层编号⑦)、全新统下组陆相冲积层(Q 4 1 al)粉质黏土(地层编号⑧ 1 )组成,该层总体透水性以极微透水为主,具相对隔水作用。

2.3 地下水补径排条件

天津平原区潜水在自然条件下侧向径流极为缓慢,局部受地形、地势、地表水影响;垂向上主要由大气降水补给、以蒸发形式排泄,体现为入渗-蒸发动态类型。

本次量测稳定自然水位(2022年6-7月),研究区潜水水位埋深一般介于0.72～2.68 m,水位高程一般介于1.27～2.34 m。目前研究区潜水含水层形成了由西南向东北的地下水流场,潜水平均水力坡度一般约为0.33‰。

3 地下水污染物运移模型建立

3.1 水文地质概念模型

依据沉积物和含水层特性,将研究区30 m以浅的地层属性从上至下概化为6层,岩层分层概化如图1所示。6个岩层依次包括1-杂填土、2-粉质黏土、3-粉粘/粉土、4-粉/泥质黏土、5-粉土、6-粘土/粉粘,每层的平均厚度分别约为1.4 m、3.8 m、5.0 m、9.7 m、4.9 m、5.0 m。其中,1-杂填土、2-粉质粘土、3-粉粘/粉土概化为潜水含水层,5-粉土概化为含水层,其余层位(4和6)为弱透水或者隔水层。

图1 研究区岩性分层概化示意图

3.2 地下水流动数值模型

3.2.1 地下水渗流数学控制方程

通过求解Richards方程用来模拟地下水在非饱和带(包气带)中的流动过程 [4] :

(1)

(2)

式中:p是压力,ρ是水的密度,g是重力加速度,C m 是比水容重,S e 是有效饱和度,S是储水系数,k s 是饱和渗透率,k r 是相对渗透率,μ是水的粘滞系数,z是高程,Q m 代表源汇项,U是速度向量。其中C m ,S e 和kr可根据van Genuchten相关模型确定:

(3)

(4)

(5)

(6)

式中:H p =p/(ρg),α是进气压力值,n、m (m=1-1/n) 和i是拟合参数,θ s 是饱和含水量,θ r 是残余含水量。

在饱和带中,可以通过求解三维地下水控制方程获取渗流过程:

(7)

式中:K为渗透系数(下标x、y和z代表流动方向);μ s 为储水率;H为水位;t为时间;W为源汇项。

3.2.2 渗流数值模型

模型西南部与河流相接触,地表水设置为给定水头边界(水位可根据时间变化);模型其余侧面由于流量和具体水文地质特性未知,根据实测水位数据将其设置为三类边界,三类边界条件参数(导水度)根据实测水位进行模型校正;模型顶部接受降雨入渗补给(平均降雨量为1.56 mm/d),考虑入渗补给系数大致为0.2～0.4;工业园区降雨入渗系数假设可以忽略,因此设置为零通量边界;模型底部设置为隔水边界。具体边界设置如图2所示。

模型在水平和垂直方向上尺度差异较大,为此采用不规则三角网格剖分格式。三角网格可以更好的捕捉到不规则的边界形状,而且能够有效表征任意不规则研究区,从而保障数值结果的准确性。模型网格剖分情况如图3所示。

图2 水流模型边界条件示意图 图3 数值模拟三维网格示意图

3.2.3 模型验证

在建立溶质模型结果之前,需要对模型结果的可靠性进行检验。经过不断调整优化水文地质参数(表1),获得较为准确模拟流场。通过提取研究区地下水位(潜水含水层)的模拟值,将其与实测数据进行对比验证(图4)。结果表明:模型模拟流场基本与测量水位吻合,一定程度上验证了数值模型的可靠性。

表1 研究区水文地质参数一览表

图4 实测水位与模拟水位对比图

3.3 地下水溶质运移模型

考虑工业园区化工厂的影响,根据实际情况需求,选取氯仿、石油烃等2种污染物作为溶质运移模拟的污染组分。研究区内污染源主要在化工厂指定场地内,污染浓度参考国家地下水质量标准(IV类)或污染物溶解度来确定。本次地下水溶质运移模拟过程主要考虑对流、弥散等物理过程以及吸附等化学过程 [5] 。

利用对流-弥散-化学反应方程来刻画地下水系统中多组分反应性物质的迁移转化过程:

(8)

式中:T为滞后系数;θ s 为介质孔隙度;C j 为物质j的浓度;D为水动力弥散张量;q为水流通量;R j 为物质j的源汇项(表征化学反应速率)。

本次研究是建立在忽略污染物自身及其与其他物质发生化学反应的基础上,所以模型需要给出的参数主要有含水介质有效孔隙度、水动力弥散系数等。由于没有现场实测数据,本次模型的有效孔隙度和弥散度采用前人总结的经验值。研究区含水层岩性主要为孔隙含水岩层,本次模拟根据不同岩层,赋值有效孔隙度;其它参数包括:纵向弥散度为0.1 m,横向弥散度为纵向弥散度的1/3。

4 地下水污染物运移模拟预测

本次模拟氯仿、石油烃污染物在释放1 000 a内的污染情况,污染源浓度参考国家地下水质量标准(IV类)或污染物溶解度来确定,即氯仿0.3 mg/L、石油烃15 mg/L。为了更好表征不同释放源位置以及浓度差异对污染物运移范围的影响,图5提供了模拟结束时污染物浓度在空间上的分布。可以看出,污染源位于模拟区东侧对地下水的影响范围更大,这是因为地下水在该处流速相对较大,且整体从模型东侧流出研究区域边界;而处于模型西部或者中部地区的污染源,由于受到复杂地下水流流场影响,污染范围相对较小。

在污染源持续不断的条件下,各类污染物的影响范围依据污染源释放浓度的1%作为判别标准,得到各类污染物在水平方向上的影响范围。可见,各类污染物影响范围随时间呈非线性增长趋势,即随时间增加,污染范围增长速率越快。最终得到各类污染物的在模拟结束时期的影响范围是:氯仿6.3×10 6  m 2 、石油烃5.7×10 6  m 2 。

图5 污染源释放1 000 a后污染范围示意图

5 结语

COMSOL Multiphysics能够较好地结合现实水文地质条件,有效地模拟自然界一些水文地质过程发生和发展过程,对进行以预防和修复治理为目的的地下水污染物数值模拟预测具有十分重要的意义。总体来看,研究区内污染物需要在百年时间量级上才可能形成较大规模扩散。从实际污染防治角度来看,氯仿、石油烃等污染物在目前污染情境下,短期内对地下水环境的污染风险较小。