地下水数值模拟模型

20世纪60年代初，随着电子计算机运算速度和容量的提高，数值模拟开始广泛地应用于大规模实际地下水流的计算。Pinder和Bredehoeft将Peaceman和Rachford提出的交替方向隐式方法用于地下水的计算，稍后又引入强隐式［13-15］。有限元法为我国数学家冯康创建。该方法于60年代后期，开始应用于地下水流计算中，如1966年Zinekiewiez把有限元用于二维稳定流计算，1968年Jevendel等进一步用有限元方法解非稳定流问题［16］；1972年引入等参有限元，1976年Capta等人用三维等参数有限元对多层地下水盆地进行了数值模拟；张宏仁和李俊亭（1979）发现了有限元求解地下水流问题时不满足局部质量守恒［17］，Neuman和Narasihuhan亦发现由于释水矩阵的非对角性所构成的缺点，因此提出了相应的改进措施［18］。
我国自1973年以来在地下水的数值模拟方面发展很快，1980年薛禹群和谢春红教授的《水文地质数值法》一书问世以来，孙讷正、陈崇希、张蔚榛、罗焕炎、陈雨孙、李俊亭、杨天行、林学钰等相继出版了多部有关地下水数值模拟的专著、教材和学术论文［19-33］。地下水数值模拟的应用已遍及与地下水有关的各个领域和各个产业部门。高校、科研院所与生产部门相结合，已运用数值模拟解决了很多国民经济建设中急需解决的各类问题，其中包括：水资源评价问题；地下水污染问题，水-岩作用和生物降解作用的模拟；非饱和带水分和盐分运移问题；海水入侵、高浓度咸水或卤水入侵问题；热量运移和含水层贮能问题；地下水管理与合理开发、井渠合理布局和渠道渗漏问题；地下水-地表水联合评价调度问题；地面沉降问题；参数的确定问题。它所涉及的地质条件多种多样，有潜水，也有承压水；有单个含水层，也有多个含水层存在越流的情况，以及种种复杂的地质构造和岩相变化的情况，由此，探讨了相应的模型概化与边界条件的处理。模型有二维的（平面的、剖面的），也有三维的，但以二维为主。水流模型有饱和的、非饱和的、饱和-非饱和的。
溶质运移模型在我国多数是处理低浓度的水质（地下水污染）问题。因此，由水流方程和对流-弥散方程分别组成的两个子模型可以独立求解，运动方程也以传统的达西定律为基础。只有少数研究海水入侵、卤水或咸水入侵和污水中高浓度污染物运移问题中，密度、粘度要由状态方程决定。此时，上述两个子模型要耦合起来求解。迭代法是解这类问题常用的解法。我国最早的三维可混溶海水入侵模型，是在20世纪80年代末期建立的。在我国这些海水入侵、卤水或咸水入侵模型以及热量运移模型、运动方程中，除了根据传统的达西定律考虑以水头梯度为基础的强迫对流外，还考虑了自然对流。卤水或咸水入侵由于浓度高，还考虑了由于粘滞性产生的切应力对水流运动的阻滞。
热量运移模型，已考虑了与热量运移有关的各种主要因素（对流、传导、热机械弥散、自然对流、水-岩间的热交换）。形变模型一般只考虑垂向变形（地面沉降或抬升），沉降计算是建立在由水流子模型求得的水头或水压强变化与水压强和含水层形变间有关关系的基础上的。
总之，我国地下水数值模拟的研究虽然起步较晚，但由于一开始在人员上就实行水文地质界与数学界的结合，在选题上一开始就结合生产实际中面临的问题开展研究，所以进展很快。在推动水文地质数值模拟工作中作出重要贡献的数学教授有：清华大学数学系肖树铁教授、南京大学数学系谢春红教授、山东大学数学系孙讷正教授、中国地质大学陈明佑教授和杨天行教授等。1978年12月首届地下水资源评价学术会议上就出现了一批理论和应用成果、以后继续沿着这条道路走下去，成绩越来越喜人，基本满足了国民经济发展建设的需要。一些成果接近或达到国际先进水平。本学科国际权威刊物《Water Resources Research》、《Journal of Hydrology》、《Advances in Water Resources》、《Ground Water》上均出现了中国学者的文章，两年一次的国际水资源计算方法会议（ICCMWR）等重要国际会议上也有了中国人的声音。在国内还成功地举办了地下水水流和污染模拟国际会议（1991），地下水和环境国际讨论会（1992）等学术会议。我国的地下水数值模拟开始走向世界。

1.2.1.1 解地下水流问题的数值方法
随着现代应用数学以及计算机技术的发展和广泛应用，数值模拟逐渐成为研究地下水运动规律、定量评价地下水资源以及模拟一些水文地质过程发生发展的主要手段。其研究范畴，由单纯研究地下水系统与自然环境系统之间的相互关系，扩大到研究与社会经济系统的相互关系；其研究内容，涉及饱和带、非饱和带和饱和-非饱和带［1］。
在计算机上利用数值法可以模拟各种复杂的水文地质条件，虽然用这种方式求出的数值解仍然是近似值，但仍能满足人们生产研究的要求［2］。和其他方法比较，数值模拟有很多优点，主要有［3］：①模拟在计算机上进行，不需要像物理模拟那样建立专门的一套设备；②有广泛的适用性，可以用于水量计算、水位预报以及地面沉降等的计算，各种复杂的含水层、边界条件、水流情况都能模拟出来；③修改算法、修改模型比较方便；④可以程序化，只要编好软件，对不同的具体问题只要按要求整理数据就能上机计算，并立即得到相应的结果。它的不足之处是不如物理模拟来得逼真、直观，且计算工作量大。
目前，解地下水流问题的数值方法有很多，但最通用的还是有限差分法（Finite Difference Method）和有限元法（Finite Element Method）。这两种方法的根本区别在于有限差分法是建立在用差商近似表示导数的基础上的，而有限元法是建立在直接求函数的近似解基础上的。除了这两种方法以外，还有特征线法（Method of Characteristics）、边界元法（Boundary Element Method）等，在此不再详述。
20世纪50年代有限差分法主要用于石油流动领域的计算，60年代中期拓宽了应用领域，用于解地下水流问题。有限差分法有许多优点：①对于简单问题（如均质各向同性含水层中的一维、二维稳定流问题）的数学表达式和计算过程比较直观、易懂；②有相应高效的算法；③对一般的地下水流问题来说解的精度比较高；④有广泛使用的商用软件，如MODFLOW、PLASM等。需要注意的是，对某些自然边界条件，有限差分法必须进行特殊处理，灵活性一般说来相对要差一些。因此，标准的有限差分法在近似不规则边界上不如有限元法方便（积分有限差分法能和有限元法一样处理不规则边界），对内部边界如断层带的处理以及模拟点源（汇）、渗出面和移动着的地下水面等，有限差分法也不如有限元法好。
有限元法于20世纪60年代后期引入地下水计算中。这种方法的优点有：①程序的统一性。有限元法对各种地下水流和溶质、热量运移问题，计算过程基本相同，程序编写比较方便，很多例子表明从解某一类问题的程序转换为解另一类问题的程序比较简单；②对不规则边界或曲线边界、各向异性和非均质含水层的处理比较方便；③单元大小比较随意，同一计算区内可以视需要采用多种单元形状和多种插值函数以适应水头、浓度等变量的变化或精度要求；④水流问题、溶质运移问题的解精度一般比较高。有限元法虽然有上述这些优点，也有一些缺陷，主要是局部区域质量不守恒，有时会影响计算结果。另一个是和有限差分法等共有的缺陷，即渗流速度、流量只能在先求出水头后，再由Darcy定律算出渗流速度，渗流速度乘以过水断面面积再得到流量［3］。这样做误差大，算不准，至今尚未彻底解决。2004年薛禹群等［4］、YE 等［5］又把数学上新出现的多尺度有限元法引入地下水领域，并得到初步应用，它不仅可以大量减少单元数，还能提高计算精度。因此，多尺度有限元是一种很有应用前景的方法。
1.2.1.2 地下水数值模拟软件
随着计算机技术的快速发展，使得复杂含水层系统中的地下水流运动及溶质运移的数值模拟变为可能。近年来，在人机交互、计算机图形学和科学可视化等技术的推动下，地下水数值模拟软件在质量上有了较大的发展和提高［6］。其中较有影响地位的有Visual MODFLOW、FEFLOW、GMS。
（1）Visual MODFLOW：由加拿大Waterloo水文地质公司在MODFLOW的基础上开发研制的Visual MODFLOW软件，是目前较流行且被各国同行一致认可的三维地下水流和溶质运移模拟评价的标准可视化专业软件系统［7］。该软件主要包括MODFLOW（水流模拟）、Modpath（平面和剖面流线示踪分析）、MT3D（溶质运移模拟）和Zone Budget（水量均衡计算）四大模块。界面设计包括三大彼此联系但又相对独立的模块，即前处理模块、计算模块和后处理模块。Visual MODFLOW以其求解方法的简单适用、适应范围的广泛及可视化功能的强大成为较有影响的地下水数值模拟软件，其使用范围越来越大。然而实践证明，它往往并不适合某些复杂的地质条件，如不饱和流、密度变化的水流（海水入侵）、热对流等棘手的问题。
（2）FEFLOW：FEFLOW是由德国水资源规划与系统研究所（WASY）历时20多年的研究，开发出来的地下水流动及物质迁移模拟软件系统［8］。软件问世以来，在理论研究和实际问题的处理上，经过了不断的发展、修改、扩充、提高，日趋完善。从20世纪70年代末至今，FEFLOW 经过了大量的测试和检验，成功地解决了一系列与地下水有关的实质性问题，如判断污染物迁移途径、追溯污染物的来源、海水入侵等，是功能较齐全的三维地下水模拟分析软件。
（3）GMS：GMS是由Brigham Young大学环境模拟研究实验室开发的较先进的、基于概念模型的地下水系统模拟软件［9］。GMS是唯一支持Tins、solids、钻孔数据、2D 或者3D地质统计学的系统，它也包括2D和3D的有限单元和有限差分模型。此外，它还封装了基于MODFLOW的水流模型、溶质运移模型MODPATH、MT3D和RT3D，以及基于有限单元法的FEMWATER模型等。其主要优点体现于：在前处理过程中，GMS软件可以采用MODFLOW 等模块的输入数据，同时MODFLOW 等模块的计算结果又可以直接导入GMS中进行后处理，实现计算结果的可视化。

一、模型计算范围及剖分

选用有限差分法建立地下水数值模型。模拟分析软件选用PMWIN（processing Modflow）和GMS（Groudwater Model System）求解地下水运动的定解问题，PMWIN是美国地质调查局开发的用于模拟和预报地下水系统的应用软件，它是一个以Modflow为核心的可以用来处理三维模型的软件（Wen-Hsing Chiang，2005）。PMWIN具有较好的数据导入界面，GMS具有较好的数据后处理可视化显示，结合两者优点进行模拟。模型计算范围，北起黄士台源，南至泾、渭河，西界起自19276km线，东界至19351km线，扣除其内不建模的部分，模型总有效面积为1513km²。以1km的均匀步长对模型进行剖分，其剖分网格实际上就是高斯-克吕格投影地图中的“公里网”。泾惠渠灌区地下水模型剖分图见图8-2。

时间剖分以自然月为时间步长。

图8-2 泾惠渠灌区地下水数值模型剖分及资源分布图 Fig.8-2 Groundwater numerical model split and resource distribution in Jinghui Canal Irrigation District

二、模型边界条件与地下水补、排要素的处理

1.侧向补给处理

模型的计算区为第四系松散沉积物潜水含水层。为简化模型，北部黄士台源洪流入渗放在模型北部边界上，其数量取多年平均值，忽略其随时间的变化。

2.降水入渗补给

根据灌区水文地质图，结合不同地形地貌单元降水入渗补给系数的取值。综合考虑包气带的降水量、士质类型、下垫面条件及地下水位埋深等因素进行了分区（Yukun Hu et al.，2010），确定出模拟区降水入渗补给系数分区图（图8-3），通过灌区1953～2000年降水资料及1988～2000年地下水93个长观井水位资料统计分析，确定降雨入渗补给系数（Yonxhui Yanx et al.，2006）。然后将所有面状、线状的源汇项数据分别换算成强度形式，然后通过叠加计算，再次换算成单个网格上强度，再以recharxe模块导入模型。

图8-3 泾惠渠灌区降水入渗补给系数分区 Fig.8-3 partition map of precipitation infiltration supply coefficient in Jinghui Canal Irrigation District

3.田间灌溉渗漏补给及渠系渗漏补给

田间灌溉渗漏补给及渠系渗漏补给是模拟区地下水两种主要补给源。灌区渠道分布基本覆盖整个模拟区，以面状补给来处理田间灌溉的渗漏补给及渠系渗漏补给。根据灌区士地利用统计资料和拥有的长系列灌溉用水量资料，结合不同灌溉定额和补给系数计算农田灌溉各时段的渗漏补给量，再将灌溉渗漏补给量平均分配到计算的单元格中。

4.井灌地下水开采排泄

根据对灌区地下水资源计算与评价结果，全灌区近年平均开采地下水资源量为1.2629×10⁸m³/a。灌区井网以灌溉渠系的斗、分渠为骨架，井排走向与潜水流向垂直或斜交，井距200～300m，浅型井占95%，中深井占2.4%，大口井占1.4%。开采量、开采动态等根据灌区灌溉年报中年度地下水取水情况统计表获得。灌区地下水开采量按照井流模块（Well）输入模型中。

5.蒸发排泄

MODFLOW中的EVT蒸发子程序包为线性蒸发模型，浅埋区地下水蒸发与埋深呈非线性关系，同时蒸发因素在地下水均衡分析中所占比重较大，利用线性模型计算蒸发量误差较大，本书采用阿里维扬诺夫非线性公式代替了EVT蒸发模块线性公式来计算蒸发量，用Visual Basic6.0在EVT模块中改写了源代码。阿里维扬诺夫非线性公式为

灌区农业节水对地下水空间分布影响及模拟

式中：R_ETMii，j为地下水面蒸发强度（m），取决于当地气象条件；R_ETii，j为潜水蒸发强度（m），随月份变化，用单位面积单位时间内水量体积表示；h_i，j，k为单元水头，或地下水位（m）；h_s，j，k为蒸发界面高程（m）；d_ii，j地下水极限蒸发埋深（m），与岩性特征有关；m为无量纲指数，该地区近似取2。

调整后的R_ET与调整前相比，精度有较大提高，同时用稳定流拟合效果较好时计算的地下水等水位线，与实测地下水埋深线进行叠加作为虚拟蒸发界面高程，利用虚拟蒸发界面高程代替实际的蒸发界面高程，然后将虚拟蒸发界面高程导入EVT模块中，反复调试拟合，再使得虚拟蒸发界面高程与计算等水位线的差值，与实测地下水埋深基本一致，这样避免了地下水流场拟合误差引起的实际蒸发量和模拟蒸发量在区域分布上的不一致，水资源模拟仿真度有所提高。

6.工业、生活地下水开采排泄

灌区附近周边乡镇截至目前仍没有实现自来水管网供水，乡镇企业和生活用水主要还是开采地下水，由于这些地下水井没有详细的统计资料，尤其是农村生活用水，基本上是每户都有一眼小型的抽水井，很难统计单井开采量，因此用农村人畜用水量定额的方法对其开采量进行估算，再按照面状负补给加入到模型中。

三、模型识别与检验

1.水文地质参数分区

水文地质参数分区依据灌区水文地质勘察、抽水试验资料，再结合模拟区的地形地貌、地质图、水文地质图等进行参数分区，水文地质参数（T、μ）采用分片常数法，其分区范围与形状，应符合地质条件与第四系沉积特征（图8-6，图8-7）。抽水试验所在参数分区，其参数值直接采用抽水试验求得的参数值。并以该分区参数作为基准参照参数，用推断类比法并参考其他单孔抽水试验数据，来估计其他分区的参数初值，待模型校正阶段进行确认。

2.模型识别与检验

模型识别与检验是地下水数值模拟及模型建立过程的一个关键环节。通常在模型识别与检验过程中，对水文地质概念模型重新认识，分析研究区水文地质条件，进一步对水文地质模型正确与否进行判断。模型识别与检验流程见图8-4。

图8-4 模型识别与检验流程图 Fig.8-4 Flow chart of model calibration and verification

模型识别与检验的优劣，同样也取决于建模过程中水文地质条件分析、模型概化等各个环节质量的优劣，识别与检验工作并不是一个调参的“数字”过程。模型检验与识别依据地下水模拟流场与地下水实际流场基本一致；模拟地下水的动态过程与实测的动态过程宏观相似；从水均衡的角度看，模拟的地下水均衡变化与实际要素基本相符；识别的水文地质参数基本符合实际水文地质条件。

本次共搜集了灌区石桥、泾阳、杨府、三渠、三原、西张、陂西、高陵、彭李、张卜、栎阳、新市、楼底、阎良共14个灌溉管理站93个地下水长观孔水位资料，每月定期观测6次，分别为1日、6日、11日、16日、21日、26日，记录观测地下水位动态，地下水位观测孔的分布基本覆盖整个灌区（图8-5）。

图8-5 泾惠渠灌区地下水长观孔分布图 Fig.8-5 Groundwater long-term observation hole distribution in Jinghui Canal Irrigation District

选取1996年1月至1998年12月的月平均水位观测数据用于模型的参数识别，1999年1月至2000年12月的月平均水位观测数据用于模型检验。

根据现状多年平均渠系渗漏补给量、田间灌溉渗漏补给量及地下水等水位线等信息，调整各分区导水系数T的数值进行拟合匹配（图8-6）。这一原则的实质是将灌区地下水循环看作“天然大型达西试验”，来调整导水系数T数值及分布，将产生不同形态流场，即等水位线分布，当模拟流场与实际流场宏观相似，即初步完成了对T的校正过程，同时也计算出了地下水蒸发量。校正给水度时，暂时固定T数值，调整各分区的给水度μ（图8-7）。给水度的大小影响地下水动态年变幅，通过调整使地下水动态年变幅与实际观测值接近。通过调整参数分区及分区参数值使两者之间的差值尽量最小，并据此来判断所用水文地质参数及分区是否合理。经反复调试参数，获得了较为满意的水文地质参数。部分观测孔的实测水位与计算水位拟合曲线如图8-8所示，拟合误差的绝对值大多数满足规范要求，可见各观测孔的水位拟合效果是较好的。

在模型调试过程中，充分利用水文地质勘探资料中所获得的各种信息及计算者对水文地质条件的判断。反复调试，直至流场及观测孔动态年变幅与实际观测值接近为止。使识别后的模型参数、地下水流场及地下水资源量之间达到较合理的匹配。

图8-6 泾惠渠灌区含水层导水系数T分区及补给排泄计算断面 Fig.8-6 partition map of aquifer hydraulic conductivity coefficient and recharge disrecharge sections in Jinghui Canal Irrigation District

图8-7 泾惠渠灌区含水层给水度μ分区图 Fig.8-7 partition map of aquifer specifie yield in Jinghui Canal Irrigation District

根据上述方法，利用1996年1月至1998年12月期间的地下水长观孔水位观测数据进行水文地质参数识别，用1999～2000年的地下水长观孔水位观测资料进行模型检验。模型计算求得的灌区典型观测孔水位变幅与实测水位变幅的拟合曲线如图8-9所示，灌区典型观测孔水位与实测水位拟合如图8-10所示。由图8-9可以看出，根据93个长观孔水位资料，剔除资料欠完善的部分观测孔，选择44个观测孔对其水位过程线进行拟合，统计绝对误差见表8-1。

图8-8 泾惠渠灌区地下水等水位线拟合图 Fig.8-8 Fitting chart of groundwater contour in Jinghui Canal Irrigation District

表8-1 拟合误差分布 Table8-1 Error analysis of groundwater level simulation

由图8-10及实际模拟过程可以看出，地下水位计算值和实测值相关系数（R²）在0.8～0.98之间变化，说明模型输入补给排泄要素及水位地质参数在该模拟区具有一定的代表性。模型拟合情况大致可以分为两类：一类是拟合情况比较好的，模型计算水位和实际观测水位相差较小，能够比较好地反映出该格点的水位动态趋势；另一类是模型计算水位值与实测水位值始终有一定的差异，但变化趋势基本保持一致。经分析，产生误差的，主要源于各源汇项的统计误差、地质资料的精度问题导致地层模拟误差等。还有一点需要说明的是，模型算法采用迭代求解，通过迭代法得到的解仅是差分方程的近似解，精度也受很多因素的影响。综上由模型计算流场和水位变化过程线拟合情况看，模拟计算结果比较真实地反映了灌区地下水流场特征，可以用其进行数值分析计算。

图8-9 泾惠渠灌区典型观测孔水位变幅拟合图 Fig.8-9 Fitting of typical observation hole water level amplitude in Jinghui Canal Irrigation District

图8-10 泾惠渠灌区典型观测孔水位拟合图 Fig.8-10 Fitting of typical observation hole water level in Jinghui Canal Irrigation District