永兴岛淡水透镜体数值模拟：水资源评估和管理的重要工具

根据建立的淡水透镜体变密度地下水流和溶质运移三维数学模型，运用三维可视化软件Visual Modflow建立数值模型，对求解区域网格、边界、参数、源汇项等进行定义，然后调用解算器进行数值运算，由输出模块得到三维可视化结果。Visual Modflow软件求解数学模型的过程为依次调用三个独立的模块：输入模块、运行模块和输出模块，实现模型输入、求解运算和结果输出。

（1）模型输入

1）模型定义

永兴岛淡水透镜体水流模拟采用饱和变密度非稳定流，溶质运移模拟为非吸附模拟，不考虑化学反应。鉴于氯离子是水体中一种常见组分，也是海水的表征性元素，采样简单，分析方法成熟且精度高；氯离子又是地下水中必然出现的化学成分，在地下水环境中比较稳定，在研究区内有明显的变化趋势和规律；国内外许多较成功的水质模拟实例中，多数都选择氯离子作为模拟因子。因此，在永兴岛淡水透镜体水流模拟计算中，选择氯离子作为溶质运移模拟计算因子。

2）坐标系的选取

以概念模型底部平面作为基准面，取直角坐标系：x轴正向向东，坐标范围为0～1 928 m；y轴正向向北，坐标范围为0～1 284 m；z轴正向向上，坐标范围为0～30 m，原点位于基准面上。坐标选取如图6.11和图6.12所示。

图6.11　平面坐标与网格剖分

3）研究区域的空间离散

研究区域的空间离散就是根据含水层水文地质特征（岩性、孔隙率和给水度等），将所研究的区域剖分成相互联系的有限个小的水力水质均衡区（即计算单元），在满足一定精度条件下，确定计算单元的形状和大小。每个单元中含水层水力水质参数视为常数。

图6.12　垂直平面坐标与网格剖分

按照这一思路，将永兴岛的外形图导入软件，作为模型剖分的底图。软件将根据模拟区域的几何形状大小和用户要求自动完成三维剖分，再根据模拟范围定义有效单元网格（计算区域以内）和无效单元网格（计算区域以外），剖分范围如图6.6所示。剖分包括平面上的网格剖分和垂向上的分层剖分。平面上剖分为50行、50列，在垂向上分为51层，其中海平面以上的范围（5 m）为一层，海平面以下的范围（25 m）剖分为50层，剖分结果如图6.11（图上按每2行×2列给出）和图6.12所示。另外，地表高程以ASCII形式输入到模型中，地表测点以外的高程采用Kriging插值法求出。

4）时间离散

水流模拟周期36 500 d，初始时间步长Δt为0.001 d，增加因子为2，最大时间步长Δt为2 d。

5）边界条件

将6.2.2节提出的边界条件输入到软件中，如图6.13和图6.14所示。求解区域上边界为定流量边界和降雨补给浓度边界；下边界为隔水边界和定浓度边界；侧向边界海平面以上区域为隔水边界和定浓度边界，侧向边界海平面以下区域为定水头边界和定浓度边界。

6）初始条件

初始水头采用统一赋值的方式输入，即模型中各点的初始水头均赋值25 m。初始浓度进行了分区处理，如图6.15所示。

图6.13　水流边界

(www.xing528.com)

图6.14　浓度边界

图6.15　初始浓度分布

7）参数赋值

①渗透系数　结合地质勘测资料和实验室研究结果，不整合深度以上的渗透系数Kx＝Ky＝Kz＝110 m/d，不整合深度以下的渗透系数参考文献取Kx＝Ky＝Kz＝500 m/d。

②弥散度　采用一维连续注入示踪剂试验方法，测得珊瑚沙样本的纵向弥散度为0.29 m。如前所述，弥散具有尺度效应，实验室数据会小于野外实际的数据，因此，实验室求得的弥散度仅为一个参考值。结合野外渗流区域的实际情况，取纵向弥散度DI＝5 m，横向/纵向弥散度比率取0.5，垂向/纵向弥散度比率取0.05。

③其他参数　其他参数均由实验室试验获得，实验室无法完全模拟实际的水文地质条件，而且存在随机性的误差，所以，其数据只能作为参考，根据试验结果，并参阅其他文献，取弹性贮水率为1×10-5，重力给水度为0.31，孔隙率为0.40。

此外，海水密度ρs取1 025 kg/m3，盐浓度cs取35 g/L；淡水密度ρ0取1 000 kg/m3，盐浓度c0取0 g/L。模型选用氯离子作为模拟因子，SEAWAT程序可自动完成氯离子浓度与海水盐浓度间的转换。

（2）求解运算

模型的数值求解调用SEAWAT 2000程序包。程序包运行设置包括模拟进程、SEAWAT Flow和SEAWAT运移设置。运行设置是模型模拟的开始，为模拟地下水流动、溶质运移以及参数估算所选择的数值程序进行参数设置。根据不同的计算需要，设置内容包括水流和运移计算器种类的选择，设置内外最大迭代次数、运移计算的初始步长和最大步长等。

1）模拟进程设置

进程设置包括变密度流VDF（Variable-Density Flow Process变密度流计算）设置、水流设置和运移设置。对于既有水流又有溶质运移的变密度流，在VDF和溶质运移IMT（Integrated MT3DMSTransport Process溶质运移计算）模拟中，隐式耦合计算比显式耦合计算稳定且精确度更高，因而采用隐式耦合。在变密度流模拟中，节点间的密度计算采用Central-in-space法，这是一个根据相邻单元长度对密度进行加权计算的方法。模拟计算的起始时间步长取0.001 d，为软件系统默认设置。在水流设置中，水流计算采用预处理共轭梯度法（PCG）来求解模型产生的联立方程组，能模拟线性和非线性地下水流条件。该方法计算收敛速度快，结果可靠。收敛性取决于水头变化标准和残差标准（计算单元间的流量差），对于SEAWAT水流模拟，水头变化标准默认值为10-6，残差标准默认值为10。在运移设置中，运移计算采用隐式共轭梯度（GCG）解算器，选择GCG解法，弥散项和汇源项都默认采用中心有限差分隐式求解，从而没有任何稳定性约束。对流计算采用中心有限差分法，孔隙率采用有效孔隙率，Courant数取0.75。Courant数是表示一个示踪粒子在同一运移时间内从各个方向，允许进入网格单元的所有单元数目，通常为0.5～1。

2）SEAWAT Flow设置

在这一设置中，除设置时间步长、初始水头、各向异性、输出控制和列表文件选择外，还需进行补给、干湿设置和模拟层设置。补给设置是指地表补给赋予的模型分层，可以是顶层、指定层和最上有效层；干湿设置是指疏干单元重新湿润的选择。在非稳定模拟过程中，如果水头下降低于网格单元底部的水位时，潜水含水层中的单元就变为干枯单元，在余下的模拟中这些单元将不参与计算；若调用干湿交替选项，可以使这些干枯单元重新赋水，这样使计算更准确。在永兴岛淡水透镜体的模拟计算中，补给赋予最上有效层，模拟层设置为潜水含水层，疏干单元可以重新湿润，由计算程序自动判别。

3）SEAWAT运移设置

设定输出时间步长，给初始浓度赋值，设定浓度转化与查看功能，实现浓度单位的转换与预览。设定输出时间步长，就是将溶质运移的结果按时间段输出，控制不同时刻浓度输出数据的精细化程度，有利于查看浓度参数随时间的变化情况。

（3）结果输出

在Visual Modflow 4.1中，按输出控制进行设置，可将不同时刻的水头、降深、流量、浓度、流线、速度等信息以等值线（面）的形式显示出来，或直接导出数据。在3D模式下，还可输出这些参数的动态视频。

设置输出控制信息后，便可进行模型的编译和运行，并输出计算结果。永兴岛淡水透镜体的模拟计算结果的输出，包括以氯离子浓度为600 mg/L界定的淡水透镜体外形、水头和浓度分布、透镜体的最大厚度和贮水量等。其中，淡水透镜体贮水量按下式计算：

式中，V为贮水量；μ为给水度；v为淡水含水层总体积；s为x、y平面上差分网格的面积；hf，i，j为差分网格i，j对应的淡水水头；面积s的计算利用AutoCAD里的“创建面域”——“属性查询”功能得到。