地下水数值模拟
[摘要]地下水数值模拟作为预测、评价地下水资源的工具在解决具体的水文地质问题过程中起着日益重要的作用。本文通过地下水模拟软件中具有代表性的GMS，阐述了它的一些基本特点，以及它在地下水模拟过程中的建模方法，并给出了一个计算实例。
[关键词]地下水；数值模拟；GMS 1 地下水数值模拟的一般步骤 1．1 水文地质条件分析
研究和了解计算区域的地质和水文地质条件，是运用数值法进行地下水资源评价的基础。根据评价区的地质、水文地质条件、评价的任务及取水工程的类型、布局等。合理地确定计算区域以及边界的位置和性质。此外，对区域水文地质条件的了解，还有助于下一步进行模型识别。为此，应查明含水介质条件、水的流动条件及边界条件等三方面。
(1)查明含水层在空间上的分布情况；含水介质厚度；含水介质透水性、储水性变化情况，做出含水层非均质分区图，即根据渗透系数K和给水度进行分区；查明主含水层与其它含水层的水力联系。对于条件复杂的地区，应进行适当的概化。
(2)查明是承压水还是无压水；是层流还是紊流；地下水流是一维，二维，还是三维。(3)区域边界定义了计算区域的范围，而边界条件的给定对于地下水资源的评价结果有着较大的影响。因而查明边界空问分布形状以及边界的性质，给出边界值，是运用数值法进行地下水资源评价的重要工作。当边界条件复杂，要给出定量数据有团难时．应通过专门的抽水试验来了解，也可以留待识别模型时来验证或修正边界条件。在应用数值法计算之前，要用均衡法对全区进行均衡计算。这样可以在总体上把握地下水的均衡情况，使数值计算结果更趋合理。把地下水的各均衡项分配到各抽水时期和各剖分单元或节点上。在均衡分析中，要特别注意与地下水位有关的均衡量的确定，如降雨人渗量、蒸发量、越流量等。 1．2 建立水文地质概念模型和数学模型
实际的水文地质条件是十分复杂的，要想完善地建立描述计算区地下水系统的数值模型是困难的。因此，应根据水文地质条件和地下水资源评价的目的．对实际的水文地质条件进行简化。这一过程称为水文地质条件的概化，其原则为：根据评价的目的和要求，所概化的水文地质概念模型应反映地下水系统的主要特征；概念模型要简单明了：概念模型要能够被用于进一步的定量描述．以便于建立描述符合研究区地下水运动规律的微分方程的定解问题。
水文地质条件的概化通常包含以下几个方面：计算区域几何形状的概化；含水性质的概化；边界性质的概化；参数性质(均质或非均质，各向同性或各向异性)的概化；地下水流状态(一维、二维或三维)的概化。
对计算区域进行剖分，是数值法的重要工作之一。对不同的计算方法，其剖分形式各不相同。剖分时应考虑各种分区界限，如参数分区。行政分区，地表水体，断层等。以便提高计算精度。剖分的疏密程度好要考虑以下因素：在重点评价区和重要开采地段应加密剖分单元；在地下水水位变化较大地区应适当加密：在水文地质条件变化较大地区也进行加密。此外，剖分时要尽量将主要开采井和拟合水位用的观测井放到节点上。 1．3 确定模拟期和预报期
根据资料情况和评价的要求确定模拟期和预测期。模拟期主要用来识别水文地质条件和计算地下水补给量，而预测期则用于评价地下水可开采量和预测一定开采条件下的地下水位。对于地下水水量评价。一般取一个水文年或若干水文年作为模拟期。这样可最大限度地避免前期水文因素对地下水系统的影响。预测期的确定主要取决于评价的目的和要求。 1．4 水文地质条件识别
为了验证所建立的数值模型是否符合实际，还要根据抽水试验的水位动态来检验其是否正确。即在给定参数、各补排量和边界、初始条件下．通过比较计算水位与实测水位，验证模型的正确性。识别既可以对水文地质参数进行识别，也可以对水文地质边界性质、含水层结构作进一步的确认。识别的准则为：计算的地下水流场应与实际地下水流场基本一致，即两者的地下水位等值线应基本吻合：模拟期计算的地下水位变化趋势应与实际变化趋势一致，即两者的地下水位动态过程基本吻合；实际地下水补排差应接近于计算的含水层储存量的变化值：识别后的水文地质参数、含水层结构和边界条件符合实际水文地质条件。满足以上准则，则可认为数值模型反映了计算区域的地下水流动规律，可用于地下水资源的评价和预报。反之，则应对水文地质概念模型进行适当修改，以达到上述要求。
检验数值模拟结果的可靠性通常是比较计算结果与实测数据(例如地下水水位)的拟合程度。通过调整参数的分区可以得到理想的拟合结果。但是单纯的追求拟合效果并不代表数值模拟结果的可靠性提高了。应当客观地认识地下水数值模拟技术的定位：地下水数值模型是定量刻画地下水系统的工具，模拟结果的好坏除了模型方法本身之外，主要取决于水文地质概念模型的合理概化(含水层结构、边界条件、补给与排泄条件等)。二者应当有机的结合在一起．才能有效的发挥地下水数值模拟技术在水文地质研究和调查中的作用。在实际工作中。对水文地质条件的认识是一个逐步提高的过程，对地下水数值模拟技术的应用应该力争达到概念模型基本正确和模型手段较合适，只有这样数值模拟的结果才具有 较高的可信度。
1．5 地下水资源评价和水位预报
经过验证的模型只能说是符合勘探试验阶段实际情况的模型，用来进行开采动态预报开寸，还应考虑开采条件下可能出现的变化。如边界条件和地下水的补给排泄条件可能会随开采发生变化。根据开采条件对模型进行修改以后．便可用来正演计算，可以解决以下问题：预报一定开采条件下。水位降深的空间分布和随时间的变化情况：计算一定期限内水位降深不超过某一限度时的可开采量；确定某些水均衡要素，计算补给量或补给资源量。求出稳定开采条件下的可开采量；进行不同开采方案的比较，选择最佳的开采方案；研究地表水和地下水的统一调度、综合利用．进行水资源综合管理。 2 基本数学模型及其解法
对承压含水和非承压含水的水流运动．以水头作为基本变量，其满足三维的有源二阶扩散型偏微分方程；对定常流动或对承压含水流动中忽略弹性释水则蜕变为二阶泊桑Poisson方程。基本方程如下：
式中 ,Kxx,Kyy和Kzz分别为对称二阶渗透系数张量K沿 x，y和z方向的主轴数值：H一地下水头，由计算所得的h值，利用Darcy公式可计算出地下水流量，利用含水层厚度及给水度或储水系数可计算出地下水储量；W一源汇项，水井、降雨人渗、蒸腾蒸发等均归人源项；Ss一对非承压含水层为给水度，对承压含水层为储水系数；T一时间。
基本方程(1)的有限差分方程可表示为：
式中， 为第m时段在单元i，j，k的水头：CV、
CR和CC为结点i,j，k和相邻结点间的水力传导系数， 为表示源汇项水头的系数，
为源汇项水量，
表示给水度，DELRj表示所有行中第j列的单元宽度，
表示单元i,j，k的垂直厚
DELCi表示所有列中第i行的单元宽度，度，
表示时间。
定解条件：对差分方程(2)，需给定初始条件和边界条件以及源项常是时间t的函数)。
初始条件：在计算域内给定
(通
，其中i,j,k为计算域 内网格单元编号；
边界条件：在已知边界上给定，其中i,j,k为已知边界 上网格单元
编号，m为计算H寸段内时间分段序号。
差分方程(2)为对角稀疏矩阵，通常用迭代 计算法求解，如强隐式、分片松弛法等。 3 GMS的基本功能
GMS是Groundwater Modeling System的英文字头缩写，由图形用户界面模块和地下水数值计算模块耦合而成，形成模型前处理(模型构造与初始及边界条件的给定)，数值计算模型(如MODFLOW、FEMWATER等)及后处理(模型和计算结果图形图像输出)结构。它的主要功能包括：水文地质数据可视化建模。如调用地理信息系统构造分析区域的地图坐标以确定计算区域、钻孔、机井、河流、湖泊等的地理位置；通过钻孔数据，建立三维地层剖面和实体；将实体的水文地质参数及所构造的边界条件调入数值分析模型。模型校验与计算。通过页面选项，数据文件输入输出方式选择计算包、计算时段、给定初始条件等：对所建模型的给定条件进行检验、核查，确定计算条件完备后实施数值计算。计算结果的图形图像显示。如用等水位线、色彩层图、流速矢量图及图像动画演示等方法显示所得出的地下水流动模拟结果。
GMS功能强大，包括所有与地下水流动有关的数值模拟问题，这里仅论及GMS5．0版中与地下水数值模型MODFLOW 有关的GUI模块的基本功能，即：
GIS(地理信息系统)：将shp文件转化为MAP文件。确定计算区域中的点(机井、监测井等)、线(河流、计算区域边界)、面(湖泊、降水、蒸发区域)的坐标值。
MAP(地图)：设定计算区域；给定点、线、面上的已知条件；将所建模型水文参数与边界条件，即方程(2)种的网格参数与源汇项值Pijk与Qijk，映射到MODFLOW 中。
Borehole(钻孔)：显示钻孔地层柱状图，生成地层剖面图和三维实体图(SOLID)；钻孔潜水面值可转换为TIN(三角形不规则网格)表示的曲面或二维散点(SCATTEK)，作为计算初始条件调入MODFLOW ：Borehole模块中的Horizons程序可生成HuF(Hydrogeologic Unit Flow)包，将地层参数转换到MODFLOW中，确定方程(2)中的计算系数。
TIN(三角形不规则网格)：该模块用来形成三角形不规则网格，用以表示空间曲面，其上的曲线与曲面可转化为MAP模块中的弧(Arc)与多边形(Polygon)；从Borehole到Solid的转换过程中，TIN网格点是插值坐标点。
Solid(实体)：该模块用来表示三维地层模型，可任意切割剖面，产生逼真的图像；三维地层参数与空间坐标可映射到MODFLOW 上，确定方程(2)中的计算系数。
Scatter(散点)：二维与三维Scatter主要功能是捅值，完成TIN、Grid间的数据传输。
Grid (差分网格)：三维Grid是计算模块，MODFLOW 的所有计算功能都可在该模块中独立实现(网格建模法)。当用前述模块构建模型时，在该模块中需选择计算包，给定时间步段、迭代计算方法，即给定初始条件。计算结果的输出显示也在该模块中。
MODFLOW2000计算框图包含有四个进程，即地下水流动分析，敏感性分析，观测井比较与参数估算进程。对方程(2)的迭代求解，提供了四种算法程序包，
即强隐式算法程序包（Strongly Implicit Procedure)，分片连续超松弛算法包 (Slice-Successive Over—Relaxation)，预置条件共轭梯度算法包
(Preconditioned Conjugate)直接求解包(Direct Solution Package)。 4 地下水模型构建与计算实例
针对北京某镇的地层、河流及机井分布情况， 上述计算实践。我们认为需在以下几个方面做构建GMS地下水数值模拟模型。建模过程如下：
(1)利用GIS模块将shp文件转化到MAP模块中。用以确定计算区域。 (2)根据机井钻探资料编制钻孔(Borehole)数据库文件。由MAP生成的机井坐标确定每一钻孔在计算区域的坐标位置。生成Borehole柱状图。
(3)在MAP模块下，给定井的抽水量、河流、边界水位变化、回补与蒸腾蒸发变化情况。将MAP文件中的数据转化到Modflow中．该步实现了有限差分方程(2)中的源汇项的转换。
(4)用Borehole中的Horizons方法生成Solids地层图。将Solids转化到Modflow中，该步实现了有限差分方程(2)中的计算系数转换。准确的地层图是保证数值计算正确的关键。Solids中插值函数的选用对生成的地层也有影响．而地层分布Horizon中id值的标定直接决定地层分布。
(5)在三维Grid下设置Modflow的初始条件，选择计算包及迭代计算方法。 (6)检查所建模型及定解条件的给定是否正确。运行Modflow进行计算。 (7)显示计算结果。因受限于计算机的储存量及计算速度。本文在计算时对地层柱状图进行简化。将土壤参数相近的层合并，将18个土壤层简化为6层。对面积为6O．4 km2的计算区域，用边长约为9 kmX9．5 km的长方形覆盖，分为6万个计算单元(1OOx100x6)。计算时有9眼井抽水，河流给定水头边界，考虑了降水回补与蒸腾蒸发的影响。
按一年计算，分为12个计算时段，井抽水量最大值发生在3-5月份，降水回补量最大值发生在6~8月份最高。
计算结果正确反映了水位变化趋势。2月份无降雨，抽水量不大，地下水位缓慢下降；5月份抽水量大，降水回补不大，地下水位迅速降低；8月份抽参数估算进程。对方程(2)的迭代求解，提供了四 量不大，降雨回补量大，地下水位回升；12月月份降雨、抽水量都不大，地下水无明显变化，但因河流水位降低，在靠近河流区域水位降低。 5 结论
GMS用于实际地区的地下水分析计算，除要深入学习地下水动力学理论、了解GMS软件结构及各模块的功能与关系外，正确地构建土壤地层模型。给定合适的边界条件尤为重要。通过上述计算实践，我们认为需在一下几个方面做好工作：
(1)加强地下水动态的监测工作，收集尽量详实的数据。
(2)正确构造计算区域地层模型。考虑含水层与参数的不确定性[4]。GMS提供了多种生成地层信息的方法，如Solids，Horizon，T—Progs和多种插值方法．需对各种方法进行适用条件与结果的比较分析，以选用最为合适的方法构建模型。
(3)单元剖分问题。若选用三维地下水模型：垂向上剖分层不能过少。否则无法描述地下水在垂向上的运动。达不到三维模型的目的；平面上不能过大，否则技术单元像纸一样薄，会引起数值技术本身的问题。
(4)对一眼井，其抽水量来自不同含水层时，如何确定每层的抽水量将对各层的地下水位分布有显著影响。这也是地下水开采管理需解决的问题。
(5)提高地下水数值模拟技术的硬件支撑水平。
(6)将地下水资源评价与规划相结合。
[参考文献]
[1]陈雨孙．地下水运动与资源评价[M]．北京：中国建筑工业出版社．1986． [2]陈崇希，唐仲华．地下水流动问题数值方法[M]．武汉：中国地质大学出版社，1990．
[3]薛禹群主编．地下水动力学原理(第二版)[M]．北京：地质出版社．1997． [4]李竞生，姚磊华．含水层参数识别方法[M].北京：地质出版社
