滨海含水层氯苯多相迁移数值模拟

Numerical Simulation of Chlorobenzene Multiphase Migration in a Coastal Aquifer

摘要

为揭示海水入侵所驱动的变密度流场对重质非水相液体（dense non-aqueous phase liquids，DNAPL）多相运移路径与相间质量交换的调控机制，本研究以氯苯为例构建变密度流-多相流-溶质运移耦合数值模型，系统模拟其入渗、再分布及溶解羽演化过程.研究结果表明：（1）NAPL相氯苯在重力主导下垂直入渗并在隔水底板形成非对称污染池；（2）咸水楔显著改变溶解相氯苯的运移路径，密度梯度驱动其沿咸淡水界面爬升，并在咸淡水过渡带积聚，且排海通量峰值较无海水入侵情景降低62%；（3）渗透系数K的增大加速NAPL相和溶解相氯苯向海洋边界迁移，同时强化溶解相在咸淡水过渡带的累积.本研究阐明了咸淡水过渡带对DNAPL迁移的截留机制，证实该区域是滨海地下水环境风险评估中不可忽视的长期次生污染源.

Abstract

To reveal the regulatory mechanism of the variable-density flow field induced by seawater intrusion on the multiphase migration pathways and interphase mass exchange of dense non-aqueous phase liquids (DNAPL), this study developed a coupled variable-density flow, multiphase flow, and solute transport numerical model using chlorobenzene as an example, systematically simulating its infiltration, redistribution, and dissolved plume evolution. The results indicate follows: (1) The NAPL-phase chlorobenzene vertically infiltrates under gravity-dominated conditions, forming an asymmetric contaminant pool above the aquitard. (2) The saltwater wedge significantly alters the migration pathway of dissolved-phase chlorobenzene, with density gradients driving its upward movement along the freshwater-saltwater interface and accumulation in the mixing zone, while the peak discharge flux to the sea decreased by 62% compared to the scenario without seawater intrusion. (3) An increase in hydraulic conductivity (K) accelerates the migration of both NAPL-phase and dissolved-phase chlorobenzene toward the marine boundary, while enhancing the accumulation of the dissolved phase within the mixing zone. This study elucidates the retention mechanism of the mixing zone for DNAPL migration and confirms that this zone represents a non-negligible long-term secondary pollution source in coastal groundwater environmental risk assessments.

Graphical abstract

关键词

滨海含水层 / 海水入侵 / 重质非水相液体 / 氯苯 / 多相流数值模拟 / 水文地质 / 环境地质.

Key words

coastal aquifer / seawater intrusion / DNAPL / chlorobenzene / multiphase flow numerical simulation / hydrogeology / environmental geology

0 引言

地下水作为滨海地区的关键淡水资源，对维持社会经济持续发展与生态系统稳定具有重要作用（Fang et al.， 2021； Hu et al.， 2024）.滨海地区地下水位浅，极易受到密集人类活动的影响，其中，以氯代烃（如三氯乙烯、四氯乙烯、氯苯等）为代表的重质非水相液体（dense non-aqueous phase liquids， DNAPL）因其高毒性、致癌性及难降解性而备受关注（李培华等，2021；Shi et al.，2023； Wu et al.， 2023；）.沿海地区往往是石油化工、港口物流及精细化工产业的集聚区，储罐泄漏或管道破裂导致DNAPL进入地下环境的风险极高（Michael et al.， 2017）.DNAPL密度大于水，进入地下后可在重力驱动下穿透包气带与含水层，直至在低渗透性隔水底板之上积聚，形成能长期释放溶解态污染物的深部污染源，对地下水水质构成持久且严重的威胁（Prommer et al.， 2002； Yang et al.， 2018； Illy et al.， 2022）.

DNAPL在含水层中的迁移是一个复杂的多相流过程，主要受毛细捕获、粘性指进及重力驱动等机制调控（Kueper et al.， 1993； Praseeja and Sajikumar， 2019）.DNAPL迁移路径与空间分布受多种因素影响，其中含水层介质的特性（如渗透率、孔隙度）起着关键作用，可导致DNAPL形成离散的高饱和度聚集区或形成广泛的污染羽（张蔚等，2013； 程洲等，2014； 邓亚平等，2015； Wu et al.， 2017a，2017b）.DNAPL沿高渗透通道优先运移，绕过低渗透透镜体，形成高度非均质的污染源区结构（Wu et al.， 2020）.此外，DNAPL自身的物理化学性质（如密度、黏度）等因素也决定了其在环境中的长期归宿与可修复性（邓亚平等，2016）.密度更高、黏度更低的DNAPL，其垂向入渗深度更深（杨宾等，2013）.

在滨海含水层中，海水入侵所导致的水动力条件改变，使不同形态污染物的迁移路径产生显著差异（Yu et al.， 2022； Zheng et al.， 2023； Gao et al.， 2024）.对于溶解态污染物，其运移主要受地下水流动控制，海水入侵会影响流速大小和方向，进而影响污染羽的对流-弥散路径（Zhang et al.， 2024）.相比之下，DNAPL的迁移存在“双路径”机制：一方面，其作为独立非水相液体的运移受重力驱动和多相流过程（如毛细捕获）控制；另一方面，DNAPL会通过溶解作用不断释放溶解相组分，形成次生污染羽，溶解相DNAPL的迁移会受到海水入侵的调控（巨梦蝶等，2025）.这种本体相与溶解相在迁移路径上的耦合机制，使得DNAPL在海水入侵背景下的环境行为远比单一溶解态污染物复杂，DNAPL在海水入侵背景下的多相运移行为却鲜有研究.

本研究采用二维剖面模型概化滨海含水层，构建了耦合海水入侵、DNAPL多相流及其溶解相运移的数值模型.通过对比存在与不存在海水入侵两种情景下DNAPL的入渗过程、相态分布及溶解相污染羽的时空演化特征，并结合渗透系数的敏感性分析，系统揭示咸淡水过渡带对DNAPL迁移与转化的控制机制.研究成果可为滨海地区地下水污染的定量风险评估与修复策略制定提供理论支撑.

1 数值模拟方法

1.1　概念模型

本研究将滨海含水层系统概化为一个理想化的二维垂向剖面（图1），设定为长500 m，高30 m.模型假设含水层为饱和均质含水层，流体和固体基质均不可压缩，水相的密度仅依赖于盐浓度.模型边界根据其水文功能设定，具体设定如下.

（1）左侧海洋边界，设定为定水头边界，水头值为25.95 m；同时定义为定浓度边界，盐浓度（C₁）恒为35 g/L.

（2）右侧淡水边界，设定为定水头边界，水头值为26.97 m；同时定义为定浓度边界，盐浓度（C₁）恒为0 g/L.

（3）顶部为定流量降雨边界，进入饱和带的降雨量为200 mm/a.在DNAPL泄漏阶段，距左边界300 m处设一点源，以4×10⁻⁴ kg/（m·s）速率注入NAPL相氯苯.

（4）底部边界设为无通量边界.

研究分为3个连续的模拟阶段：第一阶段为海水入侵阶段，模拟海水入侵直至形成稳定的流场与盐度场；第二阶段为DNAPL泄漏阶段，在第一阶段基础上，模拟DNAPL持续注入5年内的入渗、迁移和溶解过程；第三阶段为DNAPL再分布阶段，在第二阶段的结果下停止DNAPL注入，模拟其在滨海含水层长达100年的再分布过程.需要说明的是，第二阶段和第三阶段均基于第一阶段达到的稳定流场进行模拟，即假设在后续105年的模拟期内，淡水边界水头、海洋边界水头及盐度分布不随时间变化.这一稳态假设旨在聚焦于揭示变密度流场对DNAPL多相运移的核心物理机制，避免动态边界条件带来的额外复杂性.

1.2　数学模型

1.2.1　海水入侵及盐分运移模型

在海水入侵过程中，考虑含水层孔隙水为单一液相（水相），盐分在水相中的运移由对流-弥散方程控制（Yao and Zhang， 2023）：

\frac{\partial \left(\phi {\rho }_{\mathrm{w}}{C}_{\mathrm{1}}\right)}{\partial t}=\nabla \cdot \left(\phi {\rho }_{\mathrm{w}}{\mathit{D}}_{\mathrm{1}}\nabla C_{\mathrm{1}}\right)+\nabla \cdot \left[{\rho }_{\mathrm{w}}{C}_{\mathrm{1}}\frac{k}{{\mu }_{\mathrm{w}}}\left(\nabla P_{\mathrm{w}}+{\rho }_{\mathrm{w}}g\right)\right]+Q_{\mathrm{S}}{C}_{\mathrm{S}},(1)

式（1）中：D₁ 为水动力弥散系数（m²∙s^-1）；C₁为盐分浓度（g∙L^-1）；Q_s为源（汇）单位体积流量（m³∙s^-1）；C_s为源（汇）流中溶质浓度（g∙L^-1）；φ为孔隙度（-）；ρ_w为水相密度（kg∙m^-3）.

水相的密度是盐分浓度的线性函数，有以下关系式：

{\rho }_{\mathrm{w}}={\rho }_{\mathrm{f}}+\frac{\partial {\rho }_{\mathrm{w}}}{\partial C_{\mathrm{1}}}{C}_{\mathrm{1}},(2)

式（2）中：ρ_w为水相密度（kg∙m^-3）；ρ_f为淡水密度（kg∙m^-3）.

1.2.2　NAPL相氯苯入渗和迁移模型

NAPL相氯苯在含水层中的入渗和迁移遵循多相流理论.模型假设孔隙空间由水相和NAPL相共同占据，且两相不互溶、不可压缩，并处于等温稳定状态.在饱和多孔介质中，水相和NAPL相流体流速由广义达西定律和质量守恒方程控制（Feo et al.， 2023）：

\left\{\begin{array}{l}{\mathit{v}}_{\mathrm{w}}=-\frac{{\mathit{k}}_{\mathit{i}\mathit{j}}{k}_{\mathrm{r},\mathrm{w}}}{{\mu }_{\mathrm{w}}}\left(\nabla P_{\mathrm{w}}-{\rho }_{\mathrm{w}}\mathit{g}\right)\\ {\mathit{v}}_{\mathrm{n}}=-\frac{{\mathit{k}}_{\mathit{i}\mathit{j}}{k}_{\mathrm{r},\mathrm{n}}}{{\mu }_{\mathrm{n}}}\left(\nabla P_{\mathrm{n}}-{\rho }_{\mathrm{n}}\mathit{g}\right)\end{array}\right.,(3)

\left\{\begin{array}{l}\phi \frac{\partial \left({\rho }_{\mathrm{w}}{S}_{\mathrm{w}}\right)}{\partial t}-\nabla \left({\rho }_{\mathrm{w}}{v}_{\mathrm{w}}\right)=\mathrm{0}\\ \phi \frac{\partial \left({\rho }_{\mathrm{n}}{S}_{\mathrm{n}}\right)}{\partial t}-\nabla \left({\rho }_{\mathrm{n}}{v}_{\mathrm{n}}\right)=\mathrm{0}\end{array}\right.,(4)

式（3）~（4）中：φ为孔隙度（-）；S_w为水相饱和度（-）；S_n为NAPL相饱和度（-）；k_ij 为固有渗透率张量（m²）；k_r，w为水相的相对渗透率（-）；k_r，n为NAPL相的相对渗透率（-）；ρ_w为水相的密度（kg∙m^-3）；ρ_n为NAPL相的密度（kg∙m^-3）；μ_w为水相的动力黏度（Pa∙s）；μ_n为NAPL相的动力黏度（Pa∙s）；P_w为水相的压力（Pa）；P_n为NAPL的压力（Pa）；g为重力加速度（m∙s^-2）.

水相和NAPL相饱和度满足以下关系：

S_{\mathrm{n}}+S_{\mathrm{w}}=\mathrm{1} .(5)

相压之间由毛细管压力连接，有如下关系：

P_{\mathrm{c}}\left(S_{\mathrm{w}}\right)=P_{\mathrm{n}}-P_{\mathrm{w}},(6)

式（6）中：P_c为毛细管压力（Pa）.

本研究中，Brooks-Corey模型用于将毛细管压力与相饱和度联系起来（Kamon et al.， 2003）：

\left\{\begin{array}{l}{P}_{\mathrm{c}}\left(S_{\mathrm{w}}\right)=P_{\mathrm{d}}{S}_{\mathrm{e}}^{-\frac{\mathrm{1}}{\lambda }},P_{\mathrm{c}}\ge P_{\mathrm{d}}\\ S_{\mathrm{e}}=\frac{S_{\mathrm{w}}-S_{\mathrm{r},\mathrm{w}}}{\mathrm{1}-S_{\mathrm{r},\mathrm{w}}-S_{\mathrm{r},\mathrm{n}}}\\ k_{\mathrm{r},\mathrm{w}}=S_{\mathrm{e}}^{\frac{\mathrm{2}+\mathrm{3}\lambda }{\lambda }}\\ k_{\mathrm{r},\mathrm{n}}={\left(\mathrm{1}-S_{\mathrm{e}}\right)}^{\mathrm{2}}\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{2}+\lambda }{\lambda }}\right)\end{array}\right.,(7)

式（7）中：S_e为有效饱和度（-）；S_r，w为水相的残余饱和度（-）；S_r，n为NAPL相的残余饱和度（-）；λ为孔隙大小分布指数（-）；P_d为入口毛细压力（Pa）.

1.2.3　溶解相氯苯迁移模型

溶解态氯苯在水相的运移遵循对流-弥散方程（Wu et al.， 2016）：

\frac{\partial \left(\phi {\rho }_{\mathrm{w}}{S}_{\mathrm{w}}{C}_{\mathrm{2}}\right)}{\partial t}=\nabla \cdot \left(\phi {\rho }_{\mathrm{w}}{S}_{\mathrm{w}}{\mathit{D}}_{\mathrm{2}}\nabla C_{\mathrm{2}}\right)+\nabla \cdot \left[{\rho }_{\mathrm{w}}{C}_{\mathrm{2}}\frac{k_{\mathrm{r},\mathrm{w}}\mathit{k}}{{\mu }_{\mathrm{w}}}\left(\nabla P_{\mathrm{w}}+{\rho }_{\mathrm{w}}\mathit{g}\right)\right]+Q_{\mathrm{S}}{C}_{\mathrm{S}},(8)

式（8）中： D₂ 为溶解氯苯的水动力弥散系数（m²∙s^-1）；C₂为溶解氯苯浓度（g∙L^-1）；Q_s为源（汇）单位体积流量（m³∙s^-1）；C_s为源（汇）流中溶质浓度（g∙L^-1）；φ为孔隙度（-）；S_w为水相饱和度（-）；k_r，w为水相的相对渗透率（-）；ρ_w为水相密度（kg∙m^-3）.

NAPL相氯苯与溶解相氯苯之间的质量交换通过动力学传质过程实现，而非瞬时平衡假设.上式中源项来自NAPL相氯苯溶解到水相的动力学过程，其溶解速率通常采用线性驱动力模型来量化（Prieto-Estrada et al.， 2023）：

J=k\left(C_{\mathrm{2}}^{eq}-C_{\mathrm{2}}\right),(9)

式（9）中：J为氯苯溶解速率（g∙m^-3∙s^-1）；C₂^eq为溶解氯苯平衡浓度（g∙m^-3）；C₂为溶解氯苯的浓度（g∙m^-3）；k为传质速率系数（s^-1）.

当C₂C₂^eq时，NAPL相氯苯持续溶解，随着C₂趋近于C₂^eq，溶解速率J逐渐降低.

1.3　数值模型

本研究基于COMSOL Multiphysics 6.2平台，耦合“两相流达西定律”、“多孔介质相传递”和“多孔介质中的稀物质传递”模块进行求解.模型采用全耦合瞬态求解策略，结合变阶数向后差分公式与自适应时间步长控制求解精度，非线性迭代采用阻尼牛顿法.第一阶段模拟区域被离散为23 437个四边形网格.第二、三阶段模拟区域被离散为7 500个三角形网格，在Z=5 m处添加层数为20的边界层，以便更精细化处理DNAPL池的迁移过程.

模拟情景及模型参数如表1所示（Parker and Park， 2004； Page et al.， 2007； Saenton and Illangasekare， 2007）.

1.4　评价指标

1.4.1　NAPL相氯苯

通过绘制并分析注入点垂向剖面（X=300 m）及隔水底板（Z=0 m）沿线NAPL相饱和度S_n的时空变化曲线，定量评估DNAPL的垂向入渗深度、在隔水底板上的横向展布范围以及再分布过程中的收缩或聚集趋势.

在多孔介质中，当孔隙被水相、NAPL相等两种或两种以上不混溶的流体共同占据时，某一相的饱和度被定义为该相所占据的孔隙体积与孔隙总体积的比值.对于NAPL相，其饱和度S_n的数学定义为：

S_{\mathrm{n}}=\frac{V_{\mathrm{n}}}{V_{\mathrm{v}}},(10)

式（10）中：S_n为NAPL相饱和度（-）；V_v为典型单元体的孔隙体积（m³）；V_n为NAPL相氯苯占据的体积（m³）.

1.4.2　溶解相氯苯　(1) 溶解污染羽质心.在二维模拟域中，溶解相氯苯污染羽在某一时刻t的质心坐标（X_c （t），Z_c （t））由以下公式定义（Cheng et al.， 2024）：

X_c\left(t\right)=\frac{{\int }_{\mathrm{\Omega }}x\cdot C_{\mathrm{2}}\left(x,z,t\right)\mathrm{d}\mathrm{\Omega }}{{\int }_{\mathrm{\Omega }}{C}_{\mathrm{2}}\left(x,z,t\right)\mathrm{d}\mathrm{\Omega }},(11)

Z_c\left(t\right)=\frac{{\int }_{\mathrm{\Omega }}z\cdot C_{\mathrm{2}}\left(x,z,t\right)\mathrm{d}\mathrm{\Omega }}{{\int }_{\mathrm{\Omega }}{C}_{\mathrm{2}}\left(x,z,t\right)\mathrm{d}\mathrm{\Omega }},(12)

式（11）~（12）中，C₂（x，z，t）表示t时刻位置（x，z）处的溶解相氯苯浓度（kg∙m^-3）；Ω代表整个二维模拟域（-）；分母代表t时刻整个含水层中溶解相氯苯的总质量，起到归一化的作用；分子为浓度关于X坐标和Z坐标的一阶矩.

(2) 溶解氯苯质量空间分配.分别计算溶解相氯苯在整个含水层、淡水区、咸淡水过渡带内、咸水楔的累积质量随时间的变化.在二维模拟域中，某一特定区域Ω内在时刻t的溶解相氯苯累积质量M_Ω（t）计算公式为：

M_{\mathrm{\Omega }}\left(t\right)={\iint }_{\mathrm{\Omega }}^{}\phi \cdot S_{\mathrm{w}}\left(x,z,t\right)\cdot C_{\mathrm{2}}\left(x,z,t\right)\mathrm{d}x\mathrm{d}z，(13)

式（13）中，M_Ω（t）为区域Ω内在时刻t的溶解相氯苯累积质量，在二维模型代表单位厚度含水层的质量（kg）；S_w（x，z，t）为t时刻位置（x，z）处的水相饱和度（-）；C₂（x，z，t）为时刻t位置（x，z）的溶解相氯苯浓度（kg∙m^-3）.

依据盐度阈值划定积分区域，包括淡水区、咸淡水过渡带内、咸水楔.

海水盐度的10%和90%之间的区域被定义为咸淡水过渡带（Abarca and Clement， 2009； Lu and Luo， 2010）：

{\mathrm{\Omega }}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=\left\{\left(x,z\right)\left|\mathrm{0.1}{C}_{\mathrm{s}}{C}_{\mathrm{1}}\mathrm{0.9}{C}_{\mathrm{s}}\right.\right\}.(14)

在本研究中，定义盐度大于0.9C_s的区域为高盐度区（咸水楔）：

{\mathrm{\Omega }}_{\mathrm{s}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}}=\left\{\left(x,z\right)\left|C_{\mathrm{1}}\mathrm{0.9}{C}_s\right.\right\} .(15)

相应地，淡水区为含水层全部区域中除咸淡水过渡带和咸水楔以外的区域：

{\mathrm{\Omega }}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}}={\mathrm{\Omega }}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}-{\mathrm{\Omega }}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}-{\mathrm{\Omega }}_{\mathrm{s}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}} .(16)

(3) 溶解氯苯的边界通量.通过左边界的溶解相氯苯总输出通量，由对流和弥散两部分组成，是评估污染物向海洋环境释放风险的关键指标.计算公式为：

J={\int }_{{\Gamma }_S}[{\mathit{q}}_{\mathrm{w}}{C}_{\mathrm{2}}-(\phi S_{\mathrm{w}}\cdot D\cdot \nabla C_{\mathrm{2}})]\cdot \mathit{n}d\mathrm{\Gamma },(17)

式（17）中： q_w 为水相达西流速矢量（m∙s^-1）；C₂为溶解相氯苯浓度（g∙m^-3）；φ为孔隙度（-）；S_w为水相饱和度（-）； n 为边界的外法向单位矢量（-）.

1.5　敏感性分析

为评估关键水文地质参数不确定性对模型预测结果的影响，本研究对渗透系数K进行单因素敏感性分析.K值取值序列为1 m/d、2 m/d、3 m/d、4 m/d、5 m/d，其中5 m/d为基准模型K值.在基准模型基础上，仅系统改变K的取值，保持模型其他所有参数、初始条件及边界条件不变，依次运行完整的3个模拟阶段，通过对比分析第2.4节定义的关键评价指标对K变化的响应，量化其敏感性.

2 结果与讨论

2.1　DNAPL迁移与分布

2.1.1　NAPL相氯苯

DNAPL在含水层中的迁移受毛细捕获、粘性指进及重力驱动等机制调控.为揭示地下水-DNAPL两相流系统的动态演化机制，图2展示了泄漏阶段NAPL相氯苯饱和度S_n的时空演化规律.泄漏初期，NAPL相克服毛细阻力的同时在重力主导下呈指进状垂向下渗，其横向扩展相对滞后（图2a₁~2a₃）.泄漏3年后，入渗前缘抵达隔水底板（图2a₄），随后污染物转为横向铺展，并持续至泄漏结束（图2a₆）.最终，在含水层底部形成DNAPL污染池，而上部入渗路径中则残留有大量受毛细束缚的残余相.

为进一步解析S_n的时空分布特征，本研究选取泄漏点垂向剖面（X=300 m）与隔水底板界面（Z=0 m）进行分析.图2b显示，NAPL相在垂向剖面上表现为重力与毛细力协同控制的非混相驱替.注入点处维持约0.3的较高饱和度，驱动入渗锋面在泄漏0.5、1、2 a后距离隔水底板的高度分别为20 m、12 m及7 m，并于第3年时触及隔水底板.此后，DNAPL在底板运移受阻并逐渐积聚，形成厚度约5 m的污染池，其饱和度呈垂向正梯度分布，峰值达0.36，至泄漏结束时进一步升至0.40.

图2c揭示了隔水底板处的横向扩展过程.泄漏3年时，底板处监测到初始积聚（饱和度0.03）后，随后DNAPL在压力梯度驱动下向两侧铺展.然而，由于区域水力梯度由右侧淡水边界指向左侧海洋边界，加之咸水楔密度效应增强了对流拖曳，DNAPL向左侧海洋方向的扩展距离显著大于向右侧内陆方向，呈现非对称展布特征.随着持续注入，底板S_n峰值由0.29逐步升高至0.33及0.37，泄漏终止时，底板S_n达到0.4，极值中心始终稳定于注入点下方.

图3展示了再分布阶段S_n的时空演化规律.再分布初期（图3a₁~3a₃），注入停止导致垂向入渗通道残余饱和度显著降低.与此同时，在隔水底板处，DNAPL在浓度梯度与重力分量的协同驱动下持续向外扩展，但其向右侧的展布受到抑制，池体中心厚度也相应减薄.进入再分布后期（图3a₄~3a₆），污染池体的形态演化受毛细驱动与溶解作用的双重控制：一方面，污染物沿优势通道形成丝状残余路径；另一方面，相界面面积的增大会强化溶解通量导致池体边缘区域饱和度急剧衰减，呈现出边缘模糊化的特征.

为系统解析再分布阶段的相态演化机制，本研究进一步分析了垂向剖面（X=300 m）与底板界面（Z=0 m）的S_n.图3b显示，注入终止后，受毛细束缚下，垂向通道发生相重组.注入点饱和度由0.3迅速回落至残余饱和度水平.同时，底部DNAPL污染池在重力和毛细力耦合作用下持续收缩，其主体饱和度由0.4逐步衰减至0.23.图3c揭示了隔水底板处受背景水力梯度控制的非对称展布特征.经过100年的演化，下左侧锋面在压力梯度与重力的协同驱动下向前推进至180 m处，而右侧受逆向势能梯度抑制仅扩展至360 m.伴随相质量的重新分配及溶解作用的持续消耗，DNAPL池饱和度峰值同步由0.40降至0.23.

2.1.2　溶解相氯苯

溶解相氯苯浓度C₂的时空演化受NAPL自由相分布与地下水流场的共同控制.泄漏0.5年（图4a₁），溶解作用主要发生在垂向指进通道的相界面处，溶解相污染羽在含水层流场驱动下呈局部高浓度分布.随着泄漏持续进行，溶解羽随NAPL锋面向下扩展，但垂向迁移显著滞后于非水相前缘.泄漏至第3年时，尽管DNAPL已到达隔水底板，溶解相仍主要分布于底板以上区域，直至第4年迁移至底部边界（图4a₂~4a₃）.泄漏终止时（图4a₄），系统最终演变为以底板DNAPL池为核心、受水动力条件严格控制的稳定污染羽.再分布阶段初期，污染羽主要通过界面溶解在源区周边形成高浓度污染羽.之后，污染羽在对流-弥散驱动下整体向左侧迁移.约15年后（图4b₁），其前缘抵达咸淡水过渡带，并随之沿咸水楔界面发生斜向爬升直至到达地下水排泄口后迅速向海排泄（图4b₂）.

为明确海水入侵对DNAPL运移的影响，本文进一步对比分析了有/无海水入侵下污染羽质心的时空演化特征（图5a，5b）.在无海水入侵条件下，质心沿含水层中下部向海边界稳定运移，横坐标由300 m到150.56 m处，纵坐标30 m到10.12 m处，呈单调递减趋势.相比之下，受海水入侵影响下，质心横向迁移速率急剧减缓最终稳定于164.86 m处，垂向位置则停滞在约11.5 m处.

基于质心演化揭示的运移规律，区域质量分配特征（图5c）进一步证实海水入侵对污染物空间格局的调控作用.在泄漏初期（0~15年），对流-弥散主导下的高速运移促使系统溶解相总质量由1.492 kg激增至9.971 kg，且淡水区占比维持在99%以上.随着污染羽进入咸淡水过渡带（15~30年）输运与积聚，淡水区质量占比降至94.5%，而咸淡水过渡带累积质量由0.024 kg跃升至0.897 kg.到第30~105年，系统达到平衡，溶解相总质量渐增至18.561 kg，咸淡水过渡带滞留量达1.238 kg.

通过描述溶解相的入海通量，本文进一步阐明溶解相DNAPL的迁移特征（图5d）.无海水入侵时，边界通量于第16年检出后迅速攀升，至30年达1.0 kg/a，并缓慢上升至最终的1.4 kg/a.海水入侵背景下的通量演化规律与上述情景类似.排海通量于13年时初始响应，随后快速增长至0.32 kg/a.该阶段与咸淡水过渡带质量的快速累积呈显著同步性，证实咸水楔界面在截留污染物的同时，通过密度驱动效应构建了新的优势排泄路径.40~105年，通量增速放缓并终值趋于0.38 kg/a.相比于淡水背景，海水入侵情景下排海通量输出峰值削减62%.变密度流在咸水楔界面诱导形成的局部密度环流迫使溶解相沿界面斜向爬升，延长其滞留时间.同时，溶解相在过渡带中的大量积聚使通量呈现缓慢释放的特征.

2.2　渗透系数K对DNAPL迁移的影响

2.2.1　NAPL相氯苯迁移

渗透系数K对NAPL相氯苯迁移分布的影响如图6所示.随K值由1 m/d增至5 m/d，DNAPL污染池横向展布范围扩展了52.4%；饱和度峰值则由0.277衰减至0.240，且峰值中心向海洋边界偏移27.49 m（图6a）.同一时间内，NAPL相前缘位置随K增大而向海洋方向靠近；至105年，K=1 m/d和K=5 m/d时，前缘位置分别推进至221.79 m与179.33 m（图6b）.这是由于低渗透性介质中显著的毛细阻滞效应限制横向扩展，而高渗透性条件下重力驱动逐渐占据主导，强化多孔介质中的重力分异与相态重组.

2.2.2　溶解相氯苯迁移

渗透系数对溶解相氯苯迁移的影响如图7所示.随K值由1 m/d增至5 m/d，排海通量响应时间从30年提早至13年，进入准稳态的历时由65年缩短至20年（图7a），系统内溶解相总质量降幅达54.2%（图7b）.高渗透性介质一方面通过扩展自由相分布增大溶解界面，另一方面强化对流输运显著提升溶解氯苯向海洋边界的迁移.

较高的渗透系数增强密度驱动的对流效应，促使溶解氯苯质量分布向咸水区转移.淡水区质量占比从97.9%下降至93.0%，而咸淡水过渡带与咸水楔的质量占比分别提升至6.73%和0.27%（图7c，7d）.然而，咸淡水过渡带中的污染物累积呈现非单调演化趋势.K=3 m/d时，自由相分布、密度流与平流输运达到最优耦合状态，促进咸淡水过渡带中的质量积累.在更高K值条件下，增强的对流冲刷效应主导运移过程，抑制污染物在咸淡水过渡带内的滞留.

3 结论及展望

本研究通过建立并求解变密度流-多相流-溶质运移耦合数值模型，系统揭示海水入侵形成的复杂水动力与化学场对DNAPL在滨海含水层中迁移、转化及长期归宿的影响机制.主要结论如下.

（1）NAPL相氯苯在泄漏阶段受重力主导快速垂向入渗，而在再分布期转变为受毛细力与重力协同控制的横向非对称铺展.

（2）盐度梯度诱发的密度驱动对流根本性地改变溶解相氯苯的运移路径.其迫使溶解相沿咸水楔界面斜向爬升并向海洋排泄，显著降低排泄通量峰值，致使咸淡水过渡带演变为长期次生污染源.

（3）介质渗透系数K的增大显著加速NAPL相氯苯的横向迁移与溶解相输出响应，同时促使溶解氯苯质量分布向咸水区转移.

本研究作为机理探索性工作，采用含水层均质、各向同性的理想假设，旨在剥离复杂因素、聚焦揭示变密度流场对DNAPL多相运移的核心调控机制，同时忽略了温度变化、生物降解及吸附解吸等地球化学反应.在敏感性分析中，仅独立变化了渗透系数，而实际中Brooks-Corey模型中的入口毛细压力、孔隙分布指数与K存在内在的物理关联.针对模型均质各向同性假设及参数解耦分析的局限性，未来研究将进一步聚焦于非均质结构的精细刻画与复杂生物地球化学过程的耦合，并开展多参数全局敏感性分析，以系统量化参数交互对模型预测不确定性的影响.

参考文献

[1] Abarca, E., Clement, T. P., 2009. A Novel Approach for Characterizing the Mixing Zone of a Saltwater Wedge. Geophysical Research Letters, 36(6):L06601. https://doi.org/10.1029/2008gl036995

[2] Cheng, Z., Wu, Y. H., Wu, M., et al., 2024. A Higher Mass Flux from the Effected DNAPL Migration and Distribution Patterns by Brine Flow in Saturated Porous Media. Chemical Engineering Science, 295: 120161. https://doi.org/10.1016/j.ces.2024.120161

[3] Cheng, Z., Wu, J.C., Xu, H.X., et al., 2014. Study on the Migration of DNAPL under the Influence of Lenses and Surfactants. China Environmental Science, 34(11): 2888-2896 (in Chinese with English abstract).

[4] Deng, Y.P., Zhang, Y., Shi, X. Q., et al., 2015. Migration of Dense Non-Aqueous Phase Liquids in Heterogeneous Fractured Media. Advances in Water Science, 26(5): 722-730 (in Chinese with English abstract).

[5] Deng, Y.P., Zheng, F., Shi, X.Q., et al., 2016. Research Progress on the Migration Behavior of DNAPLs in Porous Media. Journal of Nanjing University (Natural Sciences), 52(3): 409-420 (in Chinese with English abstract).

[6] Fang, Y. H., Zheng, T. Y., Wang, H., et al., 2021. Experimental and Numerical Evidence on the Influence of Tidal Activity on the Effectiveness of Subsurface Dams. Journal of Hydrology, 603: 127149. https://doi.org/10.1016/j.jhydrol.2021.127149

[7] Feo, A., Celico, F., Zanini, A., 2023. Migration of DNAPL in Saturated Porous Media: Validation of High-Resolution Shock-Capturing Numerical Simulations through a Sandbox Experiment. Water, 15(8): 1471. https://doi.org/10.3390/w15081471

[8] Gao, S. B., Zheng, T. Y., Zheng, X. L., et al., 2024. Influence of Pumping Saline Groundwater on Land-Sourced Solute Transport in Coastal Aquifers. Physics of Fluids, 36(12): 126608. https://doi.org/10.1063/5.0242660

[9] Hu, S. K., Deng, Z. H., Liu, B. J., et al., 2024. Impact of Tidal Dynamics and Typhoon-Induced Inundation on Saltwater Intrusion in Coastal Farms. Science of the Total Environment, 915: 170109. https://doi.org/10.1016/j.scitotenv.2024.170109

[10] Illy, V. D., Cohen, G. J. V., Verardo, E., et al., 2022. Using 1, 1, 1-Trichloroethane Degradation Data to Understand NAPL Dissolution and Solute Transport at Real Sites. Journal of Contaminant Hydrology, 245: 103934. https://doi.org/10.1016/j.jconhyd.2021.103934

[11] Ju, M.D., Li, X., Yin, H.L., 2025. Simulation Study on Migration and Transformation of Dense Non-aqueous Phase Liquids in Groundwater Coupled with Multiphase Flow and Solute Transport. Journal of Hydrodynamics (Series A), 40(1): 149-156 (in Chinese with English abstract).

[12] Kamon, M., Endo, K., Katsumi, T., 2003. Measuring the k-S-P Relations on DNAPLs Migration. Engineering Geology, 70(3-4): 351-363. https://doi.org/10.1016/S0013-7952(03)00103-0

[13] Kueper, B. H., Redman, D., Starr, R. C., et al., 1993. A Field Experiment to Study the Behavior of Tetrachloroethylene below the Water Table: Spatial Distribution of Residual and Pooled DNAPL. Groundwater, 31(5): 756-766. https://doi.org/10.1111/j.1745-6584.1993.tb00848.x

[14] Li, P.H., Liu, H.L., Hou,S.G., 2021. Study on DNAPL Contamination Process Using Electrical Resistivity Imaging and Image Method. Progress in Geophysics, 36(5): 2185-2190 (in Chinese with English abstract).

[15] Lu, C. H., Luo, J., 2010. Dynamics of Freshwater-Seawater Mixing Zone Development in Dual-Domain Formations. Water Resources Research, 46(11): 2010WR009344. https://doi.org/10.1029/2010wr009344

[16] Michael, H. A., Post, V. E. A., Wilson, A. M., et al., 2017. Science, Society, and the Coastal Groundwater Squeeze. Water Resources Research, 53(4): 2610-2617. https://doi.org/10.1002/2017wr020851

[17] Page, J. W. E., Soga, K., Illangasekare, T., 2007. The Significance of Heterogeneity on Mass Flux from DNAPL Source Zones: An Experimental Investigation. Journal of Contaminant Hydrology, 94(3-4): 215-234. https://doi.org/10.1016/j.jconhyd.2007.06.004

[18] Parker, J. C., Park, E., 2004. Modeling Field-Scale Dense Nonaqueous Phase Liquid Dissolution Kinetics in Heterogeneous Aquifers. Water Resources Research, 40(5): 2003WR002807. https://doi.org/10.1029/2003wr002807

[19] Praseeja, A. V., Sajikumar, N., 2019. A Review on the Study of Immiscible Fluid Flow in Unsaturated Porous Media: Modeling and Remediation. Journal of Porous Media, 22(8): 889-922. https://doi.org/10.1615/jpormedia.2019024580

[20] Prieto-Estrada, A. E., Widdowson, M. A., Stewart, L. D., 2023. Numerical Modeling and Data-Worth Analysis for Characterizing the Architecture and Dissolution Rates of a Multicomponent DNAPL Source. Water Resources Research, 59(5): e2022WR034351. https://doi.org/10.1029/2022wr034351

[21] Prommer, H., Barry, D. A., Davis, G. B., 2002. Modelling of Physical and Reactive Processes during Biodegradation of a Hydrocarbon Plume under Transient Groundwater Flow Conditions. Journal of Contaminant Hydrology, 59(1-2): 113-131. https://doi.org/10.1016/S0169-7722(02)00078-5

[22] Saenton, S., Illangasekare, T. H., 2007. Upscaling of Mass Transfer Rate Coefficient for the Numerical Simulation of Dense Nonaqueous Phase Liquid Dissolution in Heterogeneous Aquifers. Water Resources Research, 43(2): 2005WR004274. https://doi.org/10.1029/2005wr004274

[23] Shi, J. X., Chen, X. H., Ye, B., et al., 2023. A Comparative Study of DNAPL Migration and Transformation in Confined and Unconfined Groundwater Systems. Water Research, 245: 120649. https://doi.org/10.1016/j.watres.2023.120649

[24] Wu, B., Li, H. Y., Du, X. M., et al., 2016. Correlation between DNAPL Distribution Area and Dissolved Concentration in Surfactant Enhanced Aquifer Remediation Effluent: A Two-Dimensional Flow Cell Study. Chemosphere, 144: 2142-2149. https://doi.org/10.1016/j.chemosphere.2015.11.005

[25] Wu, H. O., Tu, X. J., Chen, X. H., et al., 2023. A Framework for Water Supply Regulation in Coastal Areas by Avoiding Saltwater Withdrawal Considering Upstream Streamflow Distribution. Science of the Total Environment, 905: 167181. https://doi.org/10.1016/j.scitotenv.2023.167181

[26] Wu, M., Cheng, Z., Wu, J. F., et al., 2017a. Effects of Inner Heterogeneity on Long-Term DNAPL Migration in Porous Media. Environmental Earth Sciences, 76(23): 791. https://doi.org/10.1007/s12665-017-7141-4

[27] Wu, M., Wu, J. F., Wu, J. C., 2017b. Simulation of DNAPL Migration in Heterogeneous Translucent Porous Media Based on Estimation of Representative Elementary Volume. Journal of Hydrology, 553: 276-288. https://doi.org/10.1016/j.jhydrol.2017.08.005

[28] Wu, M., Yang, Y., Lin, J., et al., 2020. The Co-Effect of Heterogeneity and Solute Concentration on Representative Elementary Volume of DNAPL in Groundwater. Journal of Hydrology, 585: 124795. https://doi.org/10.1016/j.jhydrol.2020.124795

[29] Yang, B., Li, H.Y., Wu, B., et al., 2013. Study on the Fingering Front Morphology of Four NAPLs in Two-Dimensional Sandbox. Environmental Science, 34(4): 1545-1552 (in Chinese with English abstract).

[30] Yang, L. R., Wang, X. Y., Mendoza-Sanchez, I., et al., 2018. Modeling the Influence of Coupled Mass Transfer Processes on Mass Flux Downgradient of Heterogeneous DNAPL Source Zones. Journal of Contaminant Hydrology, 211: 1-14. https://doi.org/10.1016/j.jconhyd.2018.02.003

[31] Yao, J. W., Zhang, H., 2023. Comparing Darcy’s Law and the Brinkman Equation for Numerical Simulations of Saltwater Intrusion. Sustainability, 15(18): 13287. https://doi.org/10.3390/su151813287

[32] Yu, L., Zheng, T. Y., Yuan, R. Y., et al., 2022. APCS-MLR Model: A Convenient and Fast Method for Quantitative Identification of Nitrate Pollution Sources in Groundwater. Journal of Environmental Management, 314: 115101. https://doi.org/10.1016/j.jenvman.2022.115101

[33] Zhang, J. X., Lu, C. H., Zhang, C. M., 2024. Dense Contaminants Mixing into the Saltwater Wedge in Coastal Aquifers: Laboratory and Numerical Investigations. Water Resources Research, 60(7): e2024WR037452. https://doi.org/10.1029/2024wr037452

[34] Zhang, W., Shi, X.Q., Wu, J.F., et al., 2013. Influence of Permeability Spatial Variability on the Migration of Dense Non-Aqueous Phase Liquids. Geological Journal of China Universities, 19(4): 677-682 (in Chinese with English abstract).

[35] Zheng, T. Y., Gao, S. B., Liu, T., et al., 2023. Dynamic Influence of Land Reclamation on the Nitrate Contamination and Saltwater Redistribution. Science of the Total Environment, 883: 163605. https://doi.org/10.1016/j.scitotenv.2023.163605

基金资助

国家自然科学基金优秀青年科学基金项目(42422207)