铬污染迁移修复过程自适应仿真外挂系统计算与三维可视化数值模拟

Numerical Simulation of Chromium Pollution Migration Remediation Process with Adaptive Simulation Plug-In System Calculation and Three-Dimensional Visualization

摘要: 铬污染迁移修复过程数值模拟是研究铬渣遗留场地铬污染三维时空分布，优化铬污染治理方案的重要手段。由于铬渣遗留场地存在三维结构及参数的时空异质性和修复技术的综合多样性，现有的数值模拟技术无法准确地模拟复杂场地铬污染迁移修复过程，为了解决这一问题，本文提出了铬污染迁移转化与修复过程模拟引擎–自适应外挂耦合计算架构及仿真范式。以某铬渣遗留场地为例，结合场地勘探数据，实验室数据以及监测数据，利用铬污染迁移修复过程模拟引擎–自适应外挂耦合计算架构，进行铬污染迁移修复过程模拟，得到了场地铬污染三维时空分布情况。

Abstract: Numerical simulation of chromium pollution migration remediation process is an important means to study the three-dimensional spatio-temporal distribution of chromium pollution in chromium slag legacy sites and to optimize the chromium pollution management scheme. Due to the spatial and temporal heterogeneity of the three-dimensional structure and parameters of the chromium residue legacy site and the comprehensive diversity of the remediation technology, the existing numerical simulation technology can not accurately simulate the chromium contamination migration and remediation process of complex sites, in order to solve this problem, this paper puts forward the chromium contamination migration transformation and bioremediation process simulation engine-adaptive plug-in coupled computational architecture and simulation paradigm. Taking a chrome residue site as an example, combined with the site exploration data, laboratory data and monitoring data, the chromium pollution migration and remediation process simulation is carried out by utilizing the chromium pollution migration and remediation process simulation engine-adaptive plug-in coupling computational architecture, and the three-dimensional spatial and temporal distribution of chromium pollution at the site is obtained.

1. 引言

铬在自然界中存在广泛，一般以化合物的形式存在，最常见的价态是正三价(Cr(III))和正六价(Cr(VI))，具有高毒性、生物累积性和难去除性[1] [2] [3] [4]。Cr(III)在中性或碱性条件下溶解度较低，过量摄入会对人体健康造成危害；Cr(VI)表现出较强的氧化作用，其毒性是Cr(III)的数倍，具有致突变性、致癌性以及致畸胎性，不仅对人体健康构成严重威胁，导致组织损伤、呼吸道癌症、肝肾功能衰竭等严重问题，还会危及动植物生长，造成生态系统紊乱[5] [6] [7]。随着全球经济的蓬勃发展以及工业化的高速推进，金属加工和合金生产、木材和纸浆加工、石油提炼等多种工业活动逐年增加[8] [9] [10]，工业废水、废气、废渣随意排放，导致土壤和地下水中铬污染问题日益突显，严重威胁着生态环境和人类安全[11] [12] [13] [14]。为了实现可持续发展，研究土壤地下水中铬污染迁移、转化、修复过程，提出最优的铬污染治理方案迫在眉睫。

由于实时监测和采样分析的操作复杂，工作量大，成本较高，数值模拟技术以其全面性，灵活性和可视化能力成为研究土壤地下水中铬污染迁移转化修复过程及时空分布的强有力手段[15] [16] [17]。在过去的十几年中，无数学者致力于对土壤地下水中铬污染迁移修复过程数值模拟的研究[18] [19]，Visual MODFLOW、GMS以及FEFLOW等数值模拟软件应运而生[20] [21] [22]。土壤地下水中铬污染迁移修复过程可能受季节、降雨等自然因素影响，铬渣遗留场地的三维结构及参数存在时空异质性，即铬渣遗留场地的三维结构及参数不仅会随空间的变化而变化，还会随时间的变化而变化。然而，传统的数值模拟软件忽略了铬渣遗留场地的三维结构及参数存在时空异质性，无法实时修正模型参数，从而无法准确有效地模拟复杂场地的铬污染迁移修复过程。同时，传统的数值模拟软件无法对渗透反应墙、原位注入等修复过程进行模拟，也无法实现铬渣遗留场地修复技术综合多样性的过程模拟。

为了解决上述问题，本研究考虑铬渣遗留场地三维结构及属性的时空异质性和修复技术的综合多样性，提出了铬污染迁移修复过程模拟引擎–自适应外挂系统耦合计算架构及仿真范式，以数学模型为基础，数值算法为手段，模拟引擎为核心，外挂仿真为特色，建立铬污染迁移修复模拟软件，实现了土壤–地下水异质系统的行为统一和场地修复过程的多元化、集成化、通用化模拟。

2. 研究方法

2.1. 铬污染迁移修复过程模拟引擎–自适应仿真外挂系统耦合计算架构

铬污染迁移修复过程模拟引擎–自适应仿真外挂系统计算架构(图1)的本质是将复杂问题简单化，将同时涉及物理对流——弥散、氧化还原、吸附解吸、沉淀溶解、修复等多个复杂变化的铬污染迁移、转化、修复过程，分解为一个采用静态模型参数的铬污染迁移修复动力学模型和多个可变参数时空模型，以多进程同步的方式实现铬污染的多核模拟。模拟主引擎以三维地质结构数据、三维属性数据、模型参数数据以及边界条件数据为基础，构建铬污染迁移修复过程动力学耦合模型，采用有限差分的方法求解并计算铬污染迁移修复过程数值模型，并将其模拟结果三维可视化。自适应仿真外挂系统以可变参数时空模型为核心，将模型参数可变功能从模拟主引擎中剥离出来，当外部环境发生变化时，自适应仿真外挂系统根据相应的可变参数时空模型计算得到当前的模型可变参数数据，同时建立与模拟引擎之间的交互通道，实时、高效地修正土壤地下水系统中铬污染迁移修复三维动力学模型的模型可变参数，从而使得模拟主引擎的铬污染迁移修复三维动力学模型能适用于不同环境，大大增强了模型的可靠性和适应性，实现了土壤–地下水异质系统行为和铬渣遗留场地修复过程的多元化、集成化、通用化模拟。

Figure 1. Simulation main engine-adaptive simulation plug-in system coupled computational architecture

图1. 模拟主引擎–自适应仿真外挂系统耦合计算架构

模拟主引擎支持与多个自适应仿真外挂系统以松耦合的方式连接，铬污染迁移转化模拟主引擎中的模型基本参数和静态模型参数表征需进行铬污染物模拟区域的整体环境信息，而由自适应仿真外挂系统中的模型可变参数时空模型计算得到的模型可变参数表征当前时间段铬污染模拟区域的详细环境信息，二者协同工作，共同完成污染场地内复杂的铬污染迁移转化过程模拟。铬污染迁移转化模拟主引擎与自适应仿真外挂系统之间相互独立，两个系统内部存在各自的模拟内核，模拟主引擎可在无自适应仿真外挂系统连接的情况下独立运行，同时，自适应仿真外挂系统的异常不会影响模拟主引擎的正常运行，相应的模拟主引擎的异常也不会影响自适应仿真外挂系统的正常运行。铬污染迁移修复过程模拟主引擎与自适应仿真外挂系统可以在同一设备上运行，也可以在同一局域网下的不同设备上运行。

2.2. 铬污染迁移修复过程动力学耦合模型

基于铬污染场地在地下水系统中的迁移与修复特征，综合考虑物理对流——弥散、氧化还原、吸附解吸、沉淀溶解、降水补给、微生物电化学修复、抽水井修复以及原位注入修复等多个过程，构建了地下水铬污染物迁移与修复动力学耦合模型，明确了动力学模型可变参数。

$E\frac{\partial h}{\partial t}=\frac{\partial }{\partial x}\left(K_{xx}F\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}F\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{zz}F\frac{\partial h}{\partial z}\right)+\frac{Q}{\Delta x\Delta y\Delta z}-\frac{Q{}_s}{\Delta x\Delta y\Delta z}$ (1)

$\theta R\frac{\partial C}{\partial t}=\frac{\partial }{\partial x_i}\left(\theta D_{ij}\frac{\partial C}{\partial x_j}\right)-\frac{\partial }{\partial x_i}\left(q_{i}C\right)+\frac{Q{C}_s}{\Delta x\Delta y\Delta z}-f\left(t\right){\lambda }_1{\theta }_{1}C-\frac{Q{}_{s}C}{\Delta x\Delta y\Delta z}-\alpha \lambda {\overline{C}}_r$ (2)

$q=K\frac{\Delta h}{\Delta d}$ (3)

其中，E表示单位储水率；h表示水头值；K_xx、K_yy、K_zz表示x、y、z方向的渗透系数；F表示承压含水层厚度；θ表示含水层有效孔隙度；R表示阻滞因子；D表示弥散系数；C表示溶质浓度；q表示流速；t表示时间；Q表示单位时间流入地下水系统的流体体积；C_s表示单位时间流入地下水系统的流体污染物浓度；f(t)表示微生物电化学原位阻隔的作用能力，其满足f(t₀) = 1，f(∞) = 0；λ₁表示微生物电化学原位阻隔反应速率常数；θ₁表示微生物电化学原位阻隔的有效孔隙度；Q_s表示抽水速率；α表示原位注入反应系数；λ表示原位注入反应速率常数；${\overline{C}}_r$ 表示原位注入药剂浓度；∆h表示水头损失量；∆d表示渗流路径长度。

铬污染迁移修复过程动力学耦合模型中的参数均可通过自适应仿真外挂系统进行实时修正。在铬污染迁移修复模拟过程中，预设Q、C_s、λ₁、θ₁、Q_s、α、λ、${\overline{C}}_r$ 等参数值为零，若当前时间存在降水补给过程，则开启降水补给类自适应仿真外挂系统，计算当前时间的Q和C_s，并完成与模拟主引擎的连接，实时修正模型相应参数的值；若当前时间存在微生物电化学原位阻隔修复过程，则开启微生物电化学原位阻隔修复类自适应仿真外挂系统，计算当前时间的λ₁和θ₁，并完成与模拟主引擎的连接，实时修正模型相应参数的值；若当前时间存在抽水井修复过程，则开启抽水井修复类自适应仿真外挂系统，计算当前时间的Q_s，并完成与模拟主引擎的连接，实时修正模型相应参数的值；若当前时间存在原位注入修复过程，则开启原位注入修复类自适应仿真外挂系统，计算当前时间的α、λ和${\overline{C}}_r$ ，并完成与模拟主引擎的连接，实时修正模型相应参数的值。

3. 研究区概况

研究区域如图2所示，位于华南加里东褶皱带的扬子陆块东南缘，雪峰、加里东、印支–燕山等主要地壳运动在区内均留下了明显记录，其中以印支–燕山运动影响最强烈，造就了现今展示的地质结构基本轮廓。研究区第四系覆盖层相对较厚，下伏地层主要由沉积层和岩浆岩组成。本研究区域地处以丘陵山地地貌为主，西部和西南部地势较高，东部和北部较平缓。同时，本研究区属于中亚热带季风湿润气候区域，四季分明，光照充足，气候温和，降雨充沛，严寒期短，无霜期长，平均年日照1557.9小时，年平均气温17.3℃，降雨量1327.1 mm，无霜期282天，蒸发量1283.9~1340.1 mm，风速24 m/s。

Figure 2. Location map of the study area

图2. 研究区位置图

4. 铬污染迁移修复过程三维可视化模拟

4.1. 研究区水文地质概念模型建立

根据现有资料综合考虑内部结构、渗透性能、水利特征及补给排泄条件等，对场地的水文地质状况进行概化建立概念模型，构建三维地层模型，确定初始渗流场模型和初始浓度场模型，同时对边界条件和源汇项进行概化和计算。

4.1.1. 研究区概念模型及三维地层模型

根据现有资料，建立研究区铬污染迁移修复过程概念模型，如图3所示。研究区为典型的傍河型场地，铬污染主要集中在铬渣堆存区以及铬渣生产区，综合原位注入技术、抽水井技术以及微生物电化学原位阻隔技术，对研究区进行联合修复。结合钻探数据，将本研究区域污染场地内埋藏的地层概化为人工填土层，粉质黏土层，中砂圆砾层，强风化泥岩，中风化泥岩5个地层，建立研究区三维地层模型，如图4所示，其地层特征自下而上分叙如下：

1) 人工填土层：孔隙度较大，透水性佳，为强透水层，部分区域含上层滞水，主要接受大气降水补给，水量、水位因季节而异，未形成稳定连续的水位面。

2) 粉质黏土层：透水性较弱至不透水，不含水，作为隔水层存在。

3) 中砂圆砾层：呈湿至饱和状态，分选性差，孔隙发育，透水性强，为强透水层，与涟水存在水力联系，枯水季节场地内地下水向涟水补给，丰水季节涟水向场地地下水补给。其上下两层均为隔水层，具有一定的承压性，作为含水层存在。

4) 强风化泥岩：节理裂隙微弱，不透水，大多呈粉土状风化状态，作为隔水层存在。

5) 中风化泥岩：节理裂隙不发育，岩质较软，岩体较完整，岩石断面干燥，不透水，作为隔水层存在。

Figure 3. Conceptual model of chromium pollution migration remediation process in the study area

图3. 研究区铬污染迁移修复过程概念模型

Figure 4. Three-dimensional stratigraphic model of the study area

图4. 研究区三维地层模型

4.1.2. 研究区初始条件及边界条件

研究区渗流场总体流向为西北到东南，形成的流场如图5所示。将场地的西北边界设为流入边界，水头值为55~56 m，场地东南边界为流出边界，水头值为46~47.5 m。由于含水层有一定埋深，可忽略蒸散作用，考虑含水层的降雨入渗补给，场地平均年降水量为1387.4 mm，降水入渗补给系数0.179，计算得到降水入渗补给量为248.3 mm/a。

根据对研究区场地检测报告以及相关文献资料的整理分析，我们插值得到了场地地下水中铬污染物浓度的三维分布图，如图6所示。场地表层土壤(图7)的重污染区域位于三分厂、东西铬渣场和原料堆积场，铬污染浓度高达3000 mg/Kg。在降水淋滤作用下，人工填土层具有较大的渗透系数，导致其对铬污染物阻隔效果有限，铬污染物溶质向下迁移至粉质粘土层。粉质粘土层(图8)的渗透系数较小，导致铬污染物在其土层上部出现聚集现象，形成浓度较小的高峰。在长期的降水淋滤作用下，铬污染物在粉质粘土层厚度最薄处向圆砾层迁移，污染了承压含水层。为了考虑这种垂向迁移过程，在数值模拟中设置垂向污染物补给边界，其补给浓度由降水量和垂向迁移过程共同决定，会随时间的变化而变化，同时将场地边界设置为Cr(VI)溶质自由流出边界。

Figure 5. Map of the initial flow field in the study area

图5. 研究区初始流场图

Figure 6. Initial distribution of chromium contamination concentrations in groundwater

图6. 地下水中铬污染浓度初始分布图

4.1.3. 研究区水文地质参数

结合场地实验结果和相关文献资料，得到研究区域污染场地参数表，包括渗透系数、有效孔隙度、弥散系数、储水率等，如表1所示。

Figure 7. Distribution of soil chromium contamination concentrations in artificial fill layers

图7. 人工填土层土壤铬污染浓度分布图

Figure 8. Distribution of soil chromium contamination concentrations in the chalky clay layer

图8. 粉质黏土层土壤铬污染浓度分布图

Table 1. Table of hydrogeologic parameters of the study area

表1. 研究区水文地质参数表

	人工填土层	粉质黏土层	中粉砂圆砾层	泥岩层
水平渗透系数(m/d)	3.68 × 10^−^2	1.56 × 10^−^8	4.07	8.64 × 10^−^11
垂直渗透系数(m/d)	3.68 × 10^−^3	1.56 × 10^−^9	4.07 × 10^−^2	8.64 × 10^−^13
弥散系数(m^2⁄d)	8.64 × 10^−^2	8.64 × 10^−^6	7.45 × 10^−^1	8.64 × 10^−^11
孔隙度	0.78	0.64	0.46	0.21
储水率(1/m)	4.9 × 10^−^4	6.9 × 10^−^4	1.2 × 10^−^4	3.3 × 10^−^6
阻滞因子	2	10	2	2

4.2. 铬污染联合修复过程模拟

为了验证铬污染迁移修复过程模拟引擎–自适应仿真外挂系统计算架构的可行性和准确性，以研究区三维地质结构数据、三维属性数据、模型参数数据以及边界条件数据为基础，结合研究区实际情况，在研究区开展地下水铬污染迁移修复过程模拟。模拟开始时间为2020年3月1日到2023年5月1日，模拟总时长38个月，模拟步长为1天，在2020年3月1日实施原位注入修复，在2022年1月1日采取抽水井修复技术，在2022年5月1日加入微生物电化学原位阻隔技术，同时考虑土壤–地下水异质系统的行为统一，在雨季进行降水补给。修复手段位置如图9所示。

Figure 9. Location map of restoration means

图9. 修复手段位置图

研究区地下水铬污染迁移修复过程模拟开始时间为2020年3月1日，原位注入修复手段运行的时间也为2020年3月1日，首先开启自适应仿真外挂系统，选择原位注入类型的可变参数时空变化模型，导入原位注入文件并计算相关可变参数数据，选择即时修正模式，打开数据交互传输通道，等待模拟主引擎的连接，然后开启模拟主引擎，设置好相应的初始条件后开始模拟，成功与自适应仿真外挂系统建立数据交互通道，接受外挂系统传输的模型可变参数数据，对模型可变参数进行修正，完成原位注入修复模拟；保持模拟主引擎处于数值模拟状态，开启抽水井类型的自适应仿真外挂系统，导入抽水井文件，计算受抽水井影响的模型可变参数数据，选择定时修正模式，设置修正时间为2022年1月1日，建立与主引擎之间的连接，模拟主引擎接受外挂系统传输的抽水井类可变参数数据，当模拟进行到预设修正时间时，对模型可变参数进行修正，完成抽水井修复过程模拟；随后开启微生物电化学原位阻隔类自适应仿真外挂系统，构建并计算微生物电化学原位阻隔类可变参数时空模型，选择定时修正模式，设置修正时间为2022年5月1日，将可变参数数据传输至模拟主引擎，当模拟主引擎中的数值模拟进行到指定修正时间时，修正当前时间的模型可变参数，完成微生物电化学原位阻隔修复过程模拟；同时，考虑降水补给的影响，开启多个自适应仿真外挂系统，分别在雨季和旱季对模拟主引擎中的降水补给类模型可变参数进行修正，完成降水补给过程模拟；当铬污染迁移修复模拟过程结束后，模拟主引擎自动断开与所有自适应仿真外挂系统之间的连接。模拟主引擎与多个自适应仿真外挂系统协同合作，共同完成研究区铬污染联合修复过程模拟。

研究区地下水铬污染物联合修复过程模拟结果如图10所示。研究区地下水中铬污染物从整体上由西北向东南迁移，2020年3月至2020年7月属于雨季，地下水中铬污染物浓度峰值由156.44 mg/Kg增大至177.27 mg/Kg，铬污染物浓度中心位置基本保持不变；2020年8月至2021年2月属于旱季，地下水中铬污染物浓度中心向铬渣遗留场地东南角移动，铬污染物浓度最大值由178.79 mg/Kg下降至169.45 mg/Kg；2021年3月到2021年7月属于雨季，在此期间，地下水中铬污染物浓度呈上升趋势，2021年3月地下水中铬污染物浓度的最大值为168.70 mg/Kg，2021年7月地下水中铬污染物浓度的最大值为183.13 mg/Kg，铬污染物浓度中心范围有所扩大，其位置未发生明显偏移；2021年8月到2022年2月属于旱季，地下水中铬污染物浓度峰值降低，由184.86 mg/Kg变为180.59 mg/Kg，浓度中心缓慢向场地东南角移动；2022年3月至2022年7月属于雨季，在降水淋滤的作用下，土壤中的铬污染物不断进入地下水系统，但由于抽水井修复技术的作用，地下水中污染物浓度有所降低，由179.72 mg/Kg变为179.13 mg/Kg，铬污染物浓度中心位置基本不变；2022年8月至2023年2月属于旱季，地下水中铬污染物浓度中心位置发生偏移，逐渐靠近场地东南角，浓度峰值呈下降趋势；2023年3月至2023年5月属于雨季，地下水中铬污染物浓度中心范围扩大，浓度峰值由2023年3月的177.37 mg/Kg增长至182.62 mg/Kg。2020年3月1日采取原位注入修复手段后，一期场地位置铬污染物浓度成功降低至标准范围，并且有效的阻止了二期场地的铬污染物向一期场地迁移；2022年1月1日开始在研究区实施抽水井修复，抽水井布设位置的地下水铬污染物浓度显著降低，对铬污染物浓度的上升和扩散有一定的遏制作用；2022年5月1日加入微生物电化学原位阻隔修复，原位阻隔布设位置的地下水铬污染物浓度显著降低，有效地阻止了地下水中铬污染物浓度中心向下游移动。

Figure 10. Simulation results of chromium contamination co-remediation process

图10. 铬污染联合修复过程模拟结果图

5. 结论

本文针对铬渣遗留场地三维结构及参数的时空异质性和修复技术的综合多样性，提出了铬污染迁移转化与生物修复过程模拟引擎–自适应外挂耦合计算架构及仿真范式。本文研究了铬污染迁移转化和修复过程中对流–弥散–生物地球化学反应之间的耦合关系，建立了铬污染迁移修复过程动力学耦合模型，以铬渣遗留场地三维结构模型和三维属性模型为数字化基础模型，利用铬污染迁移转化与修复过程模拟引擎–自适应外挂耦合计算架构，实现了土壤–地下水异质系统的行为统一和场地修复过程的多元化、集成化、通用化模拟。本文主要研究成果如下：

1) 本文综合考虑物理对流——弥散、氧化还原、吸附与解吸、沉淀与溶解、植物与微生物反应、降水补给、微生物电化学修复、抽水井修复以及原位注入修复等多个过程，构建了地下水铬污染物迁移与修复动力学耦合模型。

2) 本文以铬污染物迁移与修复动力学耦合模型为基础，建立了铬污染迁移转化与生物修复过程模拟引擎–自适应外挂耦合计算架构及仿真范式，有效解决了由铬渣遗留场地三维结构及参数时空异质性和修复过程综合多样性导致的参数可变问题，为铬污染时空分布的研究以及铬污染场地修复方案的优化提供了技术支持。

3) 本文以某铬渣遗留场地为例，结合场地勘探数据，实验室数据以及监测数据，利用铬污染迁移修复过程模拟引擎–自适应外挂耦合计算架构，进行铬污染迁移修复过程模拟，在降水补给的作用下，粉质黏土层土壤中的铬污染物不断进入地下水系统，并在对流–弥散作用下向东南方向迁移，同时，在原位注入修复技术、抽水井修复技术和微生物电化学原位阻隔修复技术的作用下，局部地下水中铬污染浓度显著降低，有效地遏制了铬污染物的进一步扩散，模拟结果与研究区实际情况基本吻合。

基金项目

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

NOTES

^*通讯作者。