基于有限差分模型模拟地下水-湖泊相互作用的研究

龚嘉临， 朱 松， 李非里， 金赞芳

(1．浙江工业大学环境学院， 浙江 杭州 310032；2.浙江省水利水电勘测设计院， 浙江 杭州 310002)

0 引言

目前，全球几乎所有湖泊面临着富营养化及其引起的蓝藻水华问题，严重影响着湖泊的水质和生态功能[1-2]。在国内，湖泊富营养化问题也早已成为一个较为突出的环境问题，如太湖、滇池、鄱阳湖等大型湖泊均已处于富营养或重富营养状态。而且，无锡饮用水供水危机等事件的发生也表明了我国湖泊水体富营养化问题的严重。然而，尽管发达国家对湖泊的富营养化研究和治理已经经历了半个世纪，但污染湖泊的水质和生态恢复至今仍未能完全实现。这是因为湖泊的富营养化过程是一个不可逆且长期产生影响的过程[3]。首先，造成湖泊富营养化的外源有很多种，包括农业污染、工业污染、生活污水，大气沉降、地下水入湖等。其次，湖泊底泥和腐烂的水生植物等内源物质也会将富集的营养物质释放，从而影响湖泊。而在地表水源受到普遍污染的现在，人们进一步开发并利用地下水，导致浅层地下水逐渐成为氮元素的贮存库。其中，硝酸盐污染更是已成为地下水污染的重要问题[4-5]。由于地下水与湖泊紧密联系，存在着水力交换，所以无论是哪种系统发生变化，都会对关联的另一系统造成影响。

随着对湖泊生态系统的不断研究，湖泊周边地下水对湖泊富营养化问题的影响受到了国内外研究人员的广泛关注。自上世纪70年代，国外的研究人员就开始了对地下水与湖泊之间相互作用关系的研究。由于赋存于地面下的岩石空隙中的地下水在一般情况下不能被人们直接观察，所以研究人员只有通过在现场的水文地质勘察与地下水动态的持续监测才能了解区域地下水系统的存在形式及运动规律。在研究初期，研究人员对地下水-湖泊系统的研究手段主要是利用渗流计和压力计等现场监测仪器来估算浅层地下水和湖泊之间的交换量。

随着计算机技术的发展和利用，越来越多的研究人员使用可视化模型来定量化模拟区域地下水和关联湖泊之间的相互作用[6-8]。其中，有限差分模型MODFLOW作为应用最广泛的模型软件，受到了众多研究人员的青睐。本文通过对地下水水流模型MODFLOW、溶质运移模型MT3DMS和湖泊模拟方法的介绍，总结其功能与应用特点。

1 有限差分模型MODFLOW简介及特点

MODFLOW模型（Modular Three-dimensional Finitedifference Ground-water flow model）最初是由美国地质调查局（US Geological Survey）的MCDONALD等在上世纪80年代所开发出来的一款数值模拟程序，用于模拟孔隙介质中地下水流动。

MODFLOW基于离散化方法，将连续的区域空间和时间分为一系列离散的格点。在一定精度下，将区域空间离散为有限个规则单元的集合体，使得每个网格单元上的各种参数值近似为常数。在离散化之后，整个区域的计算问题就被概化为了有限个单元网格的计算。结合时间离散之后，利用有限差分方程构成线性方程组。对线性方程组联立求解之后，获得各个离散点上的水头数值近似解。这种方法可以较好地反映区域地下水的状态，且具有较高仿真度。

MODFLOW的一个最显著特点是其模块化的结构[9]。它在主程序的基础上可以整合一系列相对独立的子程序包，而这些子程序包是由具有类似功能的子程序所组成的。同时，研究人员可以根据实际情况选用其中合适的子程序包对地下水系统进行数值模拟。这种模块化结构让整个程序更容易理解、修改，并可对其添加新的子程序包，如河流子程序、模拟水平流动障碍的子程序包等。模块化的结构使得MODFLOW的功能和应用范围大大拓展。模块化结构的MODFLOW不仅仅局限于最初的模拟地下水在孔隙介质中的流动，在合理使用的情况下可以解决裂隙水的流动等实际问题。

目前，使用有限差分模型MODFLOW的软件主要有GMS（杨百翰大学环境模型实验室和美国军工水道实验室合作开发），Visual Modflow（加拿大Waterloo水文地质公司开发）和Visual Modflow flex。随着地下水数值模拟软件的广泛应用，国内外研究人员主要通过MODFLOW模拟和预测区域地下水环境的动态变化，实现对实际问题的评价与动态预测，从而制定合理的措施解决问题[10]。SONG等[11]通过构建模型来确认开采计划是否合理，从而确保巴拉素水源地的安全。马驰等[12]根据模型结果预报西华水源地地下水位，获得合理的水源地开采量，为西北地区地下水资源的评价开创了先例。尉鹏翔[13]利用Visual Modflow软件对北京某污染场区进行建模（利用MODFLOW和MT3D），分析了污染物的迁移，为该污染场区地下水评价和污染治理提供了科学依据。ELKRAIL[14]结合GIS和MODFLOW对地下水的脆弱性进行评价，用于指导松花江所在区域地下水资源的开发利用和保护。

因为MODFLOW的最终目的是为了模拟实际的水文地质情况，故在使用该模型之前需要尽可能多地收集研究区的水文地质资料。在此基础上，整个模型的构建步骤见图1。

图1 MODFLOW模型构建一般步骤

2 溶质运移模型MT3DMS

为了研究地下水和湖泊富营养化的关系，在利用MODFLOW获得水流模型的基础上，研究人员还需要根据溶质运移模型获得相应的水质模型来进一步估算地下水对湖泊的污染负荷。

20世纪90年代，ZHENG开发了模块化的三维溶质运移模型MT3D。这款模型开发完成之后就广泛地应用于地下水中的溶质运移模拟和研究。随后，可插入多重污染物组分生化反应模块程序的MT3DMS由ZHENG等[15]于1998年开发。MT3DMS的功能在MT3D的基础上得到了极大的扩散，其独特之处在于囊括了有限差分法、基于粒子示踪的欧拉-拉格朗日法以及高阶有限差分TVD法3类主要的溶质运移解算法。之所以将3种解算法综合起来是因为没有一种单独的方法可以使用所有的运移条件[16]。

与MT3D一致，MT3DMS也是以块体中心法有限差分水流模型（如MODFLOW）为背景开发的。同时，MT3DMS保留了MT3D的模块化结构。而且和MODFLOW一样的是，MT3DMS同样是模块结构。其模块化的结构可以单独模拟运移的对流、弥散/扩散、源汇问题和化学反应问题，有利于合理分配内存；处理其它运移过程与化学反应的模块可以随时插入程序，而不需要改动现有代码。因此，MT3DMS可以模拟含水层中的对流、弥散-扩散过程，也可以模拟地下水污染物质在运移过程工程的生物反应和化学反应。赵贝等[17]通过MODFLOW和MT3DMS构建了河北省东南部地区的一个垃圾填埋场的三维数学模型，并模拟预测了该区未来30 a内的特征污染物Cl-的迁移情况。模拟结果表明，污染物浓度随着离污染源的距离增大而减少，主要危害地下水流方向上的地表水体和浅层地下水；湖泊对污染物具有持续捕获作用；开采地下水会影响污染物的迁移速率和方向。GUSYEV等[18]利用MODFLOW，MT3DMS和同位素氚对新西兰的西陶波湖捕集区进行地下水年龄模拟，用于构建流域地下水年龄分布。

3 湖泊模拟方法

自然环境中的地下水和湖泊存在着水力交换，从而导致地下水影响着湖泊的水平衡、营养化水平和酸缓冲能力。而湖泊周边的地下水是作为补给源还是排泄汇，这取决于湖泊与周边地下水的水力梯度。1979年，BORN等[19]依据湖泊流域的水文循环特征，对地下水-湖泊交换定义了3种不同的类型，见图2。第1种类型是地下水穿流，湖泊的部分湖底接受地下水，而部分湖底发生渗漏；第2种类型是地下水进入湖泊，对湖泊进行补给；第3种是湖泊水对地下水进行补给。一般情况下，地下水-湖泊交换类型随季节变化而发生变化。其中，第1种类型在山区流域中出现的频率较多，且在具有恒定流向地表水体的平原地区也常出现，具体表现为湖泊在上游的入湖区接受地下水和下游的出湖区渗漏进入地下水；第2种类型常出现在沙漠地区或蒸发量远远大于降雨量的湖泊流域，当湖泊的周边地势较高时，湖泊位于流域的地势较低处，湖泊水体的出流主要依靠蒸发的作用，入流则主要来自部分降水和地下水流；第3种类型通常在地表水位主要补给的湖泊流域出现，湖泊水体的出流受地势影响，或湖底的渗透性能较好，从湖底补给的地下水能够较快地排走，从而使含水层中的地下水位维持在比湖泊水位更低的水平。

图2 地下水-湖泊的交换形式

研究地下水和湖泊的相互作用有利于了解整个系统，从而为解决实际问题提供依据。而在数值模拟中，研究人员也可以通过以下几种方式模拟地下水和湖泊的相互作用[20-22]。

3.1 湖泊区域作为高渗透系数的含水层

当进行MODFLOW模拟时，将湖泊作为构建的网格模型含水层的一部分，同时设置其所在区域的渗透系数为高渗透系数。李勇[23]和莫美仙[24]结合该方法和第2种方法，把滇池作为模型含水层中的一部分，设置湖泊所在区域单元具有高渗透系数，且将湖泊区域设置为边界，从而来模拟滇池。前者利用ZONGBUDGET模块计算出滇池东侧地下水入湖水量为8 073.8 m3/d，硝态氮和总磷入湖量分别为32.4和0.91 t/a，北侧地下水入湖水量为8 463.6 m3/d，硝态氮和总磷入湖量分别为54.72和0.832 t/a；后者计算出昆明盆地北部地下水入湖量在6.5×106～6.9×106m3/a之间，并根据地下水携带的氮、磷浓度平均值乘以入湖水量估算氮和磷的入湖量范围分别为27.13～ 406.26，1.84～ 10.04 t/a。

这种方法模拟较为便利，可以计算湖泊水位，并模拟稳态和瞬态问题。但其缺陷在于只能模拟渗漏湖泊，无法包含其它的地表水流动，且容易出现质量不守恒问题。

3.2 湖泊作为边界条件输入模型

湖泊作为边界条件输入模型即把湖泊这一系统看成是MODFLOW模拟区域的其中一个已知边界。在这种情况下，研究人员可以利用MODFLOW中的定水头边界或通用水头边界来表示湖泊。潘田[25]将研究区东北部连接的太湖定义为定水头边界，以太湖水位作为水头条件，从而实现研究区地下水与太湖的水力联系。同时，结合MT3DMS进行溶质运移模拟获得太湖流域地下水硝态氮补给太湖。结果表明，太湖周边浅层地下水补给太湖水量为0.2亿m3/a，占总排泄量的6%；太湖周边浅层地下水中的硝酸盐氮入湖负荷约为14 245 t/a，占补给太湖氮污染负荷的27.8%。郭玉川[26]将博斯腾湖作为变水头边界，湖泊水头值根据实际水位分时段赋值，通过构建模型，从而分析模拟了基于生态安全水位的不同调控方案的水均衡情况，对当地地下水位的调控方案进行优选。

该方法模拟的湖泊水位稳定且利于模拟渗漏和排水湖泊。但是，由于湖泊水位已经作为固定边界输入，在整个模型模拟过程中无法体现湖泊的水位变化且一般不包括其它地表水流。同时，当其作为定水头边界输入模型的时候，湖泊边界其实相当于无限的源汇项。

3.3 河流程序包代替湖泊模块

河流模块作为MODFLOW边界条件中的常用模块，也可以用于模拟其它地表水体。这种方法与通用水头边界条件一样都是以外部应力形式导入模型。但与通用水头边界导入不同的是，使用河流模块模拟湖泊等其它地表水体的时候可以定义水体范围，同时在结果中可以显示水域的变化情况及流场。而缺陷在于可能出现质量不守恒情况。

3.4 LAK3程序包

上文已说明MODFLOW是一种模块化的结构，这为湖泊程序包的开发和使用提供了便利。LAK湖泊程序包（LAK1,LAK2,LAK3）从上世纪90年代开发至今，已被国内外研究人员广泛应用。KROHELSKI等[27]利用 MODFLOW 和LAK3程序包来模拟威斯康星州西北部Dane镇的3个浅水渗透湖泊（Fish湖，Mud湖和Crystal湖）及其周围的浅层地下水。模拟结果表明，Fish湖水位在抽水条件下会降低；在没有持续抽水的情况下，湖泊水位又会在几年内恢复。FEINSTEIN等[28]模拟了威斯康星州西北部Waukesha镇的湖泊、井和浅层地下水的相互作用。模拟结果表明，地下水是Beaver湖最重要的入流补给部分（等于总流入的59%）；对于Pine湖和North湖来说，地下水仍然是重要的来源（分别占总流入的16%和5%），但小于降水和地表水的贡献。当在地表61 m以下以213.85 L/min速率抽取地下水5 a，对湖泊水位的影响也并不大，模拟的Pine湖的水位和出流量仅分别减少了0.91 cm和3%。但是，加入抽水速率增加到910 L/min，模拟的Pine湖的出流量明显减少，相对基本流量条件下降了14%。在我国，应用LAK3模块的情况较少。程春龙等[29]利用MODFLOW和LAK3模块建立了天津拟建宁车沾水库的耦合模型。模拟结果表明，水库水位下降速率受到当地的蒸发强度影响较大，受渗漏影响极小，年水位降幅1.16 m。

LAK3程序包作为最新的模拟湖泊程序包，具有进行稳态和瞬态模拟的功能，且可以用于模拟地表水流和溶质运移。同时，和河流程序包一样，LAK3程序包的湖泊模拟可以在结果中显示。与前面3种方法相比，LAK3输入较为复杂，包括湖泊水位、底泥层厚度、渗透系数、降雨、蒸散、径流等参数。其缺陷主要有2点：①当渗透系数较高时，模拟具有潜在的不稳定性；②在没有用户干涉的情况下，可能出现错误结果。

4 存在的问题及发展趋势

根据文献检索发现，近年来有关区域地下水-湖泊系统的三维可视化模型研究数并不少。这些研究涉及的湖泊及周边地下水的水文地质条件复杂多样。利用MODFLOW等模型的构建为区域地下水和湖泊的水资源预测规划和管理提供了重要的依据。但是，综合分析可以发现以下问题的存在。

（1）研究范围常以人为边界划定。由于实际地质条件、资料收集、工作量和计算机容量所限，研究人员在进行模拟研究的过程中往往很难找到自然边界作为边界条件，这就导致模拟区域往往比实际的地下水含水层系统要小。这种情况下，MODFLOW模型常采用人为边界。因为人为边界的引入，也导致了模型预测的困难。

（2）水文数据资料收集和验证困难。和国外相比，国内的野外实验基地没有形成规模，这也导致水文地质资料往往没有长期的监测、考证资料，无法形成数据积累。这也导致了国内构建的模型难以进一步修正优化，达不到长期预测的效果。

（3）数据处理往往以多年平均、短期观测资料为主。同样是因为实际地质条件、资料收集和工作量所限，研究人员往往只能尽可能将收集到的水文地质信息进行合理概化，以求在模拟的时候接近实际情况。因此，这也造成了根据数据资料的差异，不同研究人员的研究成果往往存在一定差异。同时，这也必然导致研究成果与实际的偏差，以至于很多情况仅仅限于科学研究。

5 结论

MODFLOW这类用于模拟实际环境状态的模型离不开各种大数据的支持。随着我国社会经济的不断发展，对于城市及周边地区地质数据的调查和完善提出了新的要求。进入21世纪以来，中国地质调查局和各地方政府合作开展的“城市地质调查项目”就是为了解决城市发展过程中的地质问题，并已在多个重要城市展开试点调查工作。这为将来的进一步地质研究提供了便利。而随着互联网的发展，大数据资料库的建立使区域水文地质资料得到了不断完善和扩充，加快了信息交流。同时，计算机技术的发展也使计算机容量不断增大，为更大范围的区域地下水模拟提供了可能。长远来看，数据资料的扩充和共享交流对模型研究具有重要意义。