基于Modflow的泾惠渠地下水动态及预测研究

黄一帆，刘俊民，姜 鹏，张 斌

（西北农林科技大学 水利与建筑工程学院，陕西 杨凌712100）

地下水是贮存于地球表面松散层孔隙或基岩裂隙中的水资源，是人类不可缺少的一种宝贵的资源，对人类生活、工农业生产及城市建设都起着重要作用。同时地下水也是影响生态环境系统的一个重要因子之一，地下水的流量大小、水位高低变化往往会影响生态系统的天然平衡状态。由于含水系统的空间变异性，求取水文地质参数和查明水文地质条件较为困难，并不能用以往的确定性方法，所以我们只能尝试选择一些不确定方法。本文欲在地下水建模预测中寻求精度更高的模型，因此利用时下最流行且预测精度最高的Visual Modflow 9.3软件来进行地下水的动态预测，以求得更贴切实际情况、预测精度更高的方案。

1 泾惠渠灌区地下水概况

泾惠 渠 灌 区 位 于 108°34′34″—109°21′35″E，34°25′20″—34°41′40″N，南北宽大约70km，总面积约1 180km2，地势总趋势是自西北向东南倾斜，海拔450～350km，坡降约为3.33%～1.66%，地表及地下径流排泄条件较好，地下水潜水平均埋深为15.8 m，最大埋深40.2m，第四季河流冲击成的黄土状亚黏土灌区广为分布，土壤耕层容重、孔隙度较适宜。灌区属大陆型半干旱季风气候，多年平均降水量512.5mm，降水时空分布极不均匀，7—9月降水量约占全年的50%～60%，年蒸发量为1 140.2 mm[1]。由于灌区受到地下水排泄及补给不均衡的影响，从而研究泾惠渠灌区地下水动态的意义重大，其从总体趋势可分为两阶段：灌区建立前及灌区建立后，本文着重研究灌区建立后。

灌区建立后，除受降水影响外，重要还是受到水事活动的影响，其中特别是灌溉引用水地下水开采等活动，这些活动综合制约并影响着地下水的动态变化，说明地下水动态与人类活动息息相关，文章分别选取泾阳县代表井256＃、254＃；三原县代表井222＃、203＃、216＃；高陵县代表井306＃；临潼区代表井12＃、6＃、9＃以及阎良区代表井11＃进行绘图分析[2]（图1，2）。

图1 泾惠渠各区各代表井水位对比

图2 泾惠渠2008年代表井水位变化过程线

由图1分析看出，1984—2008年明显地出现了地下水水位下降的现象，1984年在三原县至阎良西区地下水水位较2002年下降1.47～23.43m，泾阳县处于整个地区水位较低的地带，这与泾阳县地理及开采影响不可分割，2002年在阎良与三原县交界处出现了较为明显的地下水位下降趋势，直到2008年水位沉降加深，并在泾阳县范围埋深面积扩大并且深度加深，1984—2008年的地下水水位有3.87～27.63m的下降幅度，下降幅度变大，说明地下水位有逐年下降趋势。由图2可以看出，选取的几口代表井中虽然地区不同，但是其在2008年变化趋势基本一致，可以分析出总体趋势是冬灌之后地下水位得到回升，再下降，到了4月份春灌以后又出现了一个回升，到了夏季，地下水位开始大幅度下降，而雨季的到来使得水位又开始回升。总体来看，灌溉入渗、人工开采和降雨入渗是影响该灌区地下水动态变化的最主要的三个因素，而地质因素的影响不明显。为了进一步分析泾惠渠地下水位沉降趋势，因而对灌区地下水进行软件模拟，以此提出合理的控制方案。

经过对泾惠渠灌区近年来地下水动态现状的分析，下文将对灌区深入进行预测，基于地下水模拟软件Visual Modflow结合GIS区域分析方法，通过流域内灌区及周边地下水赋存环境、运动规律、转化条件、空间布局及动态变化研究，以建立合理评价和量化地下水资源及可持续利用的方法、调控方式和指标体系，为流域水资源优化配置和综合利用提供决策支持。

2 软件介绍

Visual Modflow是目前国际上最流行、被各国同行公认的三维地下水流模拟和溶质运移模拟的标准可视化的专业软件系统[3－4]。该系统是加拿大 Waterloo Hydrogeologic Inc.通过在 Modflow软件的应用基础上，运用现代可视化技术而开发研制的，于1994年8月首次在全世界公开发行。高度集成的各类软件包，包括用于地下水动态模拟的 MODFLOW、粒子运移轨迹和传播时间模拟的 MODPATH、模拟污染物在地下水输移过程中运动的MT3D，以及用于估计水文地质参数和优化的PEST，并具有直观的、清楚的图形交互界面。新颖的菜单结构，以便于用户对研究区进行离散化和选择有效的计算单元，确定边界条件并对参数赋值，运行及对模型进行校正，还可用等值线或颜色阴影来实现结果的可视化，真正实现了人机对话。在模型的建立及结果显示过程之中，模型网格、水文地质参数和模拟结果，都可通过剖面图或平面图显示。这个软件的最大特点是把数值模拟过程中的各个步骤紧密有效地连接起来，从建模初始、输入和修改各种水文地质参数和几何参数、模型运行、反演校正各类参数，一直到最终的结果显示输出，使整个过程从头到尾整体化、规范化。目前，这些特点是我国乃至世界上的同类软件所不具备的。

3 地下水数值模拟

3.1 水文地质概念模型

根据泾惠渠灌区含水层空间结构可知，灌区含水层主要为第四系全新统Q4及上更统Q3中粗砂、砂砾石层构成的潜水含水层。中更新统Q2及下更新统Q1黏土、亚黏土因富水性差可视为相对隔水层。由于全更新统上更新统含水层厚度较大、富水性强等特征，揭露部分深度即可满足地下水开采量，灌区地下水径流受整个第四系全新统Q4及上更新统Q3含水层控制。因此，泾惠渠灌区第四系含水层可概化为单层潜水含水层，模型中地下水径流量代表第四系全更新统Q4及上更新统Q3含水层总径流量，中更新统Q2及下更新统Q1作为相对隔水底板。研究区域的含水层是由潜水含水层和浅层承压含水层组成，因为承压含水层受到人类活动的影响很小，所以，本次研究只针对潜水含水层进行模拟研究。潜水含水层的厚度和埋深等性质是随着时间和空间的变化而变化的，在本次研究区域内，潜水含水层基本上是从东到西、从两侧到中间逐渐变薄，故把潜水含水层可概化成为非均质、各向异性非稳定流的三维模拟。因此，研究区域的地下水遵循质量守恒定律、能量守恒定律和达西定律。

3.2 地下水流数学模型

根据泾惠渠灌区水文地质概念模型，将灌区浅层地下水概化为非均质、各向异性非稳定流的三维模拟单层潜水含水层，建立研究区数学模型：

式中：K——渗透系数（m／d）；h——含水层隔水底板至自由水面的距离，即含水层厚度（m）；ε——源汇项，即单位时间、单位面积垂向流入或流出含水层的水量，流入为正，表示汇，流出为负，表示源[m3／（d·m2）]；μ——潜水含水层的重力给水度；h0——水头初始值（m）；n——第二类边界的外法线方向；q——第二类边界上单位时间、单位面积上侧向流入或流出含水层的水量[m3／（d·m2））；Ω——渗流区域；S——自由面以下含水层储水系数（1／m）；P——潜水面上的降水入渗和蒸发（m／d）；Γ0——渗流区域的上边界；Γ1——渗流区域的第二类边界；Kn——边界法线方向的渗透系数（m／d）；

3.3 模型的建立

研究采用Visual Modflow 4.3求解上述数学模型，同时求解方法采用了预调共轭梯度法即PCG（Preconditioned Conjugate－Gradient）进 行 迭 代 求解[5－7]。对包括泾惠渠灌区的水文地质图进行扫描，然后采用ArcGIS 9.3进行矢量化操作，保存为Visual Modflow 4.3可导入的shp文件并打开作为数值模型的底图。根据研究区域的现实情况和本次研究精度的要求，针对图3进行单元格剖分时参考井点分布，对井点分布密集的部分进行的单元格剖分要尽量紧密些，这样考虑就相对多些。其中共有40行，53列，共计剖分了2 120个单元，其中白色单元是有效的，绿色单元是无效的。

本次研究中所应用的泾惠渠灌区的地下水观测资料、降雨、蒸发以及灌溉量和开采量等数据比较完整的是2008年和2009年，所以用2008年1月1日作为建立地下水数值模型并进行模拟的初始时间。一般在模拟地下水动态的情况下，给水度和渗透系数是两个最主要的水文地质参数，前面提到我们将泾惠渠灌区概化为非均质的，故首先要按照不同的给水度和不同的渗透系数把研究区域划分为7个子区域，见图4，每个区域的参数初值见表1。

图3 灌区网格剖分图

图4 水文地质参数分区

表1 水力传导系数和给水度初值

3.4 源汇项处理

源汇项是影响灌区地下水动态变化的重要原因，反映着整个地下水系统的变化机理，所以，源汇项对于建立地下水数值模型是十分重要的。一般来说，处理源汇项主要有两种途径，一种是以面状强度的形式进行补给或排泄，另一种是以点状的形式参与模型的补给和排泄。但总体来说，这两种处理形式在活动的单元格上是都参与的。

以2012年为例，灌区浅层地下水总补给量约为1.585 1亿m3，其中降雨入渗补给量约占19.04%；渠渗漏补给量约占23.02%；灌溉入渗补给量约占12.89%；开采回归补给量约占37.58%；侧向补给量约占7.47%。地下水总排泄量约为2.209亿m3，其中地下水开采量约占84.71%；潜水蒸发总量约占2.90%；侧向排泄量约占2.91%；人、畜、城镇工业用水量约占9.48%。补给总量与排泄总量之差为－6 238.99亿m3，其水量是动用了浅层地下水的静储备量中的一部分，即地下水疏干水量（表2）。

表2 各水平年地下水均衡量 万m3／a

3.5 模型的识别及验证

模型的识别和模型的验证通常要反复地进行模型参数的修改和源汇项的修改，才能够得到比较理想化的模拟结果，是整个建模过程中最为重要的一个步骤。本次研究对模型的识别和验证所采用的是拟合校正法，它是通过间接方法来反求模型的参数，比起直接求解模型参数具有更好的稳定性。这种方法首先要给定模型参数的一组估计值，然后建立模型并进行计算，再对比实测结果和模拟结果，算出误差，如果模型的误差精度没有达到预先设定的要求，对原来模型的参数进行适度的调整，再一次进行计算，直到模型的计算结果和实测结果的误差达到所要求的精度时，停止迭代，现在输入模型的这一组参数认为是最优的。模型的识别验证遵循下面三条原则[8－10]：一是实际的地下水流场和模型模拟的地下水流场要保持基本一致，即在模拟期的某一时刻实际的地下水等值线和模型模拟的同一时间地下水等值线形状大致相似；二是模型模拟的地下水动态和实际过程中的地下水动态大致相似，即实测的地下水位过程线和模型模拟的地下水位过程线大致一样；三是模型最后应用的参数要于十几种的水文地质条件基本符合。

从表3可以看出，所选的7口代表井中的实测值和模拟值之间的最大均方差为0.634 2m，最小均方差为0.402 1m，最大相对误差为9.7%，最小相对误差为7.2%，而相关系数则在0.695 0到0.948 7之间，实测值和模型模拟的水位值大体一样，模型的精度较为精确。而位于模型边界的区域和人类活动较多的区域，模型的模拟精度相对较差，通过分析，主要由于以下几个原因使得部分地区的模拟精度较差。

表3 2008－2009年代表井模型检验

根据前面采用的渗透系数和给水度的初值以及对边界条件的概化和对源汇项的计算，建立水文地质三维模型，并对水文地质参数进行反复修改，最终确定了符合一定的精度要求的渗透系数和给水度，见表4。

需要说明，所建立的模型无法完全准确地反映并说明地下岩层的结构，再加上模型概化的边界条件和计算的源汇项与实际情况也不可能完全一样，所以模型的拟合结果越精确，并不能说明本次模拟就具有很高的精度与可靠性，也就是说，此次模拟允许存在个别没有达到满意拟合效果的检测点的存在。

表4 水力传导系数和给水度最终值

对模型计算结果的可靠性和代表性两个方面进行评价，进行实测地下水含水层的流场和模型模拟含水层流场的比较。对地下水的流场进行对比分析[11]，采用ArcGIS 9.3绘制模拟资料代表性最强的2009年9月的模拟地下水位等值线图，如图5所示，图中面状水位为实测2009年9月地下水位分布趋势图，图上的等值线即为模拟值。由图可知，实测和模型模拟的地下水位的总趋势基本相同，但是对于富平县、高陵县与三原县交汇处这些地区模型的模拟结果和实测值的误差比较大，尤其是靠近309、305和208这几口观测井的附近区域变化很大，富平县208号井附近由于资料相对较少，因而对总体趋势的影响就相对较大，故其误差就很大。对于观测井相对密集的地区，如灌区东西部而言，误差就相对较小，由于观测井较多，对该地区的水位定位就越准确，使得每个观测井的测量值和模拟值之间的误差不是很大，再加上每个水位的基础数值很大，所以相对误差会很小，这就最终表现为模型模拟的地下水流场和实测地下水流场在图形上基本相似。

图5 潜水层地下水实测流场与模拟流场对比（2009年9月）

4 提出预测方案并检验

假定两种不同情景的源汇项方案，将给定的源汇项数值输入已经建立验证的地下水三维模型，分别预测2013年该区域的地下水动态变化。两种动态模拟情景如下：

方案1：在加大地表水的引用量并保持当前地下水开采量不变的基础上，增大干支、斗渠衬砌率，减少渠道系统渗漏量，提高渠系利用水效率，其他的参数条件不做改变，预测2014年的地下水动态。

方案2：在研究区域灌溉面积不变的情况下，2009年后按照灌溉节水规划每年增加节水灌溉面积并提高渠系利用水效率，2014年的人畜、城镇工业用水量总用水量达到2 300万m3，地下水开采量达到15 000万m3，其他参数条件按固定值输入，预测2014年的地下水动态。

按照以上假定的两种不同方案情景，将不同的源汇项数值输入到已经建立并验证的Visual Modflow模型中[12－15]，其他条件（边界条件和水文地质参数）不做改变，然后运行模型，计算出不同情景下研究区域的地下水位。

图6 两种方案对比图

由图6可以看出，两种不同方案下预测的2014年的地下水位比起2009年的地下水位（9月水位）都有所下降，但是方案一比方案二的地下水位略高，在其他源汇项和参数条件不变的前提下，方案二比方案一增加了地下水年开采量，从而导致方案二的地下水位比方案一的地下水位稍低，但是与建模的初始水位相差不大，并且变化趋势基本一致，从而进一步验证了模型的精确性。总体来说，提高农业用水效率。提高灌溉用水有效利用系数。加强农业、工业及生活用水设备等措施，有助于发展节水型农业，为以后该地区的地下水预测和规划提供一定的参考依据。

5 研究存在的问题及建议

由于本次研究水平有限，同时加上资料的完整性和时间的限制，研究中必然存在一定的问题和不足，包括下面几个方面：

（1）研究区域为泾惠渠灌区，分析了该灌区的地下水动态类型及特征，其主要受边界条件、水文地质参数和源汇项等因素的影响，研究区域越大，其关系越复杂，影响地下水动态变化的因素也越多。

（2）研究在建模过程中只选择建立了一层模型，没有考虑其他含水层（承压含水层和弱透水层），而在不同的区域潜水含水层的厚度也不一样，为了研究方便，统一采用了一个厚度而没有考虑坡度，使得地下水的排泄和补给误差变大，影响了模型的精度。虽然在此方面不可避免地存在一定误差，但在以后的研究过程中，应尽量最优地概化含水层，从而减小误差。

（3）在处理边界时，尽管将东、西、南、北各个边界做了简单的概化，但是未能反映出其真实的自然情况，再加上数据的缺乏和更新，在一定程度上增加了模型的模拟误差。

（4）在确定水文地质参数（渗透系数和给水度）的过程中，都是根据以前的经验值或试验结果，而缺少最新的资料和数据，从而在选取参数时就有一定的误差存在，为了简化研究，把含水层概化为非均质各向异性的三维流，通过把研究区域划分为9个子区域输入不同的初始参数，而没有详细地表现出其应有的特征，进一步增大了研究误差。

本文通过对泾惠渠地下水进行模拟并利用已有资料分析了当前状况，以此提出不同方案对灌区进行地下水动态预测，并从侧面讨论了灌区地下水的合理开发问题。泾惠渠灌区是一个相对缺乏水资源的灌区，由于水资源的匮乏已成为灌区工农业发展的障碍，但多年来地下水的过量开采，使得当地地下水位急剧下降，出现很多枯井，造成大面积的地下水降落漏斗，亦导致了当地水质的恶化，所以有必要对灌区进行分析评价。为提高水资源利用率，实现灌区水资源的可持续利用，应该对灌区水资源进行优化配置，加强工农业设备，有效提高地下水资源利用率，加强调控及管理制度，提出合理开发利用地下水资源的方案措施，加强灌区管理，完善水利设备，并健全法律及提高社会节水意识等都对灌区水资源合理开发利用有着重要的意义。

[1]陕西省地下水监测资料 （2008）[M].陕西：陕西省地下水管理监测局，2009.

[2]陕西省地下水监测资料 （2009）[M].陕西：陕西省地下水管理监测局，2010.

[3]白利平，王金生.GMS在临汾盆地地下水数值模拟中的应用[J].山西建筑，2004，30（16）：78－79.

[4]丁元芳，迟宝明，易树平，等.Visual MODFLOW在李官堡水源地水流模拟中的应[J].水土保持研究，2006，13（5）：99－102，105.

[5]陈南祥，姜新慧.基于GIS与层次分析法的地下水资源分区研究[J].人民黄河，2010（11）：60－62.

[6]宋丰骥，付金霞，常庆瑞.基于GIS陕西省区域经济发展环境的综合评价[J].水土保持研究，2011，18（1）：37－41.

[7]Michael H A，Mulligan A E，Harvey C F.Seasonal oscillations in water exchange between aquifers and the coastal ocean[J].Nature，2005，436（7054）：1145－1148.

[8]郭晓东，田辉，张梅桂，等.我国地下水数值模拟软件应用进展[J].地下水，2010，32（4）：5－7.

[9]Brodie R S.Integrating GIS and RDBMS technologies during construction of a regional groundwater model[J].Environmental Modelling ＆Software，1998，14（2）：119－128.

[10]伊燕平，卢文喜，张耘，等.基于径向基函数神经网络的地下水数值模拟模型的替代模型研究[J].水土保持研究，2012，19（04）：265－269.

[11]叶剑锋，刘小勇.基于GIS技术的地下水管理系统研究[J].水土保持研究，2011，18（3）：247－251.

[12]El－Kadi A I，Oloufa A A，Eltahan A A，et al.Use of a Geographic Information System in Site－Specific Ground－Water Modelinga[J].Groundwater，1994，32（4）：617－625.

[13]郝治福.石羊河流域邓马营湖区地下水位动态变化特征及数值模拟与预报[D].北京：中国农业大学，2006.

[14]Kollarits S，Kuschnig G，Veselic M，et al.Decisionsupport systems for groundwater protection：innovative tools for resource management[J].Environmental Geology，2006，49（6）：840－848.

[15]李伟业，付强.基于GIS的三江平原地下水资源评价[J].水土保持研究，2007，14（4）：92－95.