中游区地下水流数值模拟

一、水文地质概念模型

（一）模拟范围

黑河流域中游地下水流场模拟范围，东起民乐总寨-山丹祁家店水库，西至嘉峪关市的河口-吕家庄（嘉峪关大断裂），夹峙于祁连山、龙首山、合黎山、北部夹山和南山之间的走廊平原，主要包括张掖盆地和酒泉东盆地两部分，面积约11300 km²。

（二）地下水系统结构概化及边界条件

研究区是一个只有侧向流入而没有侧向流出的闭型山间断陷水文地质盆地，其间填充的巨厚第四纪松散沉积物，构成连续的、统一的、相互叠置的、横向为盆地边界所限的地下水系统。周边地质构造或山体为天然地质边界。

在研究区东部张掖盆地，地下水自南东向北西运动，最终排泄于黑河干流，流出区外。在研究区西部酒泉东盆地，地下水自南西向北东运动，部分排泄于蒸发，部分排泄于黑河干流，流出区外。

区内潜水系统接受大气降水的入渗补给。在南部洪积扇地带接受河水和引灌渠水垂向淋滤渗漏补给。在中游冲积细土平原接受沟渠、水库和灌溉水入渗补给。在洪积扇扇缘和与之毗邻冲积细土平原带，地下水以泉的形式大量溢出而排泄。在冲积细土平原黑河河床地带，地下水以泉的形式集中排泄于黑河。在中游盆地北半部地下水埋藏较浅区，地下水以蒸发方式排泄。在人口稠密、农业集中开发区，人工开采是一种重要的地下水排泄方式。

在与洪积扇毗邻的地带，承压水系统得到了扇群带地下水侧向径流补给。在人口集中居住的地区和农业集中开发区，人工开采是一种重要排泄方式。

在黑河流域中游山前冲洪积扇及黑河河床地带皆为单一的潜水层，而在冲积细土平原地带则逐渐过渡为双层结构或多层结构的潜水-承压水统一体系。在张掖与平原堡之间的黑河河床在140 m之内都是透水性较强河流冲积物，切断了冲积细土平原北半部承压水之间的联系，从而使张掖与临泽形成两个各自独立的承压水地段。在具有多层结构含水层地带，由于下部的多层承压含水层之间隔水层较薄，且不连续，二者之间具有一定的水力联系，因此，可将下部承压含水层组视为统一含水系统。

在冲积细土平原区的潜水与承压水系统之间存在0.9～10.5 m的粘性土弱透水层，承压水水位一般高于潜水，但是在张掖地区的古城、碱滩附近潜水水位高于承压水水位，二者之间存在一定水位差，含水层间存在越流水量交换。

对于洪积扇扇群带、黑河河床地带等单一潜水层含水层地带，可在潜水系统中人为划定一薄层“弱透水层”，假设“弱透水层”的渗透系数很大，这样就可以把整个流域含水层结构概化为潜水-承压水双层结构。

在黑河流域中游南部的祁连山区、东部的民乐总寨-山丹祁家店水库一带，为二类透水流量边界。在西部的嘉峪关大断裂一线、南部的龙首山、合黎山和北部的北山一带，为二类弱透水流量边界。在张掖盆地与酒泉东盆地之间分界处，虽然两侧地下水运动方向不同，但是可近似为连续含水系统，只是导水性差。

经过以上分析，黑河流域中游区地下水系统被概化为具有大气降水、河渠水、灌溉水入渗补给，通过泉、河流、蒸发和人工开采排泄，为第二类边界条件，含水层间通过越流联系的双层结构拟三维流含水层系统。地下水流的数学模型可用式（9-3）表示。

（三）交换量处理

1.河流入渗补给量

主要分布于山前洪积扇地带，补给潜水系统。由于河水水位远高于地下水水位，河水以铅直下渗（淋滤渗漏）的方式大量补给地下水，入渗量的大小取决于河流流量、河床的岩性和地形地貌条件。计算入渗量时，参考了甘肃省第二水文地质队1967、1985年实测资料，并在模拟过程中给予识别和确定。

2.灌溉水入渗补给量

在黑河流域中游农业集中开发区，存在大量灌溉水的入渗，补给量大小取决于灌溉田地的岩性、潜水水位埋深和灌溉水量的大小。补给量按面状补给处理，在模拟过程中识别和确定。

3.渠系水入渗补给量

主要指黑河流域中游区内的引灌干渠渗漏补给，入渗量的大小由渠系衬砌类型和输水量的大小决定。较小渠系的入渗量，在灌溉水入渗中一起考虑。

4.大气降水、凝结水入渗补给量

入渗量的大小取决于降水量、潜水水位和地表浅层岩性等因素。据甘肃省第二水文地质队地渗仪实测资料表明，在荒漠区，只有地下水水位埋深小于2 m的地段才接受到降水入渗有效补给，而在灌区，地下水水位埋深小于5 m的地段，也观测到降水入渗有效补给。对于凝结水渗入补给，主要发生在地下水水位埋深小于2 m的地段。该项补给量按面状补给处理，由模型根据分区识别，确定各区入渗系数。

5.潜水蒸发量

潜水蒸发量大小主要与包气带岩性、地下水水位埋深、植被发育和气候干湿冷暖变化等因素有关，该量可按面状排泄处理。根据甘肃省第二水文地质队地渗仪观测结果，在黑河流域中游区潜水蒸发极限深度为6.5m，潜水年蒸发强度（q_e）与地下水水位埋深（g）之间关系曲线如图9-1所示，关系式为

西北内陆黑河流域水循环与地下水形成演化模式

6.越流量

根据弱透水层的岩性、厚度和上下层间地下水水位差条件，由模型识别，确定分区越流系数。

7.开采量

主要分布在人口集中地区和农业集中开发区，可按点井和面积井处理。

8.泉水溢出量

图9-1 黑河流域中游区潜水年蒸发强度与地下水水位埋深关系

泉水溢出量与泉水的出露高度、地下水水位差呈正比关系，主要参考1984、1986、1999和2001年泉水调查资料。

二、数值模拟与结果

（一）计算区域剖分

选用三角形网格进行剖分，剖分时充分考虑含水层岩性及富水性、河流的分布、水位观测孔位置、区域研究程度和土地开发利用状况、开采井位置等因素。在单层地下水系统中，共剖分成431个节点，747个单元（图9-2）；在双层地下水系统中，共剖分成862个节点，1494个单元，每层各有111个二类流量边界。

图9-2 黑河流域中游区地下水系统单元剖分图

在进行三角形网格剖分时，采用垂向切开的方法，即每个三角形对垂向两层都适用，包括三角形单元的编码，但是各节点的编码在两层中不同。网格剖分只对第一层潜水系统的三角形单元和节点编码。根据这个信息，计算机通过程序控制对下层承压水系统自动地进行新的一套编码。这样，对下部承压水系统数值模拟计算时，采用的是与上部潜水系统剖分完全不同的新编码。

（二）边界条件

黑河流域中游区地下水系统边界都是二类流量边界。根据调查资料、裂隙发育程度状况和现今山区降水情况，分段进行计算祁连山区基岩裂隙水的流入量。对于龙首山和合黎山山区基岩裂隙水的侧向径流流入量，按水利部门提供的径流深度，参照该区基岩裂隙垂向发育程度资料分段计算。对于夹山和南山基岩裂隙径流的流入量，参照合黎山资料进行处理。在嘉峪关大断裂，侧向流入量是参考酒泉盆地讨赖河水源地的勘察资料确定。

（三）模拟时期选择

模拟时期为1987年9月至1988年8月，共分为12个时间段，每个时段的步长均为1个月，潜水以计算区内24个观测孔的地下水水位作为模型识别的主要依据，承压水以计算区内21个观测孔的地下水水位作为模型识别依据。

（四）参数分区与模型识别

黑河流域中游模拟区的水文地质参数估值，采用分区法确定。根据区内水文地质勘探和现场抽水试验资料，通过数值模拟，将计算区的潜水系统和承压水系统参数初步分区，并假定每个分区具有相同的参数。水文地质参数主要有：①降水入渗系数以及地表水单长渗漏率；②潜水系统给水度与渗透系数；③承压水系统渗透系数与储水系数；④弱透水层越流系数。

模型识别要求符合下列标准：①地下水模拟流场特征与实际地下水流场基本一致；②地下水模拟动态变化过程与实际地下水动态变化过程基本相似（图9-3和图9-4）；③地下水模拟均衡变化与实际地下水均衡变化相符；④识别的水文地质参数符合实际水文地质条件。

图9-3 模型拟合阶段黑河流域中游区潜水系统初始流场平面图

图9-4 模型拟合阶段黑河流域中游区承压水系统初始流场平面图

模型的识别采用解逆法反求参数。即将参数的初值带入模型，应用上节所述的针对多层结构地下水系统数学模型，采用有限元法和波前法计算出潜水和承压水的地下水水位，将观测孔的计算水位与实测水位进行比较，其拟合程度采用评价函数来度量。

评价函数的表达式为

西北内陆黑河流域水循环与地下水形成演化模式

式中：M——时段总数；

N——观测孔总数；

W_j——各观测孔的权函数；

———i时段末j号观测孔的计算水头；

———i时段末j号观测孔的实测水位。

用试估校正法调整参数，并且根据不同时段对不同参数进行识别。当目标值E达到尽量小时，则认为达到了要求。识别后的分区参数如表9-1～表9-6和图9-5～图9-9所示。

表9-1 黑河流域中游区潜水系统渗透系数分区参数

表9-2 黑河流域中游区潜水系统给水度分区参数

表9-3 黑河流域中游区降水入渗系数分区参数

表9-4 黑河流域中游区承压水系统渗透系数分区参数

表9-5 黑河流域中游区承压水系统储水系数分区参数

表9-6 黑河流域中游区弱透水层越流系数分区参数

图9-5 黑河流域中游区潜水系统渗透系数分区图

图9-6 黑河流域中游区潜水系统给水度分区图

图9-7 黑河流域中游区承压水系统渗透系数分区图

表9-7 拟合阶段黑河流域中游区潜水系统水量均衡表（10⁸m³/a）

表9-7和表9-8表明，模型拟合期的潜水和承压水都是负均衡，地下水水位总体呈下降态势。由于基础资料序列连续性较差，特别是在酒泉东盆地的东部，可供利用的资料少，以至数值模拟结果存在一定误差。但是从拟合结果来看，控制区的观测孔地下水水位的计算值与观测值比较接近（图9-10），无论从整个拟合校正期的过程来看，还是从整个流场来看，拟合效果都较为理想。因此，可以认定所建的模型能够较好地反映黑河流域中游区地下水系统变化规律，服务于不同开采方案条件下地下水动态变化的预测研究。

表9-8 拟合阶段黑河流域中游区承压水系统水量均衡表（10⁸m³/a）

图9-8 黑河流域中游区承压水系统储水系数分区图

图9-9 黑河流域中游区弱透水层越流系数分区图

图9-10 黑河流域中游区部分观测孔地下水水位拟合曲线

三、地下水水位动态变化模拟

在以上研究的基础上，利用求得的水文地质参数，对流域中游未来不同开采条件下的地下水水位动态进行预测研究。预测方案有两种：（A）现时开采方案：不改变现有开采条件，地下水开采量保持在1987年9月至1988年8月期间水平；（B）规划开采方案：不改变现有开采井布局和其他水文地质条件，地下水（潜水和承压水）开采量增加1倍。预测时，以1999年9月初潜水和承压水系统等水位线作为初始流场（图9-11和图9-12），时间步长为1年，预测至2009年9月初潜水和承压水系统地下水水位。

方案A预测结果如图9-13和图9-14所示。从图9-13和图9-14可见，预测期内各时段都为负均衡。在张掖盆地东南部的民乐、六坝、毛城子和安阳一带，潜水水位基本保持不变。在张掖市和平原堡一带，潜水水位略有上升，附近河段地下水溢出量略有增加。其他地段，潜水水位都有所下降，平原堡到高台河段地下水溢出量均有所减少，其中临泽县城、临泽农场一带水位下降很大，年均降幅达0.6m。在酒泉东盆地，潜水水位呈上升趋势，年均升幅0.2m，高台-正义峡河段地下水溢出量有所增加。在张掖盆地的张掖市附近，承压水水位基本保持平稳。在临泽县城、临泽农场一带，承压水水位下降较大，年均降幅达0.9m。在酒泉东盆地酒泉市附近，承压水水位略有上升。其他地段承压水水位略有下降。

图9-11 黑河流域中游区潜水系统预测初始等水位线

图9-12 黑河流域中游区承压水系统预测初始等水位线

图9-13 方案A 黑河流域中游区潜水系统等水位线

图9-14 方案A 黑河流域中游区承压水系统等水位线

方案B预测结果，如图9-15和图9-16所示。从图9-15和图9-16可见，预测期内各时段也是负均衡，地下水水位呈总体下降趋势。除张掖盆地东南部民乐、六坝、毛城子和安阳一带潜水水位基本保持不变之外，张掖市西北部和平原堡一带的潜水水位略有上升，在附近河段河水溢出量略有增加。在其他地段，潜水水位普遍下降，平原堡到高台河段地下水溢出量均有所减少。其中在临泽县城、临泽农场一带，潜水水位下降很大，年均降幅达1.1 m。在酒泉东盆地，潜水水位持续上升，年均升幅0.2 m，在高台-正义峡河段地下水溢出量有所增加。在张掖盆地，承压水水位普遍下降。其中在张掖市附近，承压水水位年均降幅达0.8 m，在临泽县城、临泽农场一带承压水水位年降幅达1.2 m。在酒泉东盆地，除酒泉市附近承压水水位略有上升之外，其他地段的承压水水位呈下降趋势，年降幅达0.3m。

图9-15 方案B 黑河流域中游区潜水系统等水位线

图9-16 方案B 黑河流域中游区承压水系统等水位线