地下水模拟方法与应用软件研究进展

卢洪健，王卓然

(水利部信息中心，北京 100053)

0 引言

地下水约占地球液态淡水总量的99%，其踪迹几乎遍及全球，具有巨大的社会、经济和环境效益。目前世界上一半的居民生活用水来源于地下水，约有25%的农业灌溉来源于地下水开采[1]。尽管地下水如此重要，人们却对它普遍缺乏认识，导致地下水的价值被严重低估。近半个世纪以来，人类长期过度开采，加之气候变化的影响，导致当前世界有超过一半的地下含水系统处于超采状态，有些国家的地下水水位正在以十分危险的速度下降。地下水超采不仅严重威胁饮用水安全和依靠地下水灌溉农产区的粮食产量，而且会导致污染加剧、生态支撑能力下降及地面沉降等生态环境问题[2]。因此，准确的地下水模拟是地下水资源精细化管理与决策的重要支撑。

地下水系统由地下水含水系统和地下水流动系统组成，由于其一般不可见、受多种补排条件综合作用、受地质条件影响等特殊性，对地下水系统的研究远比地表水系统困难。自20世纪60年代以来，数值模拟开始应用于地下水计算中，地下水数值模拟的理论与方法得到了长足的发展[3-4]。准确的地下水过程模拟结果不仅可以用于地下水污染评估、地下水最优管理，还能为解决地下水资源开采安全性等问题提出决策支持。本文阐述了地下水模拟过程的理论基础，介绍了地下水数值模拟主要求解方法，比较了当前主流地下水模拟软件及其在我国的应用情况。

1 地下水模拟

地下水模拟是采用地下水数学模型，通过将概念模型转化为控制方程的形式，在特定的边界和初始条件下，用来模拟和描述真实世界的地下水流运动，并可以通过执行计算机程序来求解[5]。地下水模拟模型是一种非均匀模型，通过近似调查的地下水系统，采取简化假设，例如均匀性、各向同性、流动方向、含水层几何结构、污染物运移机制等[6]。模型的性能和效率取决于数学方程对被建模物理系统的近似程度，模型精度则取决于对地下水系统的概化和理解程度，以及数学方程推导过程中的假设[7]。

地下含水层建模过程一般包括以下步骤[8]：(1)确定表征含水层物理结构和系统条件的参数；(2)利用特定点的现场资料估算水文地质参数；(3)利用插值/外推方法估计参数的空间分布；(4)利用所有估计参数和现场数据构建概念模型；(5)利用地下水流动方程表达系统条件，建立描述概念模型的数学模型；(6)将数学模型转化为数值模型，以确定含水层响应，包括水头或污染物浓度；(7)利用数值法对生成的模型进行求解；(8)通过敏感性分析，确定校正模型时需要调整的关键模型系数；(9)利用可获取的现场数据，对模型进行校准，以更准确地预测地下水系统的变化；(10)通过验证模型，消除数值近似带来的误差；(11)采用模型评估既定管理策略对含水层恢复和地下水系统优化利用的影响。

地下水管理通常需要预测地下水流、地下水位、溶质运移和模拟自然或人为应力，地下水模拟模型被广泛用于此类预测和模拟中，这些模型通常需要解偏微分方程。数学模型可以是确定性的、随机的(统计的)或两者的结合。在随机模型中，根据发生的概率向模型提供一系列的预测，可以帮助评估系统的不确定性；确定性模型则是基于已知系统和过程的因果关系而建立，广泛用于解决区域地下水问题，可以分为解析型和数值型[9]。解析模型利用若干简化假设对地下水系统进行快速的初步分析，不能用于解决区域形状不规则、区域异质性和复杂边界条件的问题。对于这种情况，数值模型使用计算机程序来解决更复杂的问题[10]。

2 地下水数值模型求解方法

数值模型用于求解代表地下水流运动的偏微分方程，并给出近似解，模型的主要特点包括[11]：(1)模型仅在为问题定义的空间和时间域(离散值)中的指定点求解，这些点被视为整个区域的不连续状态变量；(2)描述地下水流运动的偏微分方程在某些点通过一组数学方程转化为状态变量的离散值；(3)其解是针对各种模型系数的一组指定的数值，而不是这些系数的一般关系；(4)使用计算机程序来求解大量必须同时求解的方程。目前，已发展出有限差分法(finite-difference method，FDM)、有限元法(finite element method，FEM)、有限体积法(finite volume method，FVM)、边界元和粒子追踪等数值求解方法[12]。

2.1 有限差分法

有限差分法(FDM)是根据含水层的特征和条件，将差分方程中的偏导数在小范围内用代数表达式进行变换，问题域被分割成一系列被称为节点的离散点，用一组离散的点替换连续介质，并为每个节点分配各种水文地质参数。FDM可用于时间和空间离散化，利用定义参数间时空关系的差分算子来替换偏导数，所建立的模型在每个节点上通过获取该节点上一组代数方程的解来求解，通常会采用一些迭代方法来求解简化方程[13]。

该方法利用时间步长开始时的初始条件以及时间步长期间发生的含水层抽水或回灌率来计算时间步长结束时的未知水头。因此，在每一个时间步长中需要同时求解大量的方程，这使得该模拟技术对于需要长期模拟大型含水层的地下水规划和管理目标来说非常耗时。

FDM的优点是易于跟踪包括复杂加载路径和高度非线性行为的复杂系统，易于理解和编写程序，因此是求解大型非线性地下水流运动问题的一种经济方法。然而，对于不规则几何域，FDM的使用是困难的。通常来讲，具有规则网格系统的传统FDM方法存在形状域不规则、边界条件复杂、材料非均质性等缺点[14]。

2.2 有限元法

在有限元法(FEM)中，不规则形状区域可以划分为一组具有不同尺寸或形状的单元。为了反映状态变量或参数值的变化，可以更改元素大小。采用直接法、加权残值法和变分法等方法对单元的偏微分方程进行近似，以获得一组代数方程。因变量的分段连续表示以及地下水系统的参数(可能)可以提高数值近似的精度。对于许多地下水问题，有限元法优于经典的有限差分模型。具有不规则形状区域、复杂边界条件和材料非均质性的地下水问题可以使用有限元建模，而FDM意味着复杂的插值格式来逼近复杂的边界条件[15]。

将问题域划分为若干非重叠单元，作为有限元法求解地下水流动或溶质运移问题的第一步。将问题域替换为一系列节点和离散单元或有限元网格，这些元素通过将两个或多个节点连接在一起。然后，为每个节点分配一个节点号，为每个元素分配一个元素号，元素可以具有任意大小的一维、二维或三维，在每个要素中，应规定地下水流动特征。利用地下水流动和溶质运移过程的知识绘制网格是一种有助于以合理的计算负担和可接受的精度获得地下水系统解的方法。这可以通过流动或运输过程可视化和流网来实现。可以使用不同的有限元网格类型进行建模，结果可能相似。因此，建模时没有唯一的网格类型和大小选择。

与粗网格相比，使用细网格进行有限元建模可以得到更精确的解，从而导致精度较低。精细网格有更多的节点，需要更多的计算工作来获得解决方案，因此，它是计算负担和建模精度之间的权衡。网格大小可以通过使用更细的网格重复计算来确定，并查看结果的变化有多大。在建模的第一次重复中，可以使用几个节点生成粗略的有限元网格。因此，只需很少的计算工作就可以导出一个解。在重复建模过程中，可以准备更精细的有限元网格，这需要更多的计算工作，并导致更精确的解决方案。

2.3 有限体积法

在有限体积法(FVM)中，通过将计算函数划分为控制体积，并将加权函数设置为统一于控制体积，生成了许多加权残差方程。在控制体积中，残差的积分必须等于零。问题域被划分为多个控制卷，没有重叠。将微分方程积分到围绕每个网格点的一个控制体积上。分段连续性表示网格点之间的变化，用于计算积分。结果是包含一组网格点的离散化方程。在离散化方程中得到了有限控制体积的守恒原理。质量、动量和能量等量的积分守恒在任何控制体积组和整个问题域上都是完全满足的。对于任意数量的网格点，即使是粗略的网格解也能显示出精确的积分平衡。在FVM中，不需要结构化网格，变量位于体积内，因此可以无创地应用边界条件[16]。

在FDM中，偏微分方程中的一阶导数由相邻节点的自变量值之间的差来近似，考虑节点之间的距离，并考虑两个连续时间的时间步长增量的持续时间。在FEM中，因变量和参数的函数用于评估偏微分方程(PDE)的等效积分公式。虽然每种方法都有一些优点和缺点，但由于概念和数学上的简单性，FDM通常更容易编程。在上述方法中，有两种求解偏微分方程的方法来获得因变量的网格点值。FEM中使用的一种方法是，头部变量由网格点值组成，形状函数用于网格点之间的插值。另一方面，在FDM中使用，忽略网格之间的水头变化，方程式包括水头的网格点值。

FDM在不规则含水层边界和含水层内区域参数的紧密空间近似方面具有灵活性。然而，与常规矩形有限差分网格相比，不规则有限元网格的网格生成、输入数据集的规格说明和构造要困难得多。FVM是一种将偏微分方程表示为代数方程并进行计算的方法。与FDM类似，在网格几何体上的离散位置计算值。与FDM相比，FVM的一个优点是它不需要结构化网格，在粗非均匀网格和网格移动以跟踪界面或冲击的计算中尤其强大。

3 地下水模拟软件及应用

3.1 常用地下水模拟软件

MODFLOW(Modular Three-dimensional Finite Difference Groundwater Flow Model)是是由美国地质调查局于20世纪80年代开发出来的三维地下水流数值模拟模型[17]。该软件以有限差分法为基本原理，即在不考虑水的密度变化条件下，孔隙介质中地下水在三维空间的流 动可以偏微分方程来表示。通过模拟不规则形状水流系统中的定常和非定常流，其中含水层可以是承压、非承压或承压和非承压的组合，可以模拟来自外部应力的流量，例如流入井、地表补给、蒸散、流入排水沟和流经河床的流量。该软件可以模拟特定的水头和通量边界，还能够模拟穿过模型外边界的水头相关流量，从而允许以与模型区域外水源和边界块之间的当前水头差成比例的速率向模型区域内的边界块供水。除了模拟地下水流动外，MODFLOW的应用范围已经扩展到溶质运移和参数估计等功能。

在MODFLOW 软件基础上，加拿大Waterloo Hydrogeologic Inc.应用现代可视化技术开发研制出Visual MODFLOW，于1994年8月首次在国际上公开发行。Visual MODFLOW 以其求解方法的简单实用、适应范围的广泛及可视化功能的强大正成为最有影响的地下水模拟平台环境。然而实践也证明，对于复杂的地质条件、不饱和流动、密度变化的流动( 海水入侵)、热对流等棘手的问题，Visual MODFLOW 往往并不适合。

基于有限元法的FEFLOW (Finite element subsurface FLOW system ) 软件由德国WASY 水资源规划和系统研究所于1979年开发出来[18]。FEFLOW是现有的功能最齐全最复杂的地下水模拟软件包之一，用于模拟多孔介质中饱和及非饱和地下水流与污染物的运移。FEFLOW软件具有图形人机对话、地理信息系统数据接口、自动产生空间各种有限元网格、空间参数区域化及快速精确的数值算法和先进的图形视觉化技术等特点。由于它是为满足专 门从事复杂地下水模拟工程的专家对技术的要求而设计的，对含水层分层、单元剖分、离散 点插值、数据转换、边界条件赋值、河流边界、含水层均衡项等高效处理的特点，使其适宜 于大区域地下水流模拟。

表1 国际主流地下水模拟软件

地下水模拟系统( Groundwater Modeling System) , 简称GMS, 是美国Brigham Young University的环境模型研究实验室和美国军队排水工程试验工作站在综合MODFLOW 、 FEMWATER、MT3DMS、RT3D 、SEAM3D、MODPATH、SEEP2D、NUFT、UTCHEM 等已有地下水模型的基础上，开发的一个综合性、用于地下水模拟的图形界面软件。其图形界面由下拉菜单、编辑条、常用模块、工具栏、快捷键和帮助条6 部分组成，使用起来非常便捷。由于GMS 软件具有良好的使用界面，强大的前处理、后处理功能及优良的三维可视效果, 目前已成为国际上最受欢迎的地下水模拟软件。此外，加拿大Waterloo Hydrogeologic Inc.还开发了Visual Groundwater、HydroGeo- Analyst、WHIUnSat Suite，丹麦水利研究所开发了MIKE系列模型[19]，当前国际上主流地下水模拟软件功能及特点详见表1。

3.2 我国地下水模拟软件开发与应用

我国地下水模拟模型研制和软件开发起步相对较晚，20世纪70年代初，国内高校和科研院所开展地下水数值模拟研究，应用计算机语言，如Fortran、Algol等编写一些计算程序，基本为内部使用，未公开推广应用。林学钰等[20]开展了承压含水层中二维溶质运移和弥散的微机模拟(有限单元法)研究；薛禹群等[21]研制了地下水污染模型，并耦合海水入侵、大区域地面沉降模块；武强等[22]开展了河水与地下水耦合模型，地表水、地下水和土壤水之间的水力耦合模型研究；陈崇希等[23]开发了基于多边形网格的三维地下水流有限差分模拟系统PGMS；陆垂裕等[24]开发了基于流域/区域长时间尺度水量平衡分析的分布式水文模型，具有较强的物理基础，可详细模拟大气水、土壤水、地表水和地下水之间复杂的“四水转化”过程，上述研究为我国地下水模型研究奠定了坚实基础。此外，国内很多学者采用国外主流地下水模型，或进行二次开发，针对我国地下水超采最为严重的华北平原和水资源极度短缺的西北内陆干旱区，开展了大量的地下水模拟模型应用研究，为地下水决策管理发挥了重要的基础支撑。

4 问题与展望

近几十年来，我国在地下水模拟及应用上取得很多成果，解决了很多国民经济建设中急需解决的各类地下水问题，但理论成果和原始创新相对较少；目前地下水模拟新技术、新方法的开发主要由欧美发达国家主导，我国的地下水通用模型研制及软件开发还相对落后。2021年10月21日第748号国务院令公布《地下水管理条例》，自2021年12月1日起施行，标志着我国地下水管理正式步入法治时代，这对我国地下水模拟模型和应用软件的支撑能力要求也越来越高。因此，在水利高质量发展的新阶段，为进一步支撑地下水“评价-规划-管理-保护”各环节，支撑地下水“预报-预警-预演-预案”功能建设，支撑地下水超采治理水量-水位效果联合评价等，十分必要建立全国通用地下水模型。

全国通用地下水模型建设必须有自主研发的核心计算模型，模型需适应我国大多数地区地下水模拟，可以地下水量、地下水位模拟为主，兼顾溶质运移、海水入侵、地面沉降、多相介质流动、水-热耦合模拟、水生态模拟等管理应用需求，模型应经受实际工程检验，保证计算结果稳定可靠，同时应公开源代码和设计接口，以便后续改进完善。