MODFLOW在新疆灌区地下水资源开发评价中的应用

刘 磊，董新光，吴 斌

(1.新疆农业大学 水利与土木工程学院，新疆 乌鲁木齐830052;2.新疆水利厅，新疆 乌鲁木齐830000)

新和县位于新疆维吾尔自治区天山南麓，塔里木盆地北缘，东隔渭干河与库车相望，北隔确勒塔格山与拜城县相邻，南连沙雅县，东西长136 km，南北宽91 km，全县总面积为8 823 km2，耕地面积49.35×m 104万亩。地处欧亚大陆腹地，塔里木盆地北部，远离海洋，北部和西部受天山屏障的阻隔，大西洋和印度洋暖湿气流部分可翻越帕米尔高原或天山进入本县，同时东南面又受塔克拉玛干大沙漠的干热气流影响，从而使新和县呈现暖温带大陆性干旱气候的特征。

新疆的农业属于灌溉型农业，尤其是在南疆地表水缺乏的情况下，有限的且不断减少的地表水资源已无法保证农业经济可持续发展。为改善新和县灌区耕地农业用水季节性和高峰期水资源供需矛盾，提高灌溉保证程度，保证农业经济可持续发展，科学合理的开发利用该区域地下水资源是一条有效途径。

1 研究区水文地质概况

评价区处于渭干河冲洪积倾斜平原中上部，渭干河是该冲洪积平原区第四纪松散层孔隙地下水形成的唯一来源，渭干河在出山口后形成了巨大的冲洪积扇。渭干河冲洪积扇的中上部主要为地下水的径流区，地下水径流条件由北向南，由扇顶至扇前缘逐渐变差，随着地层岩性的逐渐变细，地形坡度变缓，潜水水位埋深也逐渐变浅。在冲洪积扇的下部，由于地形坡度变缓，地层颗粒变细，从而减缓了地下水向下游排泄的速度，垂向排泄逐渐成为主要的排泄方式。

评价区地下水的补给方式主要有上游地下水侧向流入补给和区内地表水的垂向入渗补给，大气降水入渗补给鉴于区内降雨量小且次降雨量很少达到10 mm以上，在此可予以忽略。评价区地层结构简单，地层颗粒较粗大，渗透性好，地下水在接受了渭干河及上游渠系水入渗补给后，以侧向径流流入的方式补给本区。

评价区地处渭干河西岸，渠系发育，尤鲁都斯干渠、塔什艾日克干渠、依其艾日克干渠和沙雅总干渠等多条干渠均在工作区分布或穿过，同时，区内支、斗渠配套齐全，形成纵横交错的水系网，除干渠做防渗外，支、斗渠多防渗效果差或未防渗，从而使渠道渗漏成为测区地下水的另一主要补给源。此外，做为渭干河流域的古老灌区，评价区内农田分布广泛，耕作层主要为粉土和粉砂，田间灌溉水大量入渗，从而形成了本区地下水的另一重要补给源。

项目区地层岩性在314国道以北主要以砂砾石夹中粗砂、粉土、粉质粘土为主，以南逐渐变为中粗砂夹粉土、粉质粘土层，第四系松散堆积物在工作区厚度大于150 m。巨厚的松散堆积层为地下水提供了良好的赋存空间，地下水类型主要为潜水，局部存在微承压水。

在项目区东北部，含水层结构为单一的砂砾石潜水赋存区，水位埋深多为3～5 m，在县城北部和东北部大于5 m，在西南部多在1～3 m，局部地段由于受地形影响在3～5 m之间。地下水流向，总体上与地形坡降一致，即向南、西南方向径流，但局部因受引水渠及排水渠的影响而稍有改变。

地下水的径流条件主要受地形条件和地层岩性控制，测区地形平坦开阔，地势北高南低，地形坡降1.0‰ ～2.5‰。含水层岩性主要为砾石、中粗砂、细砂和粉砂，自北向南，颗粒逐渐变细，径流条件也相应逐渐变差。评价区地下水的排泄方式有潜水蒸发蒸腾排泄、地下水侧向排泄、排渠排泄和人工开采等。

2 模型的建立

2.1 模型研究范围

地下水模型的研究范围与地下水调查范围一致，在平面上呈半扇形沿东北－西南方向展布，总面积为295 km2，东北端角点坐标为(640000，4620000)，东南端角点坐标为(635500，4590000)，西北角点的坐标为(621700，4600000)(见图1)。

图1 地下水模型研究范围

2.2 含水层的概化

根据模型范围内含水层介质的成因、物质组成、空间分布规律，将模拟对象概化为非均质各项同性含水层。平面上根据含水层分布规律、结合勘探孔抽水试验资料进行了参数分区;在垂直方向大部分地段含水层变化不明显，仅作为一层，模拟深度为150 m。

2.3 数学模型

根据含水层埋藏分布、岩性及水流特征，模型区数学模型选用一层二维非稳定流数值模型，其具体数学模型表达式如下。

式中，Kxx、Kyy分别为 X、Y方向的渗透系数，L/T;H 为水头值，L;M为含水层厚度、ε为源汇项，L/T;S为给水度或弹性释水系数，潜水区取重力给水度，承压区取弹性释水系数;Ω为模拟范围;n为边界面的外法线方向;Γ为侧边界;B为底边界。

对上述地下水流数学模型，采用正交网格三维有限差分数值方法求解。为保证计算结果的收敛性，采用交替方向隐式差分格式，利用强隐式(SIP)求解差分方程。采用美国地调局(USGS)三维地下水流模拟软件包Processing MODFLOW进行数值计算。

研究区以2000年8月实测地下水位作为模拟预测模型的初始水位。利用MODFLOW自带的工具Digitizer对实测地下水等水位线图数字化，再用插值模块Field Interpolator进行插值，对每一个单元格赋初始水头值。

2.4 边界条件的处理

依据模拟区水文地质条件，平面上边界条件概化如下(见图2)，将底边界处理为隔水边界，上边界作为开放边界，考虑渠系入渗和人工开采，分别用Processing Modflow中的Recharge(面状补给)、Well(水井)程序包来处理;侧向补给边界作为定流量边界，用Well程序包来处理;侧向排泄边界用General Head Boundary(一般水头边界)程序包来处理;排水渠道用 Drain(排水)程序包来处理;蒸发用Evapotranspiration(蒸发)程序包来处理。其中AB边界为第二类补给流量边界;BC和FA边界为零流量边界;CD、DE、EF边界为第二类排泄流量边界。

2.4.1 空间与时间离散

空间上采用等间距有限差分的离散方法，在地下水模型中采用自动剖分，模拟范围内将含水层离散为1层、81行、69列，差分网格的大小为400 m×400 m，网格平均单元面积0.16 km2。模型计算区单元数(有效单元数)为1 842个，面积为 294.72 km2。

时间上以2000年1月作为初始时刻，根据地下水位的观测时间，模拟时段以月为单位，将模型调试与参数率定过程的7年分成84个时段，将2000年的初始水位运行到2007年，再以模拟好的07年的流场作为规划年的初始流场，规划年模型调试与参数率定过程的10年分为120个时段，模拟计算以天为最小单位。

图2 模型边界条件概化及规划开采井平面分布图

2.4.2 模型中数据选用

地面高程数据从1:50000比例尺的地形图上提取，然后利用MODFLOW自带的插值模块Field Interpolator插值，再赋值到各个单元格，作为地面高程。底部高程值用地面高程减去150 m。结合含水层成因、空间分布规律及勘探孔抽水试验资料进行了参数分区。依据抽水试验资料确定含水层的渗透系数和给水度。

面状补给量用Recharge程序包进行模拟，红旗农场的面状补给量为3 373×104m3/a，试验林场的面状补给量为1 685×104m3/a。侧向补给量及侧向排泄量用well程序包进行模拟，灌区边界BC、CD和DE的侧向补给量分别为316×104m3/a、773 ×104m3/a和621 ×104m3/a，灌区边界 FI、IJ和JK的侧向排泄量分别为、3 405×104m3/a、2 790×104m3/a和945×104m3/a。

地下水开采量用well程序包进行模拟。规划机井开采期每月的地下水开采量依据灌区机井开采方案中的数据进行赋值。研究区排水渠道的排水量用Drain程序包进行模拟，排水量为 2419×104m3/a。研究区蒸发量用Evapotranspiration程序包进行模拟，蒸发量为689×104m3/a。

3 模型参数率定

3.1 率定依据和方法

2007年8月实测等水位线图作为模型检验的流场图，2000年的等水位线经模型运行7年后，与2007年的等水位线拟合较好的等水位线作为规划年的初始流场。2000年到2007年评价区地下水动态观测资料作为率定含水层参数的依据。以地下水均衡计算的地下水补给量、排泄量作为模型检验的输入，以控制性勘探孔地下水位作为模型率定的检验标准。

采用自动与手动相结合的方法，通过计算水位和实际水位的拟合分析，不断地修改反复调整参数，当两者之间误差“最小时”，即认为此时的参数值代表含水层的参数。

3.2 平面流场的拟合

经调试计算，2000年的等水位线经模型运行7年后，得到2007年8月流场与实测流场比较，趋势基本一致，将2007年8月的流场作为规划年的初始流场，模拟潜水初始流场与实测流场比较见图3。计算值与实际值偏差最大为0.5 m，大部分小于0.1 m。地下水位拟合误差为0.0～0.2 m的单元数占总单元数的82%，0.2～0.5 m单元数占总单元数的18%。故认为拟合效果较好。

图3 07年8月模拟流场与实际流场拟合图

3.3 数值模拟结果与地下水均衡结果对比分析

通过数值法计算的地下水补给量和排泄量，与地下水均衡法计算的补给量和排泄量结果比较分析(见表1)，可进一步提高地下水资源的评价精度。

从全区计算结果分析:地下水均衡法计算的地下水总补给量为12 671×104m3/a，模型计算为12 511×104m3/a，相差160×104m3/a，相对误差为1.26%。地下水均衡法计算的总排泄量为13 224×104m3/a，模型计算为13 113×104m3/a，相差111×104m3/a，相对误差为0.84%。地下水均衡法计算的均衡差为－552×104m3/a，模型计算为－602×104m3/a，相差 49 ×104m3/a。

地下水模型进行参数和各项补排量调试和调整的主要依据是将2000年流场运行7年后要与2007年实测的流场基本一致，以2007年的水均衡计算的各项地下水补排量作为模型调试的初始值。在率定参数同时，调整各项补排量，使地下水模拟的初始流场与2007年的实测流场拟合最佳为止，得出评价区模型模拟各项补给和排泄量与地下水均衡法算的量拟合最佳。本次运用两种评价方法，从流场上讲，研究区地下水模拟的初始流场与2007年的实测流场基本一致;从补排关系上，水均衡法与地下水模型计算的总补给量、总排泄量和均衡结果差异不大。故认为此计算成果是可信的。

表1 地下水均衡法与数值法计算结果比较 104m3

4 应用

4.1 开采方案布置

依据《新和县布喀塔木水源地(扩建)工程报告》，水源地拟布机井50眼，计算区原水源地现状开采量为6 786.14×104m3/a，原水源地如按设计开采量开采时计算区总开采量为(6 786.14～1 880) ×104m3/a=4 906.14 ×104m3/a，计算区可开采量为6 351.81×104m3/a，则尚可开(6 351.81～4 906.14) ×104m3/a=1 445.67×104m3/a。即计算区在今后的地下水开发利用过程中，尚可开采1 445.67×104m3/a的地下水量。

设计开采量是在可开采量满足的条件下综合灌区缺水量、缺水时间及考虑本水源地置换地表水方案进行确定的。根据这50眼机井的不同的规划开采量，拟定了3种方案对评价区在项目运行5年后和10年后的水量与水位进行预测。

方案一:按照原水源地的开采井数和开采量继续运行。

方案二:根据《新和县布喀塔木水源地(扩建)工程报告》，布喀塔木水源地预计新建机井50眼，总开采量为1 000×104m3/a。原有机井的开采量保持不变，仍然为6 786.14×104m3/a，总的开采量为 7 786.14 ×104m3/a。

方案三:根据《新和县布喀塔木水源地(扩建)工程报告》，布喀塔木水源地预计新建机井50眼，总开采量为1 000×104m3/a。原有机井的开采量减小到设计开采量4 906.14×104m3/a，总的开采量为 5 906.14 ×104m3/a。

4.2 各方案模型预测

通过不同开采方案预测，方案一与方案二的水井开采量大于水源地的可开采量，运行5年后(2012年)潜水区大部分地段地下水位在现状年基础上下降0.3～2 m，集中开采区下降1～2.5 m。从方案一到方案二，随开采量增大，地下水位降幅增大，同等降幅面积在扩大。方案三开采量小于水源地的可开采量，原水源地水位在现状年基础上有所上升0.1～0.5 m，新建开采井在现状年基础上水位下降0.1～1.1 m。

各方案运行10年后(2017年)方案一与方案二的水井潜水区大部分地段地下水位在现状年基础上下降0.4～3 m，集中开采区下降2～3 m。从方案2到方案1，随开采量增大，地下水位降幅增大。方案三开采量小于水源地的可开采量，原水源地水位在现状年基础上有所上升0.2～1.0 m，新建开采井在现状年基础上水位下降0.2～1.2 m。

5 结语

对新和县老灌区通过MODFLOW软件建立的地下水数值模拟通过验证符合研究区地下水含水层的边界条件和实际真实情况，模型运行稳定可靠。利用模型对新灌区的地下水的三个开采方案进行了预测分析评价，得出按照原方案或加大开采量都超过了改研究区的可开采量，运行5年和10年后该区地下水的水量和水位都有下降趋势，故在增加开采井的情况下减少开采量是比较符合地下水资源合理利用的方案。模拟说明这一优化方案可以较好的缓解该区域地下水水资源与农业发展的矛盾，为该区域的农业经济发展和地下水开采提供了参考。

［1］Mc Donald G，Michael and Arlen W，Harbaugh.A modular threedimensional finite－difference ground－water flow model［M］.United State Government Printing Office，Washington，1998.

［2］卢文喜.地下水运动数值模拟过程中边界条件问题探讨［J］.水利学报，2003(3):33－36.

［3］吴剑锋，朱学愚.由MODFLIOW浅谈地下水流数值模拟软件的发展趋势［J］.工程勘察，2000，(2):12－15.

［4］王智，高瑾，董新光.MODFLOW在三工河流域地下水资源评价中的应用［J］.新疆农业大学报，2002，25(4):29 －34.

［5］卞玉梅，卢文喜，马洪云.Visual MODFLOW在水源地地下水数值模拟中的应用［J］.东北水利水电，2006，(3).

［6］Sun Minzheng，Chunmiao.Long－term groundwater management by a MODFLOW based dynamic optimization tool［J］.Journal of the American Water Resources Association，1999，35(1):99－111.