基于地下水数值模拟的布哈河流域地下水可开采资源量评价

彭红明，王占巍，罗银飞，袁有靖，王万平

(1.青海省环境地质勘查局，青海 西宁 810007；2.青海省环境地质重点实验室，青海 西宁 810007；3.青海省地质调查局，青海 西宁 810007；4.中国地质调查局水文地质环境地质调查中心，河北 保定 071051)

0 引 言

青海湖位于青藏高原东北部，是维系青藏高原东北部生态安全的重要水体，也是控制西部荒漠化向东部蔓延的天然屏障，同时它还是世界生物多样性保护的重要场所和环湖地区重要的水汽来源，在我国高原生态链中具有不可替代的作用[1-2]。布哈河作为青海湖最大的支流，其入湖流量占到总流量的80%，在布哈河流域内大量开采地下水势必造成河水流量的减少进而对青海湖产生一定的影响。因此，本文在论述流域水文地质概况的基础上，结合流域规划，采用数值模型的方法模拟得出合理的水资源开采量，可为当地地下水资源开发利用与环境保护提供基础的水文地质依据。

1856年达西公式的建立，开启了地下水资源的定量研究，随后斯图尔曼于1956年将数值法引入地下水资源的计算，1956年华尔顿首次利用电子计算机开展了水文地质模拟计算，之后随着计算机技术的成熟，数值模拟技术大规模地运用于地下水流动模拟和水资源计算预报。国内罗利川等[3-6]分别利用Groundwater Modeling System(GMS)、Visual Modflow等主流模拟软件在基岩山区和第四系松散堆积层中开展了裂隙岩溶水和孔隙潜水及承压水的地下水流动模拟，王平等[7-9]对天然条件下的地下水的水质演化和人为活动条件下的水污染变迁进行了模拟预测，郭帅等[10]学者借助于软件对深部地热的形成和循环径流条件进行了模拟识别，左伟等[11-13]建立地热系统的水-热模型，对不同开采或回灌条件下的地热资源量动态变化进行了模拟预测。但是目前来说公开发表的利用Groundwater Vistas[14-15]开展相关研究的工作较少。同时在布哈河流域大多数的研究成果聚焦于地表水资源，周一飞等[16-18]基于Soil and Water Integrated Model(SWIM)和Topmodel模型开展了地表水资源和地表径流的模拟和预测预报工作，在流域内所开展的地下水资源的相关研究较少。本次研究工作在开展流域水文地质调查的基础上，首先建立了研究区的水文地质概念模型和地下水流数值模型，利用其中的PEST模块，结合统测和长期观测数据对含水层参数的校准，最后预测了不同径流条件下规划水源地的开采水位降深及各项资源量变化情况。

1 研究区概况

1.1 自然地理条件

布哈河属于青海湖内陆水系，发源于祁连山脉疏勒南山，流域全长300 km，流向大致呈西北—东南方向，于刚察县汇入青海湖，全流域汇水面积约14422 km2。据河口水文站1956—2012年的长序列观测资料，年际内河水流量丰、平、枯交替出现，特枯年的年均流量为13.1 m3/s，平水年均流量为24.6 m3/s，丰水年均流量为41.2 m3/s，调查年河水流量为59.6 m3/s。

布哈河河谷属于典型的高原亚寒带草原半干旱气候，多年平均气温-1 ℃，多年平均降雨量349.5 mm，降水量多集中于5—9月，占全年90%以上。多年平均蒸发量1542.1 mm，据青海省各气象站利用的ф20与E601蒸发皿参数研究[19]结果，折算为E601的蒸发量为909.9 mm。本次开展水文地质调查的范围以布哈河河谷区为主，南北两侧略跨基岩山区，按照地貌形态及成因主要有构造剥蚀低山丘陵、侵蚀构造高山、侵蚀堆积河谷平原、剥蚀堆积山前洪积平原4种地貌类型。基岩山体由下古生界(Pz)变质砂岩、志留系(S)、二叠系(P)、三叠系(T)片岩砂岩及灰岩构成，河谷内沉积数百米的第四系沉积物。

1.2 水文地质条件

区域内地下水的分布规律和赋存条件，受到气象、水文、地形地貌、地层岩性、地质构造等诸多因素制约，具有典型的山间谷地型水文地质特征。基岩山区以层状岩类裂隙水为主，但因多分布于低山区或丘陵区，其富水性较弱。海拔3700～3800 m以上的岛状冻结区，冻结层上水发育，多年冻土(岩)为隔水层，仅在夏季消融以液态水的形式存在。布哈河及其支流沟谷沉积了较厚的第四系松散堆积物，是地下水赋存与运移的良好场所，以松散岩类孔隙潜水为主，含水层岩性以冲洪积砂砾卵石或冰水堆积的含泥质砂砾卵石为主，其厚度随沟谷的沉积环境而变化，水位埋深一般0.5～5 m不等，由河谷向两岸水位埋深逐渐增大，富水性逐渐变差。

山区基岩裂隙水及冻结层上水在接受大气降水和消冰融雪补给后，经过短暂的径流汇集，便以泉的形式排泄于沟谷中，或通过断裂以上升泉的形式溢出，汇成地表径流，是河谷地下水补给源和补给区。布哈河河谷开阔、平直，赋存丰富的地下水，主要接受河水的渗入补给。大气降水垂向渗透和两岸较大支沟的侧向补给是地下水的另一重要补给来源(图1)。据上环仓、天峻县城及天棚铁路桥三断面测流结果：上环仓至天棚西为河水入渗补给地下水，天棚西至调查区边界为地下水泄出补给河水。地下水接收补给后大致向东径流，并最终排泄于青海湖。受河谷变窄与基底隆起的影响，一部分地下水在阶地泻出形成泉，另一部分以隐伏形式排泄补给河水，此外城区的人工开采及浅埋区的蒸发排泄也是河谷地下水的排泄途径之一[20]。

图1 布哈河流域水文地质图Fig.1 Hydrogeological map of the Buha river basin

鉴于布哈河河谷第四系松散岩类孔隙潜水资源占整个地下水天然资源的绝大多数，因此本次研究以第四系地层分布的布哈河河谷作为建模范围。

河谷上游地区基本为地质资料空白区，因此将模型上游边界由水文地质调查的西边界略向下移至阳康曲汇入口，同时为确保模型预测阶段流场的变化不会波及模型边界，将模型的下游边界由调查边界向下游扩展至吉尔孟河(图2)。

图2 布哈河流域卫星影像及边界示意图Fig.2 Satellite image and boundary of the Buha river basin

2 研究区水文地质概念模型

2.1 含水层结构的概化

钻孔资料显示上部含水层为砂砾卵石层，厚度范围为2.5～102.5 m；下部含水层为泥质砂砾卵石层，厚度范围为8～150.8 m。下部含水层的渗透系数要小于上部含水层，但是二者之间无明显的隔水层，具有较好的水力联系。水平方向上河谷上游含水层的颗粒粗大，渗透系数大，富水性强。中下游河谷段地层颗粒渐细，富水性较上游稍差(图3)。因此含水层概化为单一的潜水含水层，属于二维不稳定流动模型。

图3 布哈河河谷水文地质剖面Fig.3 Hydrogeologic profile of the Buha river

2.2 边界条件

河谷两侧基岩山区地下水的富水性较差，对河谷地下水的补给作用有限，故而模型中将其概化为隔水边界，其径流量可以作为资源安全余量。

布哈河南岸的切格日唐冲洪积扇以及北岸的峻河河谷的汇水面积大，对河谷区地下水的补给量较大，在模型中将其设定为流量边界，其径流补给量根据《青藏高原柴达木盆地重点地区水文地质环境地质调查(天峻县)报告》[21]资料取值。

模型的上游边界位于快尔玛乡计算断面之上，设为流量边界，其流量由快尔玛乡断面径流量，减去该断面以上的河水入渗量。

模型下游断面设在距天棚计算断面以下较远的吉尔孟河，由于河水流量较大，且边界距离规划水源地较远，地下水开采不会影响该边界，因此模型中该边界设为定水头边界。

河谷第四系含水层以自由水面为系统的上边界，通过该边界，潜水与系统外发生垂向水量交换，如接受河水的渗漏、大气降水的入渗以及蒸发排泄。潜水含水层底部为下古生界变质砂岩，可视为隔水边界，为系统的下边界。

2.3 地下水补给和排泄项的处理

河谷区地下水主要补给项包括：河水的渗漏、上游地下水的径流流入量、降水入渗以及两侧大的沟谷中地下径流补给。主要排泄项包括：地下水开采、泉水排泄、蒸发、河流的隐伏排泄、下游断面流出的地下径流量。

河流的渗漏(排泄)：由于部分地段河水补给地下水，部分河段地下水通过河水排泄，这主要取决于河流与地下水之间的水位高程。模型中以河流模块来处理，利用河流与地下水之间的水力联系由计算单元之间的渗流来模拟：

QRIV=CRIV(HRIV-hi,j,k),hi,j,k>RBOT

(1)

QRIV=CRIV(HRIV-RBOT),hi,j,k≤RBOT

(2)

式中：QRIV为河流与地下水的交换量；hi，j，k为河流计算单元地下水水头；RBOT为河床底积物底界高程。CRIV河流含水层之间的水力传导系数，当地下水位低于RBOT时，河流渗漏量不再增加。

沟谷地下水补给：主要为南岸的切格日塘冲洪积扇和北岸的峻河河谷，在模型中以流量边界输入，根据断面径流量计算结果，其流量分别为5.32×104m3/d和5.95×104m3/d。

降水入渗：主要分布于河漫滩、Ⅰ级阶地地下水浅埋区，依据其分布范围，地下水位埋深，查取相应的降雨入渗经验系数，用入渗系数乘以研究区的平均降雨量，得到模型范围内降雨入渗量的分布，以面状补给项形式输进模型。

蒸发：主要分布在布哈河河床、沼泽及Ⅰ级阶地，蒸发量主要与潜水位的埋深、包气带岩性、地表植被和气候等因素有关，采用如下计算公式。

RET=RETM,h≥hs

(3)

RET=0,h≤hs-d

(4)

(5)

式中：RET为蒸发强度，mm/d；RETM为水面蒸发强度，mm/d；h为地下水水位标高，m；hs为地面标高，m；d为潜水蒸发极限埋深，m。

根据调查，大部分地区地表主要以粗颗粒的砂砾卵石为主，极限蒸发深度取值为4.0 m。蒸发量参数设置首先将E601蒸发量换算为大水体蒸发量，再结合相关研究[17]，进一步折算为实际蒸腾发强度，经计算后模型中设置为1434 mm/a。

地下水开采：包括自来水公司集中开采以及沿河地区零星的居民分散取水，总量约为9500 m3/d，模型中概化为排泄井。

泉水排泄：当地下水位达到地面高程时，地下水即以泉的形式排泄，在模型中以排水渠道(Drains)考虑，其计算公式如下：

Q=C(Hb-Hm),C=Kb·A/B

(6)

式中：Q为泉水排泄量，m3/d；Hb为地下水位，m；C为导水参数，m2/d；Hm为排水高程，m，即模型所在单元格的地面高程；Kb为渗透系数，m/d；A为泉水露头范围，m2；B为泉水所在地层厚度，m。

3 地下水流数值模拟模型

3.1 地下水流数学模型

根据前述的水文地质概念模型，研究区地下水流可以概化为单一潜水二维流动系统，其数学模型为：

(7)

3.2 模型求解

3.2.1 模型离散化

平面上将整个模拟区剖分为500 m×500 m的方格网，模拟面积为938 km2。由于模型上游地区钻孔和物探资料较少，建模网格放大为1000 m×1000 m。地表地形利用流域的Digital Elevation Model (分辨率23 m)来确定，垂向上根据钻孔和物探资料确定第四系基底标高，无钻孔及物探资料控制区利用软件插值确定。这样将复杂的渗流问题处理为剖分单元内简单而规则的渗流问题(图4)。

图4 模型网格剖分及边界条件Fig.4 Mesh generation and boundary condition

3.2.2 水文地质参数的初始化

主要根据地质地貌和水文地质条件，将整个模型范围划分为现代河谷区、河谷两岸阶地、支流沟谷、山前平原、山间平原，然后再从上游至下游，结合水文地质钻孔资料将每个区划分为若干个亚区，分别赋予不同的渗透系数、给水度初始值(图5)，其中河谷区河漫滩及一级阶地给水度和渗透系数最大，山前倾斜平原给水度和渗透系数最小。

图5 初始给水度参数分区Fig.5 Map showing the partition of initial water yield parameters

3.2.3 参数识别、校正

根据调查获取的河流多年平均径流量、统测点水位、地下水位长观动态、抽水试验获得的含水层导水系数与渗透系数等水文地质资料，模拟区域地下水流稳定场，以河水平均渗漏量作为稳定渗漏量，以部分统测水位数据和部分长观井4—10月的平均水位进行水位拟合，以未参与拟合过程的观测点水位埋深进行检验，采用模型自带的Parameter Estimation(PEST)模块经过反复拟合，宏观识别出含水层渗透系数和导水系数等参数的空间分布。在拟合过程中若计算水位与实际水位相差较大，则需先分析出现误差的原因，然后根据参数的变化范围进行调整，也可调整部分边界条件，直到水位拟合较好为止(图6)。拟合工作完成之后在获得最优拟合效果基础上，再利用Pilot Point软件对拟合得到的分区参数进行反演校正，最终获得模型区的渗透系数分布。

图6 模拟地下水等值线与观测等值线对比Fig.6 Comparison between the simulated contour line and the observed isoline

拟合后的水文地质参数分布规律大致符合实际(图7)，模拟得到的观测点模拟水位与实际观测水位相近(图8)，生成的初始流场与统测获得的流场比较吻合，可以较好地代表模拟区的地下水流动基本规律，同时校正后的稳定流模型中现状年与平水年两期的主要水均衡要素模拟计算值和实际调查计算值大致接近(表1)，进一步地说明所建立的水文地质模型以及赋给的参数和边界条件可以很好地代表实际情况，能够开展下一步的预测。

表1 各主要补给量和排泄量的统计结果单位(万m3/d)Table 1 Statistics of major supply and excretion (104m3/d)

图7 经Pilot Point方法识别后的渗透系数Fig.7 Permeability coefficient distribution identified by the Pilot Point method

图8 观测点观测值与模拟计算值对比Fig.8 Plot of the observed value versus the simulated value

3.3 可开采资源量预测

3.3.1 开采方案

水源位置考虑城镇发展规划、主要城镇的社会经济发展情况、地下水富水性和埋藏条件，开采技术条件等，最后确定在沿河的主要城镇快尔玛、新源镇及天棚工业发展区布设集中供水水源地。总开采量主要依据流域规划等相关资料确定，其数值要小于保障流域生态环境良性发展前提下的最大可开采量。为得到流域最大可能的地下水开采量，本次研究分别设置总开采量为30万m3/d和40万m3/d两种条件。其中方案一开采量为30万m3/d，规划在快尔玛、新源镇及天棚三处集中供水水源地，开采量分别设定为5万m3/d、10万m3/d和15万m3/d。方案二开采量为40万m3/d，在维持快尔玛5万m3/d和新源镇10万m3/d开采量不变的基础上，充分考虑天棚工业园区的发展需要，在天棚富水区另外设置三处集中式水源地开采量分别为7万m3/d、9万m3/d和9万m3/d。

3.3.2 模型预测情景

根据水文地质调查，河谷区的地下水资源大部分来自地表河水渗漏，水资源的多寡对河流径流量依赖性极高。为充分论证模拟区水资源保证程度和含水层调节能力，本次研究设置了从2014年开始连续经历30年的平均径流系列，如果开采降深小于水泵扬程，则再增加5年的枯水年进一步的验证地下水资源保障程度。

3.3.3 开采降深分析

在上述地下水流数值模型的基础上，开展预定环境下的水源地开采方案模拟预测。根据模型运算结果(表2)，水源地开采形成了区域的地下水位下降，但是形成的降落漏斗面积不大，在河流流量维持30年平均径流条件下，方案一拟设的三处水源地漏斗中心的降深分别为1.72 m、10.55 m、17.99 m，方案二拟设的三处水源地漏斗中心的降深分别为1.72 m、10.55 m、26.43 m，地下水位降深动态曲线显示基本在水源开采的前3年内即能达到稳定，反映了地下水的补给能力较强，水位稳定较快。在经历30年的平均径流后连续5年枯水期条件下拟设水源地的开采降深均有一定的增加，其中方案一拟设的三处水源地漏斗中心的降深分别为2.17 m、11.13 m、20.26 m(图9)；从开采技术条件来说，水源地水泵扬程均小于45 m，经济技术合理。方案二拟设的三处水源地漏斗中心的降深分别为2.17 m、11.13 m、43.38 m(图10)，天棚地区的开采降深增加较为显著，开采形成大于5 m的降落漏斗面积达32 km2，说明该地区的开采量大于补给能力；从开采技术条件来说天棚水源地的水泵扬程达到了48.38 m，开采井群中心最大水位降深超过了拟设水源地含水层厚度的1/2(表2)，开采经济条件不够理想。因此，认为方案二设定的40万m3/d的开采量过大，不够合理。

表2 开采条件下水源地开采技术条件汇总Table 2 Summary of technical conditions for water source exploitation under mining conditions

图9 方案一地下水降深分布图Fig.9 Distribution map of groundwater drawdown in scheme 1

图10 方案二地下水降深分布图Fig.10 Distribution map of groundwater drawdown in scheme 2

方案一条件下拟建水源地的地下水位降深动态曲线显示在三个水源地开采的前3年内即能达到稳定，反映了地下水的补给能力较强，水位稳定较快。在连续30年的平水年后连续遭遇5年的枯水年的条件下，拟设的三处水源地开采降深均有一定的增加，增加值分别为0.45 m、0.58 m和2.27 m(图11)，且在枯水期的第2年开采井内水位即能迅速恢复稳定状态，进一步证实了依赖于布哈河河谷区的地下水补给能力，在河谷区分别设置5万m3/d、10万m3/d和15万m3/d的水源地是可行的。

图11 方案一各水源地中心降深随开采时间变化曲线Fig.11 Curve of center drop of proposed water source with exploitation time in scheme 1

3.3.4 开采资源分析

鉴于方案二的开采降深较大，因此仅30万m3/d的开采资源进行分析。经计算，与多年平均径流条件下的水均衡项相比(表3)，总开采量30万m3/d时激发的河流补给量为1.8万m3/d，占总开采资源量的6%，减少的泉水及河流隐伏排泄量为20.68万m3/d，占总开采量的68.9%，由于开采导致部分地区的水位埋深增加，蒸发量减少7.47万m3/d，占总开采量的24.9%，增加的开采量主要来至于地下水泄出的减少。布哈河为青海湖最大的支流，开发利用地下水资源袭夺部分河水径流量，可能对河道的生态环境造成一定的影响。根据《黄河流域及西北诸河水资源综合规划》，青海湖水系的生态需水量约为地表多年平均径流量的84.25%，地表水的可开发利用率为多年平均径流量的15.75%[22]。按照出山口平水期平均流量计算，布哈河地表水的可开发利用量为26×104m3/d，由于目前整个布哈河流域基本处于天然状态，故此次按照30万m3/d的开采方案，虽导致地表径流量减少了22.48万m3/d，但尚在地表水可开发利用总量以内，不会影响到流域的生态环境。

表3 多年平均径流条件下方案一地下水资源均衡表(104 m3/d)Table 3 Groundwater resource balance table under the condition of annual average runoff (104 m3/d)

4 结 论

(1)在不同开采技术条件下，对拟建水源地周边水文地质条件变化的模拟结果表明，水源地开采中心水位降深能在3年内达到稳定，地下水的补给能力较强；综合考虑降落漏斗面积，抽水设备扬程等开采经济条件推荐水源地的最大开采量不宜超过30万m3/d。

(2)通过对比模拟和天然情况下地下水资源量变化情况表明，开采资源量主要来自于减少的泉水及河流隐伏排泄，其次为浅层地下水的蒸发。开采技术条件下河谷区内的河水量减少值可达74.9%，综合流域规划及河流长序列流量监测等相关资料分析认为在充分考虑生态需求的基础上，布哈河流域内地下水资源开采潜力能达到30万m3/d，该值可作为当地地下水资源可持续开发利用与环境保护的水文地质参考依据。