河套灌区井渠结合数值模拟及水资源分析预报

余乐时，朱 焱，杨金忠

(武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072)

河套灌区水资源匮乏，地下水位埋深较浅，土壤盐碱化严重，合理评价地下水资源，充分利用地下水，对于河套灌区的节水和控盐意义重大。井渠结合利用地下水与地表水共同灌溉，可减少地表用水，同时增加地下水埋深，减少蒸发损失，达到节约用水和控制土壤盐碱化的双重目的。井渠结合利用渠灌区地表水灌溉补给地下水，因此，研究井渠结合条件下，灌区地下水的采补动态过程，对维持灌区的生态平衡具有重要意义。

Visual MODFLOW(Modular Three-dimensional Finite-difference Ground-water Flow Model)是地下水数值模拟中最为普及的软件之一，国内外学者利用Visual MODFLOW研究地下水问题已取得较为丰富的成果，Zume等学者[1]利用MODFLOW分析了美国Oklahoma西北部地区抽取地下水时，河流及地下含水层的水流特征。马玉蕾[2]运用Visual MODFLOW分析了黄河三角洲浅层地下水与植被的相关关系。龚亚兵[3]建立河套盆地的地下水数值模型，分析了盆地内地下水位变化对盐碱化控制的影响。本文在全面搜集河套灌区2006-2013年的水文地质、气象、种植结构等资料的基础上，应用VisualMODFLOW软件，构建了河套灌区地下水动态数值模型，利用2006-2010年观测数据对模型主要参数进行了率定，利用2011-2013年数据对率定结果进行验证。同时，利用该验证的模型，预测了实施井渠结合灌溉方式后，区内地下水变化及水资源状况。

1 研究区概况

河套灌区位于内蒙古自治区西部，北依阴山山脉的狼山、色尔腾山、乌拉山南麓洪积扇，南临黄河，东至乌梁素海，西接乌兰布和沙漠(见图1)。地理坐标为北纬40°10′～41°20′，东经106°10′～109°30′，东西长250 km，南北宽50 km，总土地面积106.73 万hm2，现有引黄灌溉面积57.4 万hm2。整个河套灌区又分为乌兰布和、解放闸、永济、义长和乌拉特5个灌域。河套灌区属于典型的干旱—半干旱大陆性气候区，套区作物用水主要依靠引黄河水灌溉[4]。

图1 河套灌区位置

2 基于Visual MODFLOW地下水数值模型

2.1 模型概化

2.1.1 模拟范围及边界条件

本次模拟范围为整个河套灌区，利用Arc GIS在Google Earth中描绘出河套灌区范围。并划分为8.1 万个长1 000 m，宽500 m大小的网格，其中有效网格数量为3.75 万个。网格的划分是当地的水文地质条件进行划分，同时兼顾计算精度和计算成本。河套灌区南北向第一含水层厚度变化相对较大，东西方向厚度变化相对较小。因此，在划分网格时，东西方向网格长度较大，南北方向网格较小。

侧向边界条件设定为：①东边界：套区东部有乌梁素海，设为水量交换边界，套区的地下水与乌梁素海的水量交换受到乌梁素海水位的控制；②西边界：套区西部为乌兰布和沙漠，地下水水平向流动较弱，设为不透水边界；③北边界：套区北部自西向东分别为狼山和色尔腾山，乌拉山，北边界补给主要由山前侧渗补给，由于套区北部有很多抽水井，认为山前侧渗和抽水两者平衡，因此也做不透水边界处理。④南边界：黄河自河套南部流过，与南部的地下水存在密切关系，因此将其设为河流边界(三类边界)。表1展示的是黄河上离河套灌区最近的2个水文站所检测的黄河水位。模拟黄河部分时，该数据是重要的输入项。河流宽度设为350 m、河流深度平均2.5 m、河流底部厚度及渗透系数为5 m和0.5 m/d[5]。

表1 黄河水位表 m

根据河套灌区测绘资料以及地层土壤岩性特点，河套灌区内水量交换主要在第一含水层内，因此将第一含水层作为垂向模拟范围，下部边界埋深在50～280 m之间。模型垂向上分为3层，第一层为弱透水层，厚度在3～18 m之间，剩下的部分再均分为两层以便模型迭代计算。

2.1.2 源汇项

河套灌区地下水主要补给源为降雨补给和灌溉补给，两者单独计算，但是合并输入。河套灌区现有总干渠1条，干渠13条，分干渠48条，支渠372条，斗、农、毛渠8.6万条，斗渠农渠在平面上分布十分紧密。引水量较均匀的平铺在灌区面上，因此，将灌溉补给以面状补给的形式输入。由于区内渠系复杂，难以获得准确资料，且渠系所经区域也在对应的入渗分区中，不会影响总的入渗量，因此不区分渠系入渗和田间入渗，而是将两者共同计算，率定得到各分区的综合入渗系数，入渗系数不再单独进行分区。在ArcGIS中依据灌区内排水沟作为各个分区的分界，划分各渠道灌溉控制区域，如图2所示。各个灌溉控制区域的面积，渠道年均引水量如表2所示。

蒸发输入分区采用泰森多边形法则对灌区进行划分，划分结果如图3所示，各区域输入对应站点的蒸发数据。潜水蒸发系数采用解放闸沙壕渠试验结果[6]，如图4所示。

图2 降雨/灌溉入渗分区

渠名一干渠渡口渠乌拉河杨家河清惠渠黄济渠黄洋渠永济渠合济渠南边渠北边渠南三支丰济渠复兴渠义和渠通济渠长塔渠华惠渠四闸渠总计引水量/亿m35．7960．5831．9364．0020．8834．9770．1696．2841．2160．6280．0930．2884．3434．4732．7522．1256．1380．2090．73147．626面积/km21060116．8440．5673181．5893831243256175．2108138．51143827713510806158．556210139单位面积上灌溉量/m0．5470．4990．440．5950．4870．5570．2040．5060．4750．3580．0860．2080．3800．5410．3860．4170．7620．1320．1300．470

图3 蒸发输入分区

图4 潜水蒸发系数

蒸发、降雨资料采用区内的磴口、杭后、临河、五原县、乌前旗、大佘太6个水文测站的统计资料，灌溉水量则根据各渠道引水量数据。2006-2013年以来月均降雨量和月均蒸发量如图5所示。

图5 月均降雨量和蒸发量

河套灌区12月至次年5月份中旬为冻融期，水位变化机理与灌溉期不同。本文依据各分区地下水位观测数据，将相邻月份水头差值乘以由生育期和秋浇期率定获得的给水度μ，得到冻融期地层中的水量变化量，该水量变化量是降雨入渗、灌溉入渗、潜水蒸发、排水等所有源汇项的综合量，以入渗补给进行输入。

2.2 模型率定与检验

2.2.1 水位比较分析

模型利用2006-2010年河套灌区内219口观测井的水位观测数据进行参数率定，选用2011-2013年的地下水资料进行验证，以月为单位应力期，率定期共60个应力期，验证期共36个应力期。率定期和验证期灌区计算水位和实际水位的对比如图6～图8所示。

图6和图7分别为率定期和验证期各分区地下水位模拟值与实测值随时间变化对比结果，图8为率定期和验证期全灌区地下水位模拟值与实测值对比。从各个分区以及整个河套灌区的计算水位和观测水位的对比。可以看出，整个河套灌区水位的变化规律十分明显，即一年中两次上升两次下降，较好地反应了全灌区及各分区地下水位的动态变化。验证期对地下水位的模拟结果相对率定期误差较大，但全灌区的模拟结果仍然较好。

分区水位计算结果中，四闸渠、南边渠、南三支区的误差较大，主要原因有两方面：一是该部分分区主要集中在总干南部，离总干较近，用水情况比较复杂，数据上看单位面积上灌溉量较小，因此计算水位偏小。二是该部分分区面积较小，内观测井的数量不足，导致观测水位不够具有代表性，影响对比结果。在率定过程中，主要根据河套灌区200多眼观测井的总体观测数据与实测值进行率定，取总体误差最小为原则。考虑到参数合理范围，所以对该部分分区适当调大了入渗系数，无法进一步改进水位计算结果。

为评估模拟值与实测值的吻合程度，本文引入统计参数标准偏差期望SEE、均方根误差RMS、标准均方根误差NRMS以及相关系数Cor作为模型结果合理性的评价指标，如表3所示。率定期均方根误差为0.727 m，标准均方根误差为1.587%，模拟值与实测值的相关系数达到了0.99以上。验证期均方根误差为0.869 m，标准均方根误差为1.894%，模拟值与实测值的相关系数仍大于0.99，说明本文率定的模型模拟结果准确可信，可用于预测未来水位及水资源变化。

表3 地下水位计算值与观测值误差标准

本文进一步对比了灌区模拟水头分布与实测水头分布情况，结果如图9所示。

图8 全灌区率定期和验证期水位对比

2.2.2 参数结果

灌区地下水流模型中，主要参数有渗透系数、给水度、潜水蒸发系数、渠系渗漏补给系数、降雨补给系数和田间入渗补给系数。其中潜水蒸发系数和降雨入渗系数作为已知项输入，渠系补给系数和田间入渗补给系数合并为综合入渗参数进行率定。综合入渗系数率定结果如表4所示，河套灌区全区入渗系数生育期平均为0.282，秋浇期为0.345，全年的平均为0.302。给水度和释水系数如表5所示，第一层弱透水层给水度在0.02～0.03之间，弹性释水系数在0.000 001～0.000 005 m-1之间，第二、第三层给水度在0.04～0.06之间，弹性释水系数在0.000 000 5～0.000 002 m-1，该结果符合地下水文土壤参数指标[6]，结果较为可信。全区整个含水层的水平渗透系数如图10所示。

2.2.3 水均衡分析

图9 不同时期计算水头与实测水头分布对比(单位：m)

渠名一干渠大滩渠乌拉河杨家河南一支清惠渠黄济渠黄洋渠永济渠合济渠生育期0．280．430．290．310．270．270．280．280．230．25秋浇期0．330．60．350．350．350．350．330．350．330．33渠名南边渠北边渠南三支丰济渠复兴渠义和渠通济渠长塔渠华惠渠四闸渠生育期0．300．290．430．280．270．280．300．290．380．39秋浇期0．40．40．60．350．330．330．330．330．50．38

率定期和验证期各项水量的统计数据如表6所示。河套灌区地下水储量波动较小，进入灌区和排出灌区的水量大致是相等。因潜水蒸发而损失的水量率定期平均每年15.11 亿m3，验证期13.69 亿m3，是河套灌区地下水最大的输出项。灌溉入渗和降雨入渗两者作为河套灌区内地下水最大的补给来源，率定期每年向地下水补充14.20 亿m3，而验证期约13.65 亿m3。由于2012年气候干旱严重，黄河引水较少，因此验证期的入渗量和蒸发量比率定期要小，该计算结果与实际相符。黄河侧渗也是套区地下水较大的补给源，率定期平均每年0.97 亿m3左右，验证期约0.98 亿m3。乌梁素海接受河套灌区排水的同时，也直接与地下水进行水量交换，但是交换量较小。率定期的误差为2.2%，验证期为-0.9%。模型计算的各项水量总体均衡，模型准确可信。

表5 给水度及释水系数

图10 水平渗透系数(单位：m/d)

表6 水均衡分析表 亿m3

3 井渠结合模式下地下水动态预测

预测期采用井渠结合模式进行灌溉，预测时间为10 a，从2014年1月1日开始，一个月为一个应力期，共120个应力期，边界条件不变。预测井渠结合后，河套灌区内地下水变化情况。

3.1 地下水可开采区及井灌区布置

井渠结合井灌区要求地下水矿化度满足灌溉水质要求，以免地下水灌溉引起土壤盐碱化。结合地下水水质勘探资料，利用Arc GIS绘出地下水矿化度小于2.5 g/L的区域作为地下水可开采区，如图11所示[7]。本模型将典型井渠结合井灌区设置为边长为2 000 m的正方形，面积为400 hm2和边长为1 500 m的正方形，面积为225 hm2两种形式，分别布置在永济灌域和乌拉特灌域。每9个井渠结合井灌区和周围3倍于其面积的井渠结合渠灌区组成了一个井渠结合典型区，其余可开采区为井渠结合综合区。为提高典型区的计算精度，模型对该区域网格进行了加密，加密后网格大小为边长250 m的正方形。

图11 可开采区范围及井渠结合布置方式(2.5 g/L)

3.2 单位面积补给量计算

井渠结合后，非井渠结合区和井渠结合渠灌区的单位面积补给量不变。单位面积补给量主要有3类，井渠结合井灌区单位面积补给量q井、井渠结合渠灌区单位面积补给量q渠、井渠结合区综合单位面积灌溉补给量q综合。生育期与秋浇期分别进行计算，计算公式如下。

(1)

(2)

q生育期、井=q生育期、补给-q生育期、开采=

-1.5q生育期、净kα生育期、田+q降雨αp

(3)

(4)

(5)

(6)

3.3 井渠结合水文预报

井渠结合前后井灌区水位对比如图12所示，井渠结合实施后的第一年，井灌区与渠灌区内水位下降明显，随着区内地下水流场的调整，水位很快趋于稳定，达到新的平衡。井渠结合后，地下水水位振幅减小。原因是由于井渠结合前秋浇期与生育期初期，原本是水位上升时期，井渠结合后井渠结合井灌区内这两个时期的入渗水量均大幅度减小，因此上升的幅度大大减小。水位下降时期，由于入渗水量没有变化，因此降幅几乎没有变化，由此综合表现为振幅的减小。生育期间，井灌区利用地下水进行灌溉，水位开始下降，并在生育期保持较低水位，且有缓慢下降的趋势。秋浇期间，井灌区虽采用三年一灌，但是依旧补给大于损失，因此水位明显上升，并达到全年水位最大值，不过上升幅度与井渠结合前相比较小。井渠结合后的渠灌区尽管灌溉制度没有变化，但是受井灌区的影响，地下水位也有明显下降。除冻融期外，渠灌区水位始终高于井灌区，且渠灌区水位变化趋势与井渠结合前接近，生育期初与秋浇期均有明显的上升。

图12 井渠结合前后水位变化

实施井渠结合后，井渠结合典型区内地下水位分布如图13所示。冻融期由于无抽水，水头分布与从西南向东北逐渐减小，除南部区域受地形影响稍有波动外，流场整体较为均匀，无明显降深漏斗。生育期井渠结合井灌区地下水位由于地下水的开采而下降，水位低于周围的井渠结合渠灌区，形成明显的下降漏斗。受地下水整体流场影响，漏斗均不在井渠结合井灌区中心，而是稍向北或东方向偏移。周围地下水受水势影响，往井渠结合井灌区内汇流对井渠结合井灌区内地下水进行补充；秋浇期井灌区采用三年一秋浇，渠灌区一年一秋浇，灌溉水对井灌区地下水补给量的急剧减少，从而使得井渠结合井灌区地下水位也有明显下降，整体流场形态与生育期接近。

各分区的水位变化如表7所示。从计算结果可以看出，对比永济灌域和乌拉特灌域的地下水埋深变化情况可以发现，生育期乌拉特灌域的地下水降深大于永济灌域，秋浇期则相反。乌拉特单个井渠结合井灌区的面积较小，周边井渠结合渠灌区对区内的补给能力越强，井渠结合井灌区和井渠结合渠灌区的水位差距较永济灌域小。由于乌拉特灌溉面积小，距离边界较近，附近区域对地下水的补充能力有限，因此其生育期地下水反而下降更明显。井渠结合之前，乌拉特灌域秋浇期灌溉定额与永济灌域相比较小，而井渠结合后全区采用统一的灌溉定额进行计算时，其地下水降深也较小。因此，进行井渠结合区域规划时应充分考虑地域特征，并合理地选择单个井渠结合井灌区面积，使得井渠结合后，井渠结合井灌区内地下水有足够的补给来源，以避免局部区域地下水位下降过大。

图13 井渠结合典型区流场(单位：m)

灌获单个井灌区面积/hm2井灌区/m生育期秋浇期冻融期渠灌区/m生育期秋浇期冻融期井渠结合区/m生育期秋浇期冻融期永济灌域(平均降深)4000．560．890．240．340．650．230．450．770．24乌拉特灌域(平均降深)2250．640．660．380．500．640．390．570．650．39全灌区(平均降深)-------0．510．710．32

不同区域秋浇期的地下水降深均最大，冻融期降深最小。由地下水位资料表明，河套灌区生育期、秋浇期、冻融期地下水埋深分别为2.17、1.91、2.28 m，井渠结合井灌区与井渠结合渠灌区由于埋深受布置形式，地理位置等影响较大，因此取井灌区埋深均值没有太大意义，需要具体情况具体分析。但是井渠结合区的降深相对固定，可推求井渠结合后，3个时期的平均埋深变为2.68、2.82、2.50 m。

为考察井灌区与井渠结合区的水均衡状况，分别取井渠结合典型区中心的井灌区以及井灌区与对应的渠灌区作为均衡区进行分析。计算得到的年水均衡情况如表8所示，井渠结合区内地下水主要的补给来源和水量输出为入渗补给和潜水蒸发，两者分别占净补给总量和净输出总量的99.4%和96.8%。井渠结合区与区外的水量交换量较大，但净交换量很小，内部水量处于平衡状态，可见1∶3的井渠结合比是合理的。井灌区内，地下水主要输出项为潜水蒸发项，最主要来源为周围渠灌区的侧向补给。抽水与入渗合并计算，全年入渗补给量少于抽水量，因此净入渗水量为负，并成为地下水的第二大输出项。井渠结合后，侧向流动较强，井灌区周围的渠灌区地下水大量补充进入井灌区，维持了井灌区地下水的整体平衡。

表8 井渠结合区及井灌区水均衡表 万m3

河套灌区年水均衡如表9所示。可以看出，由于井渠结合后地下水位下降，年均潜水蒸发损失从井渠结合前的15.11减小至13.02 亿m3，是最主要的节水来源。黄河侧渗和湖泊侧渗量与井渠结合前相比，水量有轻微的波动，但变化量较小，并非节水对象。年均入渗补给量为从14.20 亿m3减少至12.35 亿m3，依旧是地下水最主要的补充来源。尽管入渗补给减少了，但由于蒸发量相应减少，因此总体水量依旧处于均衡状态。井渠结合后井灌区生育期采用地下水灌溉，井渠结合前对应区域均采用黄河水灌溉。井渠结合后井灌区秋浇期采用黄河水灌溉，但是频率调整为三年一灌，由此可以计算出两个时期的节水量，从而得出井渠结合后每年可节约黄河水4.09 亿m3左右，节水效果显著。

表9 河套灌区井渠结合前后水均衡表 亿m3

井渠结合后井灌区生育期采用地下水灌溉，井渠结合前对应区域均采用黄河水灌溉。井渠结合后井灌区秋浇期采用黄河水灌溉，但是频率调整为三年一灌。由此可计算出井渠结合后每年可节约黄河水4.09 亿m3左右，节水效果显著。

本研究在井渠结合区和非井渠结合区的界面处布置示踪粒子，初期和末期粒子位置如图14所示。粒子的运动结果表明，界面随时间将向非井渠结合区偏移，但移动速度很小，因此不会出现由于咸水入侵而导致淡水区变咸的情况。

图14 示踪粒子运动

4 结 语

本文利用Visual MODFLOW建立了河套灌区地下水数值模型，经过模型率定与检验，从观测井水位对比效果、流场拟合情况、土壤参数质量、水量均衡分析结果来看，所建立的河套灌区地下水流模型基本达到了精度要求，符合区 内实际的水文地质条件，较好地反映地下水系统的动态特征，可用于进行河套灌区的地下水资源评价和地下水流场演化的趋势性预测。

计算表明，井渠结合后，灌区内地下水逐渐下降，稳定后井渠结合区地下水位生育期、秋浇期以及冻融期分别平均下降0.51、0.71、0.32 m，埋深变为2.68、2.82、2.50 m。井渠结合后井灌区内地下水位降深与单个井灌区的面积大小有关，面积越大则下降幅度越大，同时与井灌区位置有关，若井灌区地下水在附近缺乏足够的水平补充来源。

采取井渠结合灌溉制度后，井渠结合区内水量基本平衡，表明井渠结合比1∶3是较为合理的。同时，在满足区内作物用水的条件下，总干渠首可减少引水约4.09 亿m3，节水作用明显。

□

[1] Zume J, Tarhule A. Simulating the impacts of groundwater pumping on stream-aquifer dynamics in semiarid northwestern Oklahoma, USA[J]. Hydrogeology Journal, 2008,16(4):797-810.

[2] 马玉蕾. 基于Visual Modflow的黄河三角洲浅层地下水位动态及其与植被关系研究[D]. 陕西杨凌：西北农林科技大学, 2014.

[3] 龚亚兵. 河套盆地地下水数值模拟及盐碱化水位控制研究[D]. 北京：中国地质大学(北京), 2015.

[4] 内蒙古自治区水文地质队. 土壤盐渍化水文地质条件及其改良途径的研究[R]. 呼和浩特：内蒙古自治区水文地质队，1982.

[5] 刘廷玺. 黄河内蒙古段干流水量损失研究(灌溉期)[R]. 呼和浩特：内蒙古自治区水利厅水资源专项, 2015.

[6] 王亚东. 河套灌区节水改造工程实施前后区域地下水位变化的分析[J]. 节水灌溉, 2002,(1):15-17.

[7] 何 彬, 赖 斌, 毛 威,等.基于GIS的河套灌区井渠结合分布区的确定方法[J]. 灌溉排水学报, 2016,35(2):7-12.

[8] 王璐瑶, 彭培艺, 郝培静,等.基于采补平衡的河套灌区井渠结合模式及节水潜力[J]. 中国农村水利水电, 2016,(8):18-24.