沿海地区经济发展导致需水量迅猛增长，为保持经济发展速度，弥补地区用水量的不足，地下水被过度开采。同时全球变暖导致海平面上升^[1-3]。在这两方面因素的双重作用下，滨海地区含水层中天然的咸淡水动态平衡系统被破坏，咸淡水界面向内陆地下淡水运移，引发了海水入侵。海水入侵严重威胁滨海地区的地下淡水资源^[4-6]。

海水入侵是一个复杂的环境变化过程，受多种自然条件和人类活动的影响。其核心影响因素是过量开采地下水造成水动力平衡的破坏^[7-10]。海水入侵的影响因素包括地质、构造、岩性、含水层渗透性、含水层补给条件、含水层在海底方向上的延伸状况、大气降水等^{[4, 11-12]}。建立海水入侵的数学模型是研究海水入侵机理和预测未来的有效手段。国外对海水入侵关注较早，也开发了一系列数学模型，用以模拟预测海水入侵的范围、持续时间等^[13-17]。目前海水入侵的研究已经从单一化研究向综合化研究方向发展，从简单的定性分析向定量研究方向趋近，从概化含水层结构、地质条件建立概念模型，最终用数学的思想解决问题^[18-20]。尤其在海水入侵方面的基础理论认识、数学模型的建立和完善、预测和综合防治等方面已逐渐成熟。但现有海水入侵模型均基于饱和含水层的模拟（如SEAWAT）^[21-24]，忽略了地表水和地下水的联动作用。在气候变化和人类活动影响下，地表水率先响应，进而对地下水产生影响^[25]。因而，地表水和地下水耦合模拟可全面准确地描述海水入侵过程，及在全球变化条件下，入侵过程对水资源的影响。在海水入侵模拟理论及方法研究的基础上，本文提出了海水入侵模拟及预测模型——变密度水流与溶质运移模型（VFT3D）。本文系统介绍了VFT3D模型，并提供了应用案例分析。

1.   VFT3D模型

1.1   模型理论

VFT3D模型地下水溶质迁移部分是多组分反应迁移模型PHT3D（A reactive multicomponent transport model for saturated porous media）的变密度版本。PHT3D在恒定密度系统中通过溶质迁移模型MT3DMS（A modular three-dimensional multispecies transport model）和地球化学反应模型PHREEQC-2求解对流-弥散-反应迁移问题^[26]。VFT3D将MT3DMS和HEIFLOW（Hydrological-Ecological Integrated watershed-scale FLOW model）结合求解非反应性多组分溶质在地表水-地下水耦合系统变密度条件下的迁移过程。

1.1.1   变密度地下水流中非反应性多组分溶质迁移

经典的恒密度水流模型采用体积守恒建立水流方程，但对于变密度地下水模型，体积守恒不再适用。因此，采用质量守恒建立变密度地下水流方程，并将地下水的密度定义为地下水中溶解物质浓度的函数，忽略温度对地下水密度的影响。Guo等^[27]推导出包含等价淡水水头的变密度地下水流动方程：

式中：$K$——渗透率；

${S}_{\mathrm{f}}$——以等价淡水水头表示的储水率；

$\rho$——地下水的密度；

$\theta$——有效孔隙度；

$C$——溶解物质的浓度；

$q_{\rm{s}}$——单位体积含水层补给或者蒸发的体积流量；

$\stackrel{-}{\rho }$——源、汇项中的溶解物质的密度；

ρ_f——淡水密度；

h_f——等价淡水水头。

通过地下水密度耦合了密度驱动的地下水流动方程和溶质运移方程。在无反应的模拟中，线性状态方程（流体密度为浓度的线性函数）为：

用Richards方程描述非饱和带-饱和带联立的流动过程^[28]：

式中：C(h)——容水度；

h——水头。

变密度条件下的水流和非反应性溶质运移模拟模型通过耦合流模型和溶质运移模型MT3DMS建立，称为MT3DMS2。

由于地下水流动方程式中右边有描述含水层中浓度随时间改变的偏微分项，所以地下水流流动会引起浓度的重新分配，进而影响流场。因此，地下水流动和溶质运移是一个耦合的过程，两个控制方程必须耦合求解。多孔介质的溶质运移方程包括水动力弥散项、地下水对流项、源汇项以及反应项：

式中：${\boldsymbol{D}}$——水动力弥散系数张量；

${C}_{{\rm{s}}}$——通过补给进入或者蒸发离开的溶质的浓度；

q——达西流速；

$R$——反应项。

前3项由MT3DMS2代码求解。在反应模拟中，第4项由耦合代码中基于PHREEQC-2的部分计算。采用空间中心加权的隐式有限差分格式求解对流项。

1.1.2   变密度多组分溶质反应性迁移

为了考虑反应过程，使用了多组分反应迁移模型的变密度版本。模型通过连续算子分裂技术将地球化学反应模型PHREEQC-2与变密度条件下的溶质运移模型MT3DMS2耦合。与早期版本的反应迁移模型相比，PHREEQC-2通过计算H和O元素的质量平衡计算pH和pE，因此，溶质运移模拟中包括总氢、总氧。

电荷平衡（或更准确地说是电荷不平衡）也是一种运移属性。但是，在关键情况下，如果模拟中应用的所有水溶液（无论是初始条件还是边界条件）都完全电荷平衡，则可以获得更可靠的模型结果。

第n个水相（流动相）化学组分的反应迁移方程

（MT3DMS求解）：

对于非流动相，例如矿物：

式中：${v}_{\alpha }$——$\alpha$方向上的孔隙流速；

${D}_{\alpha \beta }$——水动力弥散系数；

${C}_{n}^{{\rm{s}}}$——源/汇通量浓度；

${r}_{{\rm{reac}},n}$——由化学反应产生的源/汇速率；

${C}_{n}$——第n个水相（流动相）化学组分的浓度。

式中：${c}_{n}$——第n个水相（流动相）化学组分的摩尔浓 度（未络合）；

${n}_{{\rm{s}}}$——与水相组分络合的溶解形式物质的数量；

${Y}_{m}^{{\rm{s}}}$——第m种络合物中水相组分的化学计量数；

${s}_{m}$——第m种络合物的摩尔浓度。

2.   模型结构

VFT3D耦合HEIFLOW、MT3DMS、PHREEQC-2模型功能，实现变密度条件下多组分溶质反应性迁移。

HEIFLOW是一个基于物理过程的分布式生态水文模型^[29]，能够模拟完整的水循环过程、植被生态过程和水资源利用过程。其中水循环过程可以进一步分为地表部分（蒸散发、入渗、坡面流等）、地下部分（非饱和带水分运移、饱和带水流等）和水体部分（河道汇流、湖泊动态变化等）；植被生态过程主要包括植物生长和水分利用的动态，由生物量、LAI等变量的动态变化体现；水资源利用过程主要是指流域内用于农业、工业生活的地表引水和地下水开采过程。HEIFLOW通过松散耦合方式与MT3DMS2连接，为VFT3D提供水流运动信息，实现了完整水循环过程的海水入侵的模拟。

VFT3D包括一个主程序和一系列相对独立的子程序包，子程序包包括多个模块和子程序。VFT3D的基本程序结构见图1。整个模拟过程可分为一系列应力期（stress periods）。在每一个应力期内，所有的外应力（stress），如抽水量、蒸发量等保持不变。每个应力期又可再分为若干个时间段（time steps）。通过对有限差分方程组的迭代求解，可以得到每个时间段结束时的溶质浓度。模拟包括4大循环：应力期循环、时间段循环、运移反应段循环以及迭代求解循环。

图  1  主程序流程图

Figure  1.  Main program flow chart

下载: 全尺寸图片 幻灯片

在反应迁移模型中，MT3DMS2和PHREEQC-2之间的耦合是通过几个新的子程序实现的，这些子程序支持2个程序之间的通信，并调用PHREEQC-2数据库计算由于反应产生的化学传输。MT3DMS2的代码作为驱动程序（主要程序），PHREEQC-2的接口程序包（the PHREEQC-2 interface package， PHC）作为MT3DMS2和PHREEQC-2之间的接口，其结构与MT3DMS2中其他程序包的结构相同（例如对流、弥散、反应程序包），因此只需对MT3DMS2的源代码做很小的修改。MT3DMS2采用Fortran 90，PHREEQC-2采用C语言，PHC程序包可以选择Fortran 90或C语言。

反应迁移模型中虽然添加了新的子程序，但仍保留了MT3DMS2和PHREEQC-2的主要模型结构。若干子程序已经编写完成，在MT3DMS2主程序循环时可以被调用（图1）。子程序PHC1AL为MT3DMS2中PHC程序包使用的数组分配内存。子程序PHC1RP初始化几个变量并调用函数MCRP，该函数读取反应迁移模型中特定的输入并为PHREEQC-2准备输入数据，同时也调用PHREEQC-2的函数计算溶液的初始组分、相态组合、交换量和表面。子程序PHC1SS确保通过MT3DMS SSM软件包定义的任何溶液在PHREEQC-2中均进行组分计算，得到总氢和总氧。子程序PHC1RU调用PHREEQC-2计算时间步长内的反应，并确定在随后的反应步骤中单元格是否需要被激活。子程序PHC1BD计算化学反应转移对总质量平衡的贡献。最后，PHC1OT子程序控制输出（写为未格式化的二进制或者ACSII文件）。建立并运行VFT3D模型的工作流程见图2。

图  2  VFT3D模型的工作流程

Figure  2.  Workflow of VFT3D

下载: 全尺寸图片 幻灯片

3.   案例模拟

3.1   概念模型

本案例描述了发生在承压沿海含水层二维截面的海水入侵过程。其概念模型为：承压含水层区域长2 km，高1 km，宽1 m，上下顶板和底板为隔水边界，左边界为给定流量边界（Well），流量为10 m³/d，右边界为定水头边界。承压沿海含水层的二维截面初始以淡水饱和，右边界为定浓度边界，设置为淡水；左边界为给定溶质通量边界，设置为海水。海水以给定流量从左边入侵，向右边的淡水边界排泄。模型区域以及相关水动力学参数取值见图3。模型区域在空间上离散为100层、200列、1行，网格间距为10 m$\times$10 m，共20000个有效单元格，模型模拟期为60000 d。

图  3  模型区域及参数取值

Figure  3.  Model domain and parameters used

下载: 全尺寸图片 幻灯片

本案例考虑变密度多组分溶质反应性迁移过程，对流体密度产生影响的组分包括Cl⁻、Na⁺、SO₄^2-、Mg²⁺、Ca²⁺、K⁺、${\rm{HCO}}^{-}_3$。对于变密度过程来说，溶质浓度与流体密度相关线性状态方程的斜率($\partial \rho /\partial C$)是一个非常重要的参数，该参数可理解为每个组分浓度变化对应的密度变化。由于水溶液组分之间的离子相互作用，密度与多组分浓度之间的关系十分复杂。Monnin^[30]基于Pitzer的离子相互作用模型开发了VOPO程序。该程序能够联系密度与各组分浓度，计算得到的海水密度误差为0.00002 g/cm³。但当流体密度较高时，计算密度误差约为0.0004 g/cm³。本次研究用VOPO程序检验密度和浓度之间的经验线性关系。基于VOPO程序理论，密度和浓度之间的经验线性关系表述为：

式中：${\varepsilon }_{i}$——描述第i个组分的浓度对流体密度影响 的系数，即其是每个组分的密度-浓度耦 合项。

忽略密度对流体动力学黏度的影响，因为它主要受温度变化的影响，而不受密度变化影响^[27]。模拟过程中给定的咸水（海水）、淡水的各组分浓度值（包括pH值）以及描述各组分浓度对流体密度影响的系数${\varepsilon }_{i}$均来自Mao等^[31]的研究，${\varepsilon }_{i}$值利用VOPO计算程序确定，咸水（海水）、淡水的各组分浓度值以及对应的${\varepsilon }_{i}$值见表1。结果表明，由于离子相互作用的影响，密度对阴离子的浓度变化最为敏感。

表  1  咸水（海水）、淡水的各组分浓度值以及对应的${\varepsilon }_{i}$值

Table  1.  Concentrations of ions in sea water and fresh water and values for ${\varepsilon }_{i}$

	pH	Cl−	Na+	${\rm{SO}}^{2-}_4$	Mg2+	Ca2+	K+	${\rm{HCO}}^-_3$
淡水/(mol·L−1)	7.93	0.0296	0.0205	0.0067	0.0017	0.0092	0.0008	0.0003
海水/(mol·L−1)	7.95	0.5295	0.4521	0.0265	0.0483	0.0100	0.0092	0.0004
${\varepsilon}_{i}$值/(g·mol−1)		35.1	6	116.8	9.6	21.7	12.2	54.4
注：pH值无量纲。

下载: 导出CSV 

| 显示表格

国内外的多项研究表明阳离子吸附交换作用是海水入侵、咸淡水混合过程中不可忽视的一种地球化学作用。吴吉春等^[32-33]指出海水入侵过程中，咸淡水不是简单的机械混合，而是同时受到了水岩间阳离子交替吸附作用的影响，其中以Na⁺−Ca²⁺交换为主，其次为Mg²⁺−Ca²⁺交换，K⁺−Ca²⁺交换很少。因此，模拟过程中主要考虑的阳离子吸附交换作用及反应常数见表2。模拟区域介质的水平渗透系数为10 m/d，接近砂的渗透系数。阳离子交换容量（cation exchange capacity，CEC）受多种因素的影响，沉积物颗粒越细，CEC值越大。砂的粒径较黏土类偏大，且不考虑其它吸附剂（例如：有机胶体）等的影响，因此，阳离子交换容量初始值可以设置的偏低，为0.0001 mol/L。

表  2  阳离子交换作用及反应常数

Table  2.  Ion-exchange reactions in the column

阳离子交换作用反应公式	lg Ka
Na++X−=NaX	0
Mg2++2X−=MgX2	0.6
Ca2++2X−=CaX2	0.8
注：a反应常数来自于PHREEQC-2的数据库。

下载: 导出CSV 

| 显示表格

3.2   结果与讨论

3.2.1   地球化学反应对变密度流溶质运移过程的影响

SEAWAT和VFT3D模拟计算的地下水水位结果见图4（a）。两种模型计算所得的水位结果差别不明显，左下角海水入侵处的水位偏低，右边淡水定水头边界处的水位偏高，随着时间的推移，海水整体呈向下、向右流动。

图  4  SEAWAT和VFT3D模拟计算的水位结果对比和Cl⁻、Na⁺迁移过程

Figure  4.  Comparison of the hydraulic heads and migration of Cl⁻ and Na⁺ calculated with SEAWAT and VFT3D

下载: 全尺寸图片 幻灯片

Cl⁻化学性质较为稳定，化学反应对Cl⁻的迁移过程几乎没有影响，可以看作惰性示踪剂，其在含水层的迁移过程可以反映溶质的物理迁移特征和水动力弥散特征。SEAWAT和VFT3D模型模拟的Cl⁻迁移过程见图4（b）。Cl⁻向右侧迁移的过程不均匀，区域下方向右迁移更快，这与水位结果显示的地下水向右下流动的趋势相符。区域左侧侵入密度较大的海水，区域一开始以密度较低的淡水饱和且右侧边界为以淡水设置的定浓度边界，因此受到密度流的影响，区域上部Cl⁻向右侧迁移受阻，迁移十分缓慢，随着深度加深，迁移变快。60000 d时，尽管一半以上的区域Cl⁻浓度较高甚至接近海水中Cl⁻的浓度，但是区域右上角的浓度偏低，接近淡水中Cl⁻的浓度。Cl⁻迁移过程中，从高浓度向低浓度过渡带宽度偏窄，显示出较大的浓度梯度。${\rm{SO}}^{2-}_4$和${\rm{HCO}}^-_3$的迁移过程与Cl⁻基本一致，SEAWAT和VFT3D模拟计算的阴离子运移过程差别不大，推测阴离子在迁移过程中受地球化学反应的影响不大。

尽管变密度流多组分反应迁移模拟中考虑了Na⁺的阳离子交换作用，但是VFT3D模拟计算结果与SEAWAT模拟计算结果差别不大，见图4（c），推测Na⁺在迁移过程中，受阳离子吸附作用的影响不太明显。受密度流的影响，区域上部迁移受阻，Na⁺在区域下部向右迁移明显快于上部区域的迁移过程。到60 000 d时，一半以上的区域Na⁺浓度较高甚至接近海水中Na⁺的浓度，区域右上角的浓度偏低，接近淡水中Na⁺的浓度，过渡带偏窄，浓度梯度较高。K⁺和Mg²⁺表现出与Na⁺较为一致的迁移特征，且Na⁺、K⁺和Mg²⁺的迁移特征与上述水位结果显示的地下水流动的趋势以及阴离子的迁移特征一致。

在不考虑反应的SEAWAT模拟计算中，Ca²⁺迁移过程与其它组分的迁移过程基本一致（图5）。但是，VFT3D模拟计算的Ca²⁺迁移过程明显与其它组分的迁移特征不同（图6）。受密度流的影响Ca²⁺在区域下部向右迁移明显快于上部区域的迁移。到60 000 d时，左侧大半的区域Ca²⁺浓度接近海水中的浓度，区域右上角的浓度接近淡水中的浓度。与其它组分迁移过程中的过渡带宽度特征一致，Ca²⁺迁移过程中的过渡带偏窄。但是，其它组分在过渡带中从左向右浓度逐渐变低，从海水浓度过渡为淡水浓度，表征了咸淡水的机械混合作用过程。而Ca²⁺在过渡带中的浓度高于其他区域的Ca²⁺浓度，且高于海水中的浓度，表现出从左向右，海水Ca²⁺浓度向过渡带中心方向升高，之后又降低至淡水中Ca²⁺的浓度。随着时间的推移，过渡带中尤其是过渡带中心处Ca²⁺浓度变高。Ca²⁺浓度在过渡带中心的富集是海水入侵过程中阳离子吸附交换作用的结果。地下淡水中，土壤颗粒上负电荷的饱和主要是靠Ca²⁺，海（咸）水中土壤颗粒上的负电荷达到饱和则主要是靠Na⁺。当海（咸）水入侵地下淡水时，土壤颗粒表面释放了吸附的Ca²⁺，从而导致Ca²⁺浓度升高。因此，在模型区域海水向右入侵流动的过程中，区域左侧介质中的Ca²⁺被解吸交换，致使Ca²⁺浓度增加且高于海水中Ca²⁺浓度，并随着水的流动在过渡带中心富集，这与李梦娇^[34]的研究中海淡水驱替实验的结论一致。

图  5  SEAWAT模拟计算的Ca²⁺迁移过程

Figure  5.  Migration of Ca²⁺ from the calculation using SEAWAT

下载: 全尺寸图片 幻灯片

图  6  VFT3D模拟计算的Ca²⁺迁移过程

Figure  6.  Migration of Ca²⁺ from the calculation with VFT3D

下载: 全尺寸图片 幻灯片

3.2.2   阳离子交换容量对海水入侵过程中溶质迁移的影响

考虑含水介质为砂，阳离子交换容量值为0.0001 mol/L情况下，阳离子吸附交换作用对Ca²⁺的迁移过程产生明显影响。但是砂的CEC值偏低，对于有机胶体含量高且黏土矿物以蒙脱石、伊利石为主的细颗粒含水介质来说，其CEC值偏高，推测CEC值的变化对溶质运移过程产生影响，因此，将CEC值升高至0.01 mol/L，探讨不同CEC值条件下的阳离子吸附交换作用对溶质运移的影响。

高CEC值条件下，VFT3D模拟计算的水位结果见图7（a），其与低CEC值条件下的水位结果，图4（a），具有相似特征。地下水呈现向下、向右的流动趋势。高CEC值条件下，Cl⁻的迁移过程见图7（b），Cl⁻以及其它阴离子迁移特征与低CEC值条件下阴离子迁移特征，图4（b），一致。受变密度流的影响，较下部区域，上部区域阴离子迁移过程较为缓慢，迁移至60 000 d，左侧大部分区域被海水咸化，阴离子浓度较高接近海水浓度，仅右上角小部分区域维持淡水浓度，过渡带浓度从左到右逐渐降低，过渡带偏窄。随着CEC值的升高，阴离子迁移过程变化不大，表明其受阳离子吸附交换作用的影响不明显。

图  7  高CEC值条件下VFT3D模拟的水位结果和Cl^-迁移过程

Figure  7.  Spatial distribution of hydraulic heads and migration of Cl^- from the calculation using VFT3D with high concentration of CEC

下载: 全尺寸图片 幻灯片

受高CEC值影响的阳离子迁移过程见图8。与阴离子过程相似，阳离子迁移过程也表现出上部迁移慢于下部迁移的特点。但是，高CEC条件下阳离子迁移过程的过渡带表现出不同的特点。对于Na⁺和K⁺来说，高CEC值导致其过渡带宽度随着时间的推移逐渐变宽。对于Mg²⁺和Ca²⁺来说，受阳离子吸附交换作用的影响，过渡带中Mg²⁺和Ca²⁺浓度高于其他区域Mg²⁺和Ca²⁺的浓度，且过渡带中心的浓度最高，随着时间的推移，过渡带中Mg²⁺和Ca²⁺浓度升高，过渡带以及高Mg²⁺和Ca²⁺浓度区域逐渐变宽。在低CEC值条件下，仅有Ca²⁺在过渡带中富集，Mg²⁺并没有富集表现。由于CEC值变高，Mg²⁺受阳离子吸附交换作用的影响开始表现出来，且过渡带中心的Mg²⁺和Ca²⁺浓度也高于低CEC条件下过渡带中心Mg²⁺和Ca²⁺浓度。迁移过程中，K⁺迁移特征较为特别，不同于低CEC条件下K⁺与其它组分同时迁移的特点。高CEC条件下，随着时间的推移，K⁺在后期的迁移较其它组分来说逐渐变慢，该现象说明在海水向区域右侧入侵过程中，可能偏向于将介质中的Mg²⁺和Ca²⁺释放而吸附K⁺，从而导致Mg²⁺和Ca²⁺在其过渡带区域浓度升高、K⁺在对应区域浓度降低，表现出运移滞后的特点。当运移发生至60 000 d时，高浓度K⁺仅覆盖了左侧一半的区域，污染面积较其它组分偏小。

图  8  高CEC值条件下VFT3D模拟的Na⁺、K⁺、Mg²⁺和Ca²⁺迁移过程

Figure  8.  Migration of Na⁺, K⁺, Mg²⁺ and Ca²⁺ from the calculation using VFT3D with high concentration of CEC

下载: 全尺寸图片 幻灯片

4.   结论

本文详细描述了海水入侵变密度水流与溶质反应运移模型VFT3D及其应用。

（1）VFT3D模型求解反应性多组分溶质在变密度条件下的迁移过程，对复杂条件下咸-淡水过渡带进行准确描绘。同时，解决了SEAWAT等软件无法刻画具有强烈地表水-地下水交互作用的沿海区域地下水动态变化的海水入侵过程的难题，实现了地表水-地下水对海水入侵的协同控制模拟。

（2）案例中VFT3D 和SEAWAT结果对比发现，受阳离子吸附交换作用的海水入侵过程中，水位变化以及阴离子（Cl⁻、${\rm{SO}}^{2-}_4$、${\rm{HCO}}^-_3$）迁移过程受地球化学反应的影响不明显，但地球化学过程（阳离子吸附交换作用）对阳离子（Na⁺、K⁺、Mg²⁺和Ca²⁺）运移过程产生明显影响。高CEC值导致更多Mg²⁺和Ca²⁺在过渡带富集，K⁺运移滞后。地球化学作用引起过渡带中组分浓度发生变化，导致地下水化学类型发生转变，其对海水入侵过程的影响也是不容忽视的。

因此，考虑变密度和复杂反应过程，才能更加准确地描述输运现象，从而有利于地下咸水治理工程的科学实施。

VFT3D模型可应用于沿海地区水资源、水环境及水利工程的模拟和管理之中，有效提高沿海地区海水入侵的模拟精度、综合提升水资源管理水平、促进海岸带区域地下水资源利用及经济发展和生态需求的协调发展。