中国石油、化工等重点行业遗留场地土壤-地下水系统（土-水系统）中有机污染场地多、范围广、面积大、污染重^[1-4]。有机污染物大多微溶于水，常以非水相流体（non-aqueous phase liquids，NAPLs）形式存在，在全生命周期污染过程中NAPLs以多相形式（自由相、溶解相、吸附相和气相）赋存于土-水介质中，且缓慢而持续地释放难以降解的高毒性物质，对公众健康和生态安全造成极大威胁^[5-6]。场地土-水系统污染修复和风险管控是国家重大战略需求，然而当前自然修复缓慢、人工修复难度大，已引起行业的高度重视^[7]。针对国家加强土壤-地下水生态环境保护与污染风险管控的需求，为实现场地土-水系统中有机污染物时空分布与演化变化过程刻画^[8-9]、健康风险评估和生态环境治理与防控^[10-11]，数值模拟与可视化是其中关键环节和量化工具，亟需开展场地土-水系统中有机污染物时空分布模拟与可视化研究，以提升场地污染防控与管理能力^[12-14]。

场地有机污染物在非均质土-水系统中形成复杂的多相流动体系，其时空分布驱动机制复杂，涉及“温度场-水动力场-化学场”多场耦合、“气-水-油”多相态、场地复合型污染的多组分以及有机污染物在多介质中复杂的化学衰减过程。目前美国劳伦斯伯克利实验室开发的TMVOC（the Simulator for Multiple Volatile Organic Chemicals）程序在有机污染物多相分配和迁移转化模拟中应用最为广泛，其重点考虑多相流动力学机制^[15-17]，然而，未考虑土-水介质中有机污染物多场耦合驱动作用，场地尺度模拟预测精度低，难以满足污染风险精准管控的需求。亟需在多相流模拟技术基础上构建多相流场-温度场-化学场多场耦合的数值模型，提高刻画精度，这是科学、规范开展土水污染绿色低碳修复管控的前提，亦具有很强的挑战性。

场地环境信息化是环保行业创新驱动下的一个热门研究领域，通过数据库技术和地理信息技术相结合，将传统场地污染调查与监测数据信息化和可视化^[8]。许亚飞^[18]基于ArcGIS Server架构构建持久性有机污染物污染场地信息管理系统，实现了污染场地信息的采集与编辑、查询与统计、空间地图服务、危害等级评价、优先治理名录筛选、场地档案资料管理、报表统计分析、可视化输出及系统维护等功能；毕天平等^[19]基于网络地理信息系统WebGIS（Web Geographical Information System）开发突发性环境污染公共服务平台，实现了环境污染信息和扩散效果可视化，并可根据公众自身安全的风险源分布情况及特点进行风险评估及模拟分析；尤其浩^[20]开发了基于WebGIS的土壤污染评价与三维可视化系统，开展土壤污染动态图表统计、污染评价与污染场地三维可视化，为土壤污染防治、预警、修复治理以及相关政策提供帮助和支持。当前系统的主要功能是实现场地污染现状监测数据的集成管理、可视化和现状污染评价等。然而，土-水系统中污染物的时空分布是变化的，只有在充分掌握污染过程和变化趋势的条件下，方可实现污染风险的科学评价和绿色低碳修复管控。随着互联网技术发展，WebGIS技术和网络图形库WebGL（Web Graphics Library）图形可视化技术功能愈发强大，亟需在现有的场地环境信息系统功能基础上嵌入有机污染物驱动过程的数值模型，开发集“场地调查与监测数据集成管理、有机污染物时空分布图形化建模和污染全过程场景可视化表达”功能于一体的浏览器/服务器（Browser/Server，B/S）架构的场地土-水系统有机污染物时空分布模拟与可视化系统（以下简称“可视化系统”），以提升场地有机物污染的综合管理与决策能力^[21-23]。

1 材料与方法 1.1 系统设计

可视化系统采用WebGIS地理信息系统为基础平台，前端采用超文本标记语言HTML（Hyper Text Markup Language）、层叠样式表CSS（Cascading Style Sheets）、JavaScript编程语言、Cesium可视化技术进行页面管理和可视化，后端嵌入多源数据三维空间插值、有机污染多场耦合模型等内核计算程序，采用Spring框架技术处理用户发送的请求。系统设计“数据管理与对接”、“有机污染物时空分布图形化建模”和“场景可视化”三大功能模块，实现有机污染场地土-水系统“调查和监测-模拟评价-三维可视化”全链条设计与智慧化管理。可视化系统主体功能和结构分如图 1所示。

1.1.1 数据管理与对接模块

系统采用WebGIS数据管理可视化系统为基础平台，嵌入超文本传输协议HTTP（Hypertext Transfer Protocol）标准的应用体系和开源的客户/服务器关系型数据库管理系统PostgreSQL，实现服务器端和客户端（浏览器）信息互联和发布。数据管理与对接模块功能的通用性体现在：①通过自主定义数据表头格式可实现场地污染调查和监测多源数据（包括钻孔、地层岩性、土壤 & 地下水参数、多组分监测浓度、水位等）的存储、查询与管理；②数据格式自由转换，可实现底层数据与图形化建模的对接；③对接场景可视化模块，实现多类型数据的平台展示。

1.1.2 有机污染物时空分布图形化建模模块

系统前端采用模型-视图-控制器MVC（Model View Controller）编程模式开发通用的在线创建数值模型的图形用户界面，作为独立模块嵌入系统平台，实现模型创建与基础数据实时交互，建模过程图形化和界面化的便捷编辑和输入；系统后端嵌入自主开发的场地土-水系统有机污染物时空分布多场耦合模拟内核程序（详见1.1.3）。

1.1.3 场景可视化

采用Cesium三维渲染引擎对空间信息进行二维地图/三维实体可视化渲染和场景可视化，可视化对象包括场地基础调查与监测多源异构数据、在线建模过程数据、多情景多场耦合模拟预测结果等，以支撑面向场地污染调查、监测、污染运移过程的动态模拟与仿真。

1.2 土-水系统有机污染物时空分布多场耦合模拟技术

土-水系统中有机污染物的赋存状态受地层非均质性控制，呈现“水-气-NAPLs多相态共存、多介质空间变异、持续释放”等特征，典型的NAPLs迁移转化过程包括：①NAPLs自由相在重力作用下的多相流动；②NAPLs溶解相的对流-弥散（图 2）。NAPLs迁移转化还受到温度和化学作用的影响，化学作用包括NAPLs相吸附、挥发和降解等生物化学衰减作用，温度会改变NAPLs理化参数和化学衰减过程的控制参数，进而影响NAPLs迁移转化过程。因此，NAPLs在土-水系统中的多相分配和迁移转化过程本质是多相流场-温度场-化学场多场耦合驱动的结果。

基于此，本文在多相流模拟程序的基础上，开发土-水系统有机污染物时空分布多相流场-温度场-化学场多场耦合模型，并作为内核程序嵌入可视化系统，以期实现场地有机污染物时空分布的模拟与仿真。多场耦合控制方程如下：

1.2.1 多相流动力场连续性方程

考虑土-水系统中水和NAPLs的流动过程。

$-\left[\frac{\partial v_{\beta x}}{\partial x}+\frac{\partial v_{\beta y}}{\partial y}+\frac{\partial v_{\beta z}}{\partial z}\right]=S_\beta \frac{H_\beta}{\mathrm{t}}$

(1)

${v_\beta } = - k\frac{{{k_{\gamma \beta }}}}{{{\mu _\beta }}}\left( {\nabla {p_\beta } - {\rho _\beta }g} \right)$

(2)

${p_\beta } = p + {p_{c\beta }}$

(3)

式中，x，y，z为各变量（如渗流速度）在三个坐标轴上的分量；v_β为β相的流动速度，L·T^–1；S_β为β相饱和度；t表示时间，T；H表示水头，L；k为总渗透率，L²；k_γβ表示β相的相对渗透率（处于0到1之间）；μ_β为β相黏滞系数，M·（L·T）^–1；ρ_β为β相密度，M·L^–3；g重力矢量，L·T^–2；p_β为β相流体压力，M·L^–3，可表示为β相中压力p和毛细压力p_cβ之和。

1.2.2 温度场连续性方程

考虑土-水系统中温度的变化和热量传输过程。

$\begin{aligned} C \frac{\partial T}{\partial \mathrm{t}}= & \frac{\partial}{\partial x}\left(\lambda_x \frac{\partial T}{\partial x}\right)+\frac{\partial}{\partial y}\left(\lambda_y \frac{\partial T}{\partial y}\right)+\frac{\partial}{\partial z}\left(\lambda_z \frac{\partial T}{\partial z}\right)- \\ & \mathrm{C}_{\mathrm{w}} \frac{\left(v_x T\right)}{\partial x}-\mathrm{C}_{\mathrm{w}} \frac{\left(v_y T\right)}{\partial y}-\mathrm{C}_{\mathrm{w}} \frac{\left(v_z T\right)}{\partial z}+f_1 \end{aligned}$

(4)

式中，x，y，z为各变量（如渗流速度）在三个坐标轴上的分量；T为温度，K；C和C_w分别为多孔介质和水的热容量，J·（L³·K）^–1；v为流动速度，L·T^–1；$\lambda$为热动力弥散系数，J·（L·T·K）^–1；f₁为单位时间单位体积含水层内由于生化反应等作用引起的温度变化。

1.2.3 化学场连续性方程

考虑NAPLs溶解相的对流-弥散和NAPLs相吸附、挥发和降解等化学衰减作用。

$\begin{aligned} \frac{\partial S}{\partial t}= & \frac{\partial}{\partial x}\left(\mathrm{D}_x \frac{\partial S}{\partial x}\right)+\frac{\partial}{\partial y}\left(\mathrm{D}_{\mathrm{y}} \frac{\partial S}{\partial \mathrm{y}}\right)+\frac{\partial}{\partial \mathrm{z}}\left(\mathrm{D}_{\mathrm{z}} \frac{\partial S}{\partial z}\right)- \\ & \frac{\left(v_x S\right)}{\partial x}-\frac{\left(v_y S\right)}{\partial y}-\frac{\left(v_z S\right)}{\partial z}+f_2 \end{aligned}$

(5)

式中，x，y，z为各变量（如渗流速度）在三个坐标轴上的分量；S为溶质的浓度，M·L^–3；D为水动力弥散系数，L²·T^–1；v为流动速度，L·T^–1；f₂为单位时间单位体积含水层内由吸附、挥发、溶解和降解等作用引起的溶质质量的变化，本研究采用综合化学衰减系数X来表征。

$\mathrm{X}=\left(\sigma_1 K_{\mathrm{V}}+\sigma_2 K_{\mathrm{S}}+\sigma_3 K_{\mathrm{A}}+\sigma_4 K_{\mathrm{bio}}\right) \cdot M$

(6)

$K_{\mathrm{V}}=\mathrm{VP}_{\mathrm{A}}+\mathrm{VP}_{\mathrm{B}} \cdot T+\mathrm{VP}_{\mathrm{C}} \cdot T^2+\mathrm{VP}_{\mathrm{D}} \cdot T^3$

(7)

$K_{\mathrm{S}}=\mathrm{SOL}_{\mathrm{A}}+\mathrm{SOL}_{\mathrm{B}} \cdot T+\mathrm{SOL}_{\mathrm{C}} \cdot T^2+\mathrm{SOL}_{\mathrm{D}} \cdot T^3$

(8)

$K_{\mathrm{A}}=K_{\mathrm{oc}} \cdot f_{\mathrm{oc}}$

(9)

式中，M为污染物质量分数，mg·kg^–1；T为温度，℃；K_V为挥发系数（VP_A、VP_B、VP_C、VP_D为挥发系数经验参数）；K_S为溶解系数（$\text{SO{L_A}}$、$\text{SO{L_B}}$、$\text{SO{L_C}}$、$\text{SO{L_C}}$为溶解系数经验参数）；K_A为吸附系数；${K}_{oc}\;\text{n}\;{f}_{oc}$为土壤介质中有机质质量分数/水与有机物之间的分离系数；σ表示各化学参数的权重。

1.2.4 多相流场-温度场-化学场之间的耦合关系

① 多相流场与温度场通过对流弥散产生的热量迁移进行耦合，从数学方程的角度则利用流速v进行变量传递，利用黏滞系数μ和温度之间的量化关系实现多相流与温度的耦合；②多相流场与化学场则通过对流产生的溶质运移进行耦合，从定量刻画的角度则是利用流速v进行变量传递，并利用溶质浓度S、密度ρ和黏滞系数μ之间的量化关系实现多相流与溶质运移之间耦合；③温度场与化学场方程通过综合化学衰减系数X和温度之间的函数关系进行耦合，耦合关系式和关键参数如图 3所示。

1.2.5 模型求解

在多相流模拟源程序基础上，采用序列耦合多孔介质温度场模块和化学场模块求解有机污染物迁移多相流-温度-化学多场耦合控制方程，程序设计流程如图 4所示。

建模流程包括刻画地层结构、模型分层与网格剖分、分区与材料属性赋值、模型初始条件和边界条件设置、模型参数和计算参数设置、模拟结果输出与可视化等。可视化系统采用MVC编程模式开发在线创建模型的图形用户界面，以实现以上建模过程的图形化操作。

2 结果与讨论 2.1 场地概况

应用场景选择南方某遗留化工污染场地，曾从事甲醛溶液、乙酸乙酯、乙酸丁酯的生产。钻孔揭露的浅部地层可概化为五层：①素填土层：土质松散，渗透性较大，为地块非饱和-饱和潜水含水层，厚度0.5~4 m；②粉质黏土层：渗透性较小，为相对弱透水层，地块内连续分布；③粉砂层和④中砂层：渗透性相对较大，含水量丰富，为地块主要含水层，连续分布于地块中部低洼的冲沟附近，厚度1~4 m；⑤砂岩层：基岩体渗透性与第②层相仿，为相对隔水层，地块内连续分布。

经场地土-水污染初调和详调与分析，邻苯二甲酸二（2-乙基己基）酯（di-（2-Ethylhexyl）Phthalate，DEHP）在平面上已形成连片污染，局部回填型污染源导致场地北部与东北部区域形成高浓度污染物聚集。图 5联合展示了场地土-水系统地层结构和DEHP污染物空间分布，NAPL相污染物在垂向上主要分布于第①~③层，溶解相污染物主要分布于第③~④层。DEHP是一种轻非水相有机污染物（LNAPL），化学式为C₂₄H₃₈O₄，用作聚氯乙烯（PVC，Polyvinyl chloride）等塑料制品的增塑剂，2017年世界卫生组织国际癌症研究机构将其列为2B类致癌物清单中，本文选取DEHP为场地特征污染物开展模拟与可视化应用与验证。

2.2 系统数据管理与可视化

可视化系统设置了项目管理、钻孔管理、监测点位、模型展示、模拟预测、二维展示、三维展示和统计分析八大模块。底层数据库可对场地钻探、多源监测数据进行存储管理，并通过二维地图和三维实体进行二维/三维展示。二维地图可视化即以百度/天地图卫星地图为底图展示场地边界、多类型监测点位、二维等值线等；三维实体采用Cesium渲染引擎展示钻孔柱状信息、地层结构、三维流场/污染羽等，并可弹窗形式显示点、线、面属性信息。系统界面及可视化效果如图 6所示。

2.3 土-水系统有机污染物时空分布在线建模与分析

基于系统开发的图形用户界面组件，在场地地质结构模型构建的基础上，按照网格生成与加密、地层分区、材料属性添加、初始条件设置、边界条件设置、参数设置、工况计算的流程在线创建土-水系统有机污染物时空分布多相流场-温度场-化学场多场耦合模型，模拟有机污染物时空分布和预测污染演化趋势。基于场地条件和调查数据，模型的相关条件和参数设置如下：

2.3.1 边界条件和初始条件设置

对于饱和带与非饱和带的初始条件分别进行如下设置：

$\left\{\begin{array}{l} \text { 气相饱和度: } S_{\mathrm{G}}=0.3(\text { 非饱和带 }) S_{\mathrm{G}}=0(\text { 饱和带 }) \\ \text { NAPLs相饱和度: } S_O=0 \\ \text { 水相饱和度: } S_{\mathrm{w}}=1-S_{\mathrm{G}}-S_O \\ \text { 溶解相、挥发相NAPLs浓度 }: c_1=0 c_2=0 \\ \text { 温度: } T_0=25{ }^{\circ} \mathrm{C}(\text { 饱和带与非饱和带 }) \end{array}\right.$

模型的边界条件如下设置：

① 多相流场

$\left\{\begin{array}{l} \text { 上边界: } \nabla S_{\mathrm{G}}=0, \nabla S_{\mathrm{W}}=0 \\ \text { 左右边界: } \nabla S_{\mathrm{G}}=0, S_{\mathrm{WI}}=0.3(\text { 非饱和带 }) \\ S_{\mathrm{G}}=0, S_{\mathrm{W}}=1(\text { 饱和带 }) \\ \text { 下边界: } S_{\mathrm{G}}=0, \nabla S_{\mathrm{W}}=0 \\ \text { 泄漏边界: } V_{\mathrm{G}}=3 \times 10^{-5} \mathrm{~kg} \cdot \mathrm{~s}^{-1} \end{array}\right.$

② 化学场

$\left\{\begin{array}{l} \text { 上边界: } S_0=0 \\ \text { 左右边界: } S_0=0 \\ \text { 下边界: } S_0=0 \\ \text { 下边界: } S_{\mathrm{G}}=0, \nabla S_{\mathrm{w}}=0 \\ \text { 泄露边界: } V_0=3 \times 10^{-5} \mathrm{~kg} \cdot \mathrm{~s}^{-1} \end{array}\right.$

③ 温度场

$\left\{\begin{array}{l} \text { 上边界: } T=\text { 各季节平均温度 } \\ \text { 左右边界: } \nabla T=0 \end{array}\right.$

2.3.2 模型参数设置

多场耦合模型中含水介质、地下水、气、DEHP理化参数设置如表 1所示。介质的水动力参数（渗透率、孔隙度）来自场地水文地质勘察报告；土、水、气理化性质参数（比热容、导热系数、密度等）以及DEHP的化学衰减（吸附、挥发、溶解等）作用的初始参数和多相扩散系数参考TMVOC说明书^[16]和相关文献^[6，24-29]，其中DEHP化学衰减作用参数进一步通过数值模型拟合反演而得。

2.3.3 模拟分析

基于以上设置的模型边界条件、初始条件和模型参数，经多场耦合模拟得到场地土-水系统中水动力场、温度场、化学场（综合化学衰减系数）和DEHP浓度分布（图 7）。采用拟合优度R值作为评估指标：

$R=\frac{\sum\limits_{i=1}^N\left(\mathrm{cal}_i-\mathrm{cal}\right)\left(\mathrm{obs}_i-\mathrm{o}\overline{\mathrm{b}}\mathrm{s}\right)}{\sqrt{\sum\limits_{i=1}^N\left(\mathrm{cal}_i-\mathrm{cal}\right)^2} \sqrt{\sum\limits_{i=1}^N\left(\mathrm{obs}_i-\mathrm{o}\overline{\mathrm{b}}\mathrm{s}\right)^2}}$

(10)

式中，cal_i为计算值，obs_i为观测值，$\mathrm{c}\overline{\mathrm{a}}\mathrm{l}$与$\mathrm{o}\overline{\mathrm{b}}\mathrm{s}$表示计算值与观测值的平均值。采用场地多点位、多期监测数据，包括2020年6月完成的场地初调和详调数据及后续三期污染监测数据（时间分别为2021年6月、2022年9月和2023年3月），对模型进行拟合校准，整体模拟拟合优度R达0.91，模拟误差低于30%（图 8）。

土-水系统中DEHP空间分布模拟结果显示，地下水位位于第③层粉砂层底部，DEHP为LNAPL污染物，NAPL相污染物在垂向上主要分布于第①~③层土壤中，第③层粉砂层和第④层中砂层为渗透性相对较好的含水层，DEHP从第③层向下迁移，NPAL相在第③层底部赋存，DEHP溶解组分在第③~④层中扩散。DEHP于土壤中的总体积为3.3 m³，总质量为3.2×10³ kg，平均浓度为96.2 mg·kg^–1；DEHP于地下水中的最高浓度5 μg·L^–1，总质量为1.3×10^–2 kg，DEHP在重力及溶解作用下有0.04‰进入地下水中。

考虑水动力条件、温度变化和化学衰减作用等进行多种情景模拟，利用模拟结果与场地实测值进行对比，进一步探究多相流场-温度场-化学场多场耦合对DEHP迁移影响，识别主控因素。

水动力条件：模型校正后地下水流场如图 7a所示，由北向南平均水力梯度为0.5%，DEHP运移的优势通道与地下水流方向相同，地下水流速增大则流体动压增强，进而促进LNAPL相流体水平运移。利用模型进行敏感性分析，结果表明水力梯度每增加0.5%，污染物最大浓度点会沿地下水流方向每年平均增加运移4.68 m。

温度场：温度改变了有机污染物理化参数及其驱动过程（如吸附、挥发）控制参数，进而改变了有机污染物在土-水系统中的多相分配和迁移转化过程。依据场地多年月平均地表温度资料，设置月平均温度（1月（14.5℃）→4月（24.0℃）→7月（30.5℃）→9月（29.0℃）→12月（16.0℃））为上边界定温条件，模拟得到土-水系统温度场空间变化（图 7c）。在温度场耦合作用下，模拟结果显示在1月、4月和7月土-水系统中DEHP分别有5.8%、7.6%和11.3%的NAPL相转化为气相，3.2%、6.1%和7.4%的NAPL相转化为溶解相。

化学场：采用综合化学衰减系数来表示土-水系统中有机污染物吸附、挥发、溶解化学作用场，化学场与吸附、挥发、溶解作用控制参数和污染物浓度呈正相关关系。在化学场耦合驱动作用下，模拟期末土-水中DEHP在吸附作用下NAPL相中有50.3%成为残留相，其余NAPL自由相的衰减率达到18.7%，，其中挥发作用占11.3%，溶解作用占7.4%。

逐层模拟分析化学场、温度场和水动力场对DEHP迁移转化的影响，对比场地实际监测数据，多场耦合模型相对于分别忽略化学场、温度场和水动力场三种模拟情景，模型精度分别提升了13.05%、7.48%、4.66%（图 8），原因是场地的水动力条件相对较弱，且四季温差变化相对较小，而场地平均温度相对较高，在温度耦合作用下土-水系统中有机污染物的综合化学作用活跃。因此，考虑该应用场地实际条件，对DEHP污染物迁移转化的影响程度由高到低依次为化学场、温度场、水动力场。

2.4 场地有机污染物运移可视化

基于前述拟合校正的多相流-温度场-化学场多场耦合模型，采用现状水文年的水动力、温度和化学作用条件作为预测方案，预测时长20年。可视化系统可实现不同模拟场景下有机污染物时空分布与变化过程的可视化，图 9展示了土-水系统中有机污染物时空分布的预测结果，大图为模拟预测5年末，3张小图分别为1年末、10年末和20年末土壤中DEHP的空间分布。

场地东西向存在两个不同的DEHP污染源区，地下水流整体由场地北侧向东南和西南方向流动，污染物迁移转化模拟预测结果显示，在水动力场驱动下DEHP呈向东南和西南方向迁移扩散的趋势，土壤中DEHP含量在吸附、挥发和溶解综合化学作用下呈现不断衰减的趋势。场地土-水系统中DEHP时空分布预测结果的总体趋势合理，符合场地有机污染多相流-温度场-化学场多场耦合驱动机制，然而预测的精度尚需持续开展长效监测来验证。

3 结论

本文基于WebGIS开源架构和Cesium三维可视化技术，嵌入自主研发的场地土-水系统有机污染多相流场-温度场-化学场多场耦合模拟程序，开发B/S架构的场地土-水系统有机污染物时空分布模拟与可视化系统。可视化系统应用于南方某有机污染场地土壤-地下水中DEHP特征污染物时空分布与变化多场耦合模拟与可视化，总体拟合优度R达0.91，模拟误差低于30%。基于多场耦合模型进一步预测评估DEHP的衰减过程，得到模拟周期内NAPL相被吸附为残留相的比例（50.3%）和NAPL自由相中由于挥发（11.3%）和溶解（7.4%）作用而对应发生衰减的比例。由于现场参数获取困难，模型采用的土、水、气、NAPL的理化性质参数参考相关文献，导致场地尺度的模拟结果存在一定的不确定性，然而DEHP时空分布预测结果的总体趋势合理，符合场地有机污染多相流-温度场-化学场多场耦合驱动机制，通过可视化系统不同场景下有机污染物时空分布模拟预测和可视化，可支撑管理者科学、规范地开展土水污染绿色低碳修复管控。可视化系统可进一步拓展污染风险评估模块、污染管控模块、修复方案设计与效果对比分析等模块，以更全面和系统地推广应用于污染场地环境信息化与智慧化管理。