某铁合金厂地下水铬污染运移三维动力学过程数值模拟分析

Numerical Simulation Analysis of Three-Dimensional Kinetic Process of Chromium Pollution Transport in Groundwater of a Ferroalloy Plant

摘要: 地下水重金属污染运移规律对污染治理十分重要，重金属污染数值模拟是了解地下水重金属污染运移过程的重要手段，污染物运移模拟结果贴合实际变化规律很大程度上依赖于地下水重金属动力学运移模型、参数和算法的设置，国内尚缺乏可根据场地污染情况自定义数学模型与算法的地下水模拟软件。本文以铬污染为例，基于现有的勘察工程数据和实验得到的污染场地模型参数，构建铬污染运移三维动力学数学模型，自编模块实现有限差分离散化方法构建并求解三维数值模型，模拟六价铬在污染场地中的三维运移过程，为分析铬污染物运移规律提供技术支撑。

Abstract: The transport law of heavy metal pollution in groundwater is very important for pollution control. Numerical simulation of heavy metal pollution is an important means to understand the transport process of heavy metal pollution in groundwater. The simulation results of pollutant transport are consistent with the actual change law, which largely depends on the setting of dynamic transport model, parameters and algorithm of heavy metals in groundwater. There is still a lack of groundwater simulation software which can customize the mathematical model and algorithm according to the site pollution in China. We take chromium pollution as an example, construct a three-di- mensional kinetic mathematical model of chromium pollution transport based on the existing survey engineering data and experimental model parameters of the contaminated site, and construct and solve a three-dimensional numerical model based on a self-programmed module to implement a finite difference discretization method to simulate the three-dimensional transport process of hexavalent chromium in the contaminated site, providing technical support for the analysis of chromium pollutant transport law.

1. 引言

地下水是重要的饮用水水源和战略资源，地下水的污染防治也一直是地下水管理的重要环节 [1] 。而近年来，随着工业生产过程中废渣、废水随意堆放 [2] [3] [4] ，导致地下水重金属污染含量超标严重影响下游河流水质，并对人体健康造成较大的隐患 [5] [6] 。重金属污染数值模拟，是以某一时刻的地下水重金属污染信息为起点，根据污染物运移规律模型和数值离散化方法构建数值模型并量化模型参数 [7] [8] [9] [10] ，从而监测一定时间段内的污染物浓度时空分布。虽然当前Visual Modflow、GMS等国外的三维数值模拟软件十分成熟 [11] [12] [13] [14] ，但针对重金属污染实际模拟问题，国内尚缺乏可根据场地污染情况自定义数学模型与算法的自主知识产权软件。

研究区是一处典型的重金属铬污染区域，大量的铬渣未经处理直接堆积。在降雨及渗流作用下，渗入地下水中，该过程三价铬在地下水中转化为毒性极强的六价铬，造成周围区域的土壤地下水严重污染 [15] [16] ，因此了解地下水铬污染空间分布、模拟污染运移过程对于污染场地的修复与治理具有重大意义。本文以重金属铬为例，收集与整理污染场地降雨、流速等各类水文地质信息和污染源分布，以地下水渗流模型和溶质运移模型为基础，实现三维有限差分数值模型离散化求解方案，模拟重金属铬在地下水中的三维动力学运移过程，为了解地下水重金属污染过程、预测未来污染运移趋势提供技术支撑。

2. 研究区水文地质概况

研究区属中亚热带季风湿润气候区域，年平均气温17.1℃，日照时数1558.7小时，降雨量1326.8毫米，蒸发量1284.1毫米。属华南湘赣丘陵区、加里东褶皱带，其周围为红岩及第四纪堆积物，是中生代断裂凹陷的产物。古河流冲积物散布在二、三级阶地，上为亚粘土，下部砂砾层。近代冲积物为灰色、灰黄色亚砂土，厚2~5 m，呈微酸性。地下水主要来源于上游的湘江支流补给和大气降水，并以孔隙水和岩溶水状态，由西北往东南方向的涟水河排泄。

研究区地层从上到下分为四层(图1)：1) 杂填土：具有上层滞水性质，并接受大气降雨、地表径流与厂房排水补给，灰黑色为主，由炉渣、建筑垃圾等构成，厚度3~4米；2) 粉质粘土：黄褐色，成分以粘土为主，含砂10%，韧性中等，厚度2~4米，弱–无透水层，部分含有上层滞水；3) 圆砾：中等孔隙水，接受大气降水和人工排水的越流补给以及地下水的侧向径流补给，黄褐色，含砾50%左右，粒径以20~40 mm为主，厚度3~4米；4) 风化泥质粉砂岩：灰色为主，岩石风化强烈，厚度达9~10 m，弱透水，含少量裂隙水裂隙面，泥质充填。

Figure 1. Stratigraphic model of the study area

图1. 研究区地层模型

3. 研究方法

本文进行了工程地质数据处理，从工程地质剖面数据中提取地层及水位信息并构建场地三维模型，基于地下水渗流方程和溶质运移方程构建铬污染三维运移过程耦合数学模型，为高效地获取污染模拟结果，采用三维有限差分方法离散化耦合模型 [17] [18] [19] ，从而得到三维有限差分数值模型，采用共轭梯度法和稳定双共轭梯度法求解数值模型形成的稀疏线性系统，将以上数值模拟过程使用C++语言编写为重金属污染数值模拟模块，并验证和测试程序地正确性和稳定性，最后在分析影响模拟效果的关键因素后，模拟计算污染三维运移过程，结合场地污染情况验证并分析污染物运移规律。

3.1. 数学模型

目前，在重金属污染数值模拟研究中，常采用对流–弥散迁移理论 [20] ，考虑土壤地下水性质和对流、弥散、流体源/汇、吸附等作用，以地下水渗流方程和溶质运移方程为核心构建铬污染三维运移过程耦合模型 [21] [22] [23] [24] 。

$S\frac{h^t-h^{t-1}}{\Delta t}={\displaystyle \underset{i=1}{\overset{D}{\sum }}\frac{\partial }{\partial x_i}\left(k_i\frac{\partial h^t}{\partial x_i}\right)+q_s}$ (1)

$\theta R\frac{\partial C}{\partial t}={\displaystyle \underset{i=1}{\overset{D}{\sum }}\frac{\partial }{\partial x_i}\left(\theta {\displaystyle \underset{j=1}{\overset{D}{\sum }}{D}_{ij}\frac{\partial C^t}{\partial x_j}}\right)}+{\displaystyle \underset{i=1}{\overset{D}{\sum }}\frac{\partial }{\partial x_i}\left(q_i{C}^t\right)+{\lambda }_1\theta C^t}$ (2)

式中，t为时间(T)；S为承压含水层单位储水量(1/m)；q_s为流体源/汇项(mg/(m³*s))；F为含水层厚度(m)；K_xx/yy/zz为x/y/z方向的渗透系数(m/s)；D_xx/yy/zz为x/y/z方向的弥散系数(m²/s)；q_x/y/z为x/y/z方向流速(m/s)，θ为地层介质有效孔隙度，无量纲；R为延迟因子，无量纲；ρ_b为沉积物干密度(mg/m³)，Cs为浓度源汇项(mg/s⁻¹*m⁻³)，λ₁为吸附化学反应的速率系数，无量纲；C为溶解项浓度(mg)。

3.2. 数值模型

对数学模型数值离散化，得到三维有限差分渗流场和浓度场数值模型，联立每个有效节点的线性代数方程，即得到每个模拟时间步下的线性方程组：

${\displaystyle \underset{k=1}{\overset{N}{\sum }}{h}_k^{\left(t\right)}\left(A_{hk}^{\left(h\right)}+B_{hk}^{\left(h\right)}\right)}=C_h^{\left(h\right)}-{\displaystyle \underset{k=1}{\overset{N}{\sum }}{h}_{hk}^{\left(t-1\right)}{B}_{hk}^{\left(h\right)}},l=1,\cdots ,N$ (3)

$A_{hk}^{\left(h\right)}=-{\displaystyle \underset{i=1}{\overset{D}{\sum }}{k}_i\frac{\partial H_k}{\partial x_i}},B_{hk}^{\left(h\right)}=-\frac{S_s}{\Delta t}{H}_k\left(x\right),C_l^{\left(h\right)}=q_s$ (4)

${\displaystyle \underset{k=1}{\overset{N}{\sum }}{C}_k\left(A_{ck}^{\left(C\right)}+B_{ck}^{\left(C\right)}+C_{ck}^{\left(C\right)}\right)}=D_{ck}^{\left(C\right)}+{\displaystyle \underset{k=1}{\overset{N}{\sum }}{C}_k^{\left(t-1\right)}{A}_{ck}^{\left(C\right)}},l=1,\cdots ,N$ (5)

$A_{ck}^{\left(C\right)}=\theta \left(\frac{R}{\Delta t}+{\lambda }_1\right)H_k\left(x\right),B_{ck}^{\left(C\right)}=\theta {\displaystyle \underset{i=1}{\overset{D}{\sum }}{\displaystyle \underset{i=1}{\overset{D}{\sum }}{D}_{ii}\frac{\partial H_k}{\partial x_i}}},C_{ck}^{\left(C\right)}={\displaystyle \underset{i=1}{\overset{D}{\sum }}{q}_i\frac{\partial H_k}{\partial x_i}},D_{ck}^{\left(C\right)}=M.$ (6)

式中，N表示有限差分格网节点的个数，h^t表示第t + 1时刻有限差分格网中心点的水头值，h^t−1表示第t − 1时刻有限差分格网中心点的水头值，C^t表示第t时刻有限差分格网中心点处的浓度值，C^t−1表示第t − 1时刻有限差分格网中心点处的浓度值。

3.3. 数值求解技术

针对重金属污染三维运移模型产生的稀疏线性方程组求解问题，数值求解技术采用共轭梯度法(求解渗流方程)和双稳定共轭梯度法(求解溶质运移方程)。

共轭梯度法是一种路径优化的梯度迭代下降方法，将最优解x*和初始值x₀的差表示为共轭向量的线性组合(式(1-28))，共轭梯度法在搜索共轭向量p_n时只需要利用上一个共轭向量p_n−1，每一次迭代的新搜索方向负梯度方向和上一个搜索方向的线性组合(式(1-29)、(1-30))。

$x^{\ast }-x_0={\alpha }_0{p}_0+{\alpha }_1{p}_1+\dots +{\alpha }_{n-1}{p}_{n-1}$ (7)

$p_k=-r_k+{\beta }_k{p}_{k-1}$ (8)

其中α代表搜索步长，p代表搜索方向。

${\beta }_k=\frac{r_k^T{r}_k}{r_{k-1}^T{r}_{k-1}}$ (9)

稳定双共轭梯度法作为一般线性方程组的求解方法，对线性方程组求解问题中矩阵的正定对称性质无要求，是对共轭梯度法的优化，具备收敛速度快、稳定性高、不需要任何外来参数等的特点。

3.4. 模拟效果的关键因素分析

重金属污染数值模拟效果受污染物运移数学模型的模型参数、数值模型的离散化精度和数值求解技术类型及收敛条件等因素的影响。

污染物运移数学模型的模型参数包括渗透系数K、弥散系数D、有效孔隙度θ、储水率S (一般默认为0.0001)、延迟因子R等，根据土壤性质可以划定模型参数的大致范围，并现场采样得到准确的模型参数值(表1)。

针对数值模拟结果细部的精度要求，有限差分数值模型的网格离散化数量越多，精度越高，但计算耗时更长，对数值求解技术的性能要求更高。重金属污染场地内一般将网格边长设置为2~10 m，本文的有限差分三维网格边长dx、dy、dz分别为10 m、10 m、0.5 m。

不同的数值求解技术求解数值模型产生的不同线性方程组，方程组的对称正定性质影响数值求解技术的选择，如将共轭梯度法应用于求解非对称正定矩阵，则会导致结果错误。同时，数值求解技术的绝对误差和相对误差收敛条件影响共轭梯度法和稳定双共轭梯度法的迭代次数和求解效率，一般将绝对误差设置为1e−3，相对误差设置为1e−5。

Table 1. Model parameter ranges in soil groundwater

表1. 土壤地下水中模型参数范围表

4. 三维场地污染羽分析

4.1. 模型初始条件和边界条件设置

本研究采集场地土壤地下水中用于污染运移模拟的模型参数，包括渗透系数、有效孔隙度、弥散系数、储水率等，模拟主要地层为粉质粘土层(上层滞水)和圆砾层(含水层)，其中泥岩层近似为不透水层，如表2所示。

Table 2. Main model parameters of chromium pollution transport process

表2. 铬污染运移过程主要模型参数表

场地实测的土壤地下水性质参数和场地三维模型决定初始渗流场分布(图2)，水头补给主要为降雨形成的上层滞水，根据涟水水位设定两个定水头边界，场地西北侧定水头为54~56 m，东南侧定水头边界46~48 m，水流主方向为西北至东南(图3)。

重金属的初始污染区域位于研究区中部，在土壤中的污染浓度峰值超过3000 mg/L，且主要集中在杂填土层的下部和粉质粘土层的上部，随着降雨下渗作用，少部分污染物溶质向含水层移动，粉质粘土层中的大量六价铬在长期的降雨淋滤作用下向圆砾层渗透(图4)，在圆砾层上表面形成源源不断地污染补给(图5)，六价铬溶质边界设定为自由流出边界。

4.2. 铬污染运移三维动力学模拟

为验证自编实现的地下水铬污染运移有限差分算法的正确性，根据实际场地三维模型、初始条件、边界条件和土壤地下水性质参数，采用有限差分法模拟铬污染场地(仅模拟含水层)的三维时空变化过程。模拟时间从2018年3月15日到2020年3月15日，单位模拟步长为1天，其中在第366天模拟挖除杂填土层和粉质粘土层中含有六价铬的土壤，观察含水层中铬污染物的变化。

Figure 2. Three-dimensional distribution schematic diagram of initial water

图2. 初始渗流场三维分布图

Figure 3. Schematic diagram of flow boundary conditions

图3. 水流边界条件示意图

Figure 4. Schematic diagram of hexavalent chromium contamination in the lower part of the powdered clay layer

图4. 粉质粘土层下部六价铬污染示意图

Figure 5. Three-dimensional distribution schematic diagram of initial concentration

图5. 初始浓度场三维分布图

4.2.1. 六价铬三维动力学运移模拟结果

渗流场模拟结果显示，初始时刻在场地中部污染汇入处设置降雨补给(图6)，补给量根据降雨量1326.8 mm/a和蒸发量1284.1 mm/a计算，水头补给带动水流运动加快，并产生向四周以及下方扩散的趋势，随着迭代步数增加，三维渗流场内水头可以保持稳定，水流平缓均匀流动，呈现西北向东南方向的地势更低处运移，水头值范围在40~60 m之间。三维渗流场模拟的是含水层的水流运动，因此在垂直方向上无明显差异。

Figure 6. Schematic diagram of three-dimensional seepage field rainfall recharge location

图6. 三维渗流场降雨补给位置

4.2.2. 三维浓度场动力学模拟结果

六价铬污染模拟结果显示(图7)，模拟结果每30天记录一次，在降雨作用下污染物的不断下渗，粉质粘土层下部与圆砾层上部交界处堆积量逐渐增加，圆砾层中的污染物在对流弥散作用下随水流方向往东南角的涟水河迁移，其中，对流起主导作用。在前365天中污染峰值从1.63 mg/L到2.45 mg/L逐渐增大，且污染中心无移动，下游污染量逐渐增加，随着带有铬污染的杂填土和粉质粘土层挖除后，第420天时污染中心向东南侧场地边界运移100 m，此时污染羽峰值为1.46 mg/L，第480天时污染中心运移至距离东南侧边界120 m附近，污染羽峰值为1.17 mg/L。在第540天时，污染中心迁移到场地边界，污染羽峰值浓度为1.00 mg/L。第720天时仅在场地边界处存在小范围浓度超标区域，浓度羽峰值为0.30 mg/L。土壤中六价铬在地表附近降水下渗作用主导，运移以垂向为主，水平运移速率较小；到达含水层后，水平运移速率明显增加，含水层中的污染物浓度先是逐渐增大，当地表附近的污染物停止向下运移且已存在的大部分污染物全部进入含水层后，在对流弥散作用下污染物浓度开始逐渐减小并逐渐向场地东南角下游边界迁移，整体上符合场地中六价铬的实际运移规律。

Figure 7. Three-dimensional dynamic simulation process of chromium pollution

图7. 铬污染三维动力学模拟过程

5. 结论

本文以某铁合金厂铬污染为例，开发一套可自定义数学模型、参数及算法的、具备国产自主知识产权的地下水重金属污染数值模拟软件，数值模拟流程包括数学模型构建、数值模型构建与求解和污染模拟结果分析，本文主要研究成果如下：

1)本文针对重金属污染三维动力学运移模型，自编实现有限差分数值模拟算法，为地下水重金属污染运移数值模拟提供技术支撑。

2)结合工程勘测数据、污染场地三维信息和有限差分数值求解算法模拟污染场地六价铬的运移规律，场地水流总方向为由西北至东南，六价铬从粉质粘土层随着降雨下渗到圆砾层上表面，在圆砾层中对流作用下往东南方向场地边界外移动，符合六价铬在场地内的三维运移规律。

3)本文所使用的数值模拟算法支持模块扩展，因此后续可以加入有限元等数值模拟算法以及多重网格等数值求解技术。

基金项目

国家重点研发计划课题项目(2019YFC1805905)。