海河流域地下水开采初期数值模拟及水量平衡分析（张晓明,薛丽娟,张奇,李建新）

关键词:地下水模型;MODFLOW;开采初期状态;海河流域

中图分类号: TV213. 1 　　　　　　　文献标识码: A

海河流域平原位于中国北方华北地区,研究区总面积13. 1万km2。自上世纪70年代,海河流域人口和社会经济规模迅猛发展,加之连年干旱少雨,流域内浅层地下水资源已经成为社会发展和经济建设的主要供水水源[ 1 ] ,因此海河流域地下水超采日益严重。2000年海河平原浅层淡水总超采量为45亿m3 ,超采总面积已达59550km2 ,占平原区中有淡水区面积的63. 1%。由于长期超采地下水,已经导致唐山北部、北京西南、河北太行山前平原部分地区含水层被疏干[ 2 ] 。河北太行山前的浅层水漏斗已经发展成为彼此相连的复合漏斗群,成为我国浅层地下水第一大漏斗。海河流域地下水可开采利用的潜力已经十分有限,可利用前景暗淡[ 3 ] 。

地下水的过度开采已带来一系列生态环境问题,如地面沉降、湿地湖泊萎缩、河流干涸、海水入侵、咸水位下移等等[ 4 ] 。加强地下水管理,是目前海河流域水利所面临的一项重要工作。国内外很多学者专家提出了恢复地下水的各种措施,包括人工回灌、跨流域调水等,但是这些工程措施并不能完全解决现存的严重问题[ 5 ] 。海河流域地下水资源管理的核心是保证地下水资源长期持续利用,地下水模型系统的建设将是管好地下水的关键技术之一。模型可以模拟地下水漏斗的发展过程,并能预测未来年份地下水的开采程度,从而为水资源规划提供了保障。地下水模型在水资源管理中作用已得到了充分的认可[ 6 - 8 ] 。

在海河流域平原区,由于对初始水文地质环境的重视程度不够,且70年代以前的地下水补给、排泄、水位等方面的观测资料较少,以往的地下水研究大多以70年代后期地下水环境为研究对象,缺少对开采初期水文地质环境的探求。因此,大多数模型工作者以70年代后期地下水超采期的某一时刻为初始条件,建立了非稳定流模型。这个初始解不是真正的稳定条件下的初始值,所建的非稳定流模型也不能反映地下水环境变迁的全貌。另外,地质结构的研究一直是开展地下水研究的难点之一,但是由于本地区面积辽阔,地层分布差异性较大,尚没有一个能全面反映该区域地质概况的地层模型。本文采用美国地质勘查局MODFLOW模型,以地下水模拟系统GMS[ 9 ]为平台,较为理想地解决了地层模型建设的难题,并建立了地下水稳定流模型,模拟海河流域平原区开采初期的水文地质条件。研究目的主要有,全面认知本地区水文地质条件和初始的水量均衡场,在计算机上模拟区域地层结构和水量均衡关系,为非稳定流模拟提供初始水头条件、地质参数和其它计算条件,为该地区的水生态环境修复和水资源规划提供参考数据。

1 浅层地下水系统开采初期状态

海河流域平原区按成因类型可划分为三个地貌单元:山前倾斜平原区、中部冲积平原区和滨海平原区。除个别地点有基岩出露外,绝大部分为第四纪松散物质覆盖。华北平原第四系是一套几何形态复杂的多种沉积类型交叉叠置的含水岩系,含水层岩性、结构、厚度等具有水平变化规律。在山前平原含水层呈扇状结构,扇轴含水层岩性以砾石卵石为主,厚度大;扇间含水层粒度变细,厚度变薄。在中部平原含水层逐渐过渡为湖相沉积穿插河流沉积的舌状结构,含水层岩性以中细砂为主,厚度在靠山前平原方向变薄,向滨海方向又略变厚。向东部、南部的滨海平原含水层又过渡为湖积的岛状结构,含水层岩性以粉细砂为主,厚度又变薄[ 10 ] 。

开采初期的浅层地下水流场与地表地形坡降相一致,地下水位等值线与地表等高线基本相似。地下水的流场也是从北、西、南三面汇入渤海湾。以1964 年的地下水埋深为参照,分析得出各个分区的原始地下水位如下:滏阳河- 子牙河一线以北的北部平原,降雨相对于南部平原少,河网不如南部平原密集,河流多为季节性河流,不能常年保持河水位,因此北部平原的地下水位比南部平原低。在山前区,地下水埋深大约为2m (燕山平原) 、3m (太行山前平原) ;中部平原地下水位埋深大约为1. 5m;滨海平原小于1m。滏阳河- 子牙河一线以南的南部平原,雨量相对充沛,河网纵横,地表水充沛,地势较为平缓、退水慢,因此地下水位普遍高于北部平原。山前区地下水位埋深大约1. 5 - 2m;徒骇马颊平原、黑龙港平原埋深小于1. 5m,局部地区多有泉水,盐碱地发育。平原区地下水在大量开采以前,第一含水层组和全淡水区的第二含水层组的水位多年动态变化与降水过程具有同步变化特点,属于降水入渗补给～蒸发径流排泄型。浅层地下水的补给来源主要是降水入渗补给,其次是山前侧向补给、河道渗漏补给以及黄河沿岸侧渗,排泄主要是潜水蒸发、向河流排泄以及排泄入海[ 11, 12 ] 。

2 地下水数值模型

2. 1 水文地质概念模型

本文采用的地质钻孔勘探井共计292 眼, 井深度150m～720m,主要描述了岩性、水位、单位涌水量、渗透系数、水温及水化学特征。这些资料来自于水利部天津勘测设计院地质勘探队和河北省地质局。钻井分布(图1) 。应用ACCESS软件建立了钻孔资料数据库[ 13 ] ,并在GMS中建立地层模型,根据地下水模型的目的和GMS对数据的要求,对钻孔资料作了必要的概化: 1)井深裁剪。

本次研究以埋深120m以上的浅层地下水为研究对象,此深度外的勘探资料不作分析; 2)岩性分类简化。原始的近300眼的勘探资料描述了40多类岩性。有些岩性水文特征相近,也有些描述与本次地下水模拟关系不大。在不影响地下水模拟前提下,对岩性作了概化处理,最后归纳出6类岩性,即:表土、细砂、中粗砂、卵砾石、粘土、基岩,用于地层模型的建立。对某些缺乏钻孔的区域,依据地质剖面资料[ 14 ]人为补充一些勘探井资料。GMS建立的地层模型(图2) 。典型的地质剖面图(图3, 4) (剖面位置在图1中作了标注) 。

2. 2 数值模型

将平原区地下水边界分为三类: 1)侧渗边界:包括山前侧渗边界(起始于秦皇岛,终止于河南焦作以南)和黄河内滩侧渗边界(起于河南武涉,向东止于山东东阿) ,在模型中定义为通用水头边界,模型根据边界内外的水头差,计算越过边界的水量。2)海岸线边界:即东部渤海海岸线边界,在模型中定义为定水头边界,赋予平均海平面高程。3)黄河边界:即南部黄河边界(起始于山东东阿,终止于渤海湾) ,在模型中定义为河流边界,模拟黄河与地下水之间的水交换。

由于研究区面积较大,考虑模型的复杂程度以及计算机运行的速度,将模拟区剖分为96行, 92列, 20层,共计176640个单元,网格尺寸为7000m ×7000m,根据单元所处的位置和所包含的岩性, GMS计算各单元的岩性物理参数值。本文以1965年作为海河流域开采初期状态模拟的年份,文中没有收集到的数据均采用多年平均值代替。

2. 2. 1 降雨入渗补给计算

降雨是地下水补给的主要水源,在海河流域平原区13万km2 区域内均匀分布了141个具有代表性的雨量站,每站控制面积平均为921km2 ,基本控制了该区域降雨资源。由于没有收集到1965年的降雨资料,模型采用多年平均降雨量进行计算,该值代表了较长年份系列的多年平均值。将雨量插值到离散的网格单元中,并与该单元的降雨入渗系数的乘积作为该单元的降雨入渗量。在海河流域平原区接受降雨入渗补给的主要岩性有:砾石、粗砂、中砂、细砂等主要分布在北京、保定、石家庄等山前区;亚砂土、亚粘土主要分布在中部平原,淤泥质亚粘土、粘土主要分布在滨海区。以野外土工试验给出的不同地区不同岩性降雨入渗值[ 15 ]为参考,在地层模型相应岩性和对应分区中进行赋值。经计算,海河流域平原区年降水入渗补给量137亿m3 / a,该结果与海河流域水资源综合规划得出的多年平均值136亿m3 / a基本吻合。

2. 2. 2 潜水蒸散发计算

通过查阅1964和1965年水文年鉴,得到了各省市水位观测站的年均蒸发量观测值。研究区蒸发量计算主要用GMS中蒸发软件包通过地表高程,水位埋深以及蒸发极限埋深计算所得。

当计算潜水水位等于或者高于地表时[ 16 ] : ET Rate =Max ET Rate (1)

潜水水位埋深等于或者是低于极限埋深时: ET Rate = 0 (2)

潜水水位处于二者中间时,按埋深线性变化。

根据海河流域平原区的特点,以地下水自由面5m埋深作为潜水蒸发极限埋深,即自由面大于该埋深的区域不发生蒸散发。最大蒸发率(自由面位于地面时的蒸发率)以水文年鉴和水资源规划成果为依据。

2. 2. 3 山前侧渗量计算

山前侧渗指山区地下径流补给平原区的地下水量。山前侧向补给量及下游侧向流出量计算公式为[ 16 ] :

(3)

式中:Q侧入- 计算时段侧向流入量(亿m3 ) ; Q侧出- 计算时段侧向流出量(亿m3 ) ; K - 含水层加权平均渗透系数(m /d) ; I - 计算断面水力坡度( ‰) ; L - 计算断面长度(m) ; H - 含水层平均厚度(m) ; T - 计算时段长度( d) 。

模型直接采用了水资源规划多年平均研究成果,分段赋值到模型中。山前多年平均侧渗总量约17亿m3。

2. 2. 4 河流渗流补给量

依据收集的1965年平原区主要河流汛期的月平均观测水位值可知,平原区内开采初期的地下水位大都高于河流水位,河流附近渗透系数较大,有利于地下水排泄。采用MODFLOW河流计算子模块, 河流与含水层交换水量计算公式为[ 16 ] :

Qriv =Criv ( hr - h) (4)

式中,Qriv为河流与含水层之间交换流量; hr 为河流水位; h为河流附近网格的地下水位; Criv为河流与含水层相互作用的水力传导系数:

Criv = KrivLrivW /M (5)

式中, Kriv为河底渗透系数; Lriv为网格单元内的河道长度;W为网格单元内河道宽度;M为河底沉积物的厚度。

2. 2. 5 模型识别与校正

模型识别主要遵循以下原则: 1)模拟的地下水流向与实际观测的地下水流向基本一致; 2)模拟等水头线与实测值相吻合; 3)率定的水文地质参数在合理范围内。根据以上原则,采用试估- 校正法,对块状结构的华北平原地下水数值模拟模型中不同的岩相赋与不同的水文地质参数,利用1959年地下水观测数据(由于没有收集到1965年的地下水观测数据,而1959年与1965年均处于开采状态初期,地下水一般接受降雨补给,流动较平缓,所以地下水自由面等值线比较接近,没有较大的变化,因此选择1959年的地下水等值线作为模拟依据) ,对研究区地下水系统进行了识别,模拟期内潜水拟合结果(图5) ,各种岩性的水文地质参数(表1) 。

3 水量平衡分析

表2和表3分别列出了研究区稳定流条件下地下水总的水量平衡和三级区水量平衡分量。从表2中可知,研究区地下水系统总补给量为156. 6亿m3 ,总的排泄量为156亿m3 , ,拟合均衡的效果很好。其中,降雨补给量为136亿m3 / a,占总补给量的86. 8%;河道渗漏补给量为3. 1亿m3 / a,占总补给量的2. 0%;山前侧渗量为17. 5亿m3 / a,占总补给量的11. 2%。由此可以看出,对于浅层含水层来说,主要的补给资源来自降雨入渗补给。在排泄量中,蒸发量为50. 1亿m3 / a,占总排泄量的31. 9% ;河道排泄量为88. 5亿m3 / a,占总排泄量的56. 5% ;海水排泄15. 2亿m3 / a,占总排泄量的9. 7% ;边界侧渗出1. 8亿m3 / a,占总排泄量的1. 1%。由此可以看出,在地下水未开采状态下,河道排泄量是研究区的主要排泄量。另外,值得说明的是,由于没有收集到1965 年的降雨数据,所以采用了多年长时段降雨量的平均值,平均值可能要比1965年实际的降雨量大,因此计算出的地下水补给河流量也相应变大。表3表示的是三级分区水量平衡值,由于篇幅的限制,本文只列出了降雨补给量和蒸散发量以及与实测值的相对误差。

综上所述,从模拟期等值线拟合、地下水系统均衡分析来看,研究区地下水稳定流数值模拟模型基本上反映了海河流域平原开采初期状态下的地下水流动规律和特征,验证了边界条件和初识条件的准确性,符合研究区实际的水文地质条件,故可利用稳定流模型模拟的数据进行地下水开采量的模拟。

4 讨论与结论

海河流域目前地下水已经严重超采,形成了大片的降落漏斗,为帮助该区域地下水的管理,有必要对目前的地下水状态开展模拟工作,而本文通过对开采初期地下水的模拟建立的稳定态模型,可以较清楚地了解60年代平原区地下水的地质条件、流场形式、主要补给排泄量等,可以为下一步非稳定态模型的建立提供计算条件和初始水头条件。研究成果主要包括以下几方面:

(1)大区域地层模型的建立。通过对200多眼水文地质勘探钻孔的概化和分析,建立了平原区地层模型,该模型可以详细、定量说明平原区的地下含水层分布、埋深、储水条件,为地下水资源评价提供基础信息。该研究表明了在较大区域范围内和复杂的水文地质条件情况下地层模型建立的可行性以及必要性。

(2)开采初期地下水模型的建立。模型考虑了平原区地下水开采初期降雨入渗补给、山前侧渗、地下水与河流间水交换、地下水潜水蒸发、地下水排泄入海等条件,基本反映了当时地下水的补给、排泄条件。模型揭示,海河流域平原区开采初期状态下地下水水流的总体趋势是自北、西、三个方向指向渤海湾,吻合于地形坡度走向。模拟结果对下一步非稳定流的建立提供了边界和初始条件。

(3)地下水水量平衡分析。水平衡模拟分析结果显示,平原区总补给量156. 6亿m3 ,其中,山前侧渗17. 5亿m3 ,降雨补给136亿m3 ,河道入渗补给3. 1亿m3 ;总排泄量为156亿m3 ,其中,河道排出88. 5亿m3 ,蒸发50. 11亿m3 ,排泄入海15. 2亿m3 ,黄河内滩侧渗排出1. 8亿m3。模型展示了该区域自然状态的地下水补、排条件及水量,为进一步认识该区域地下水条件提供有力的数据支持,同时,也为水资源管理、水环境的修复提供科学依据和计算工具。

参考文献

[1] 费宇红,张光辉,曹寅白,等. 海河流域平原浅层地下水消耗与可持续利用[ J ]. 水文, 2001, 21 (6) : 11 - 13.

[2] 水利部海河水利委员会. 海河流域水资源及其开发利用情况调查评价[ R ]. 天津:水利部海河水利委员会, 2000, 102.

[3] 费宇红,曹寅白,张光辉,等. 海河流域平原地下水消耗与环境地质问题[ C ]. 九五全国地质科技重要成果学术交流会, 2001, 640 - 643.

[4] 王秀杰,练继建. 海河流域超采地下水引起的生态环境效应及其生态恢复对策[ C ]. 2003年第一届全国水力学与水利信息学学术大会, 2003, 119 - 121.

[5] 韩再生. 海河流域地下水资源保护[ J ]. 水文地质工程地质, 2001, 28 (1) : 28 - 29.

[6] 陈崇希. 防止模拟失真,提高仿真性"是数值模拟的核心[ J ]. 水文地质工程地质, 2003, 30 (2) : 1 - 4.

[7] 张奇. 数值模型在地下水管理中的应用[ J ]. 水文地质工程地质, 2003, 30 (6) : 72 - 78.

[8] 薛禹群, 吴吉春. 面临21 世纪的中国地下水模拟问题[ J ]. 水文地质工程地质, 1999, 26 (5) : 1 - 3.

[9] Groundwatermodeling system, Training manual[M ]. Beijing Normal University: 2002, Sep tember: 16 - 19.

[10] 水利部海河水利委员会. 海河流域水资源综合规划[ R ]. 天津:水利部海河水利委员会, 2002.

[11] 张光辉,费宇红,刘克岩,等. 海河平原地下水演变与对策[M ]. 北京:科学出版社, 2004.

[12] 国土资源大调查项目:华北地下水可持续利用开发前景专题报告[ R ]. 石家庄:中国地质科学院水文地质环境地质研究所, 2005.

[13] SURPAC Vision地质数据库操作手册[ R ]. 澳大利亚:澳大利亚国际软件公司, 2002.

[14] 河北省地质局[R ]. 石家庄:河北省水文地质局, 1981.

[15] 华北地区地下水资源评价75 - 57 - 01 - 03专题组. 华北地区及山西能源基地水资源研究专题报告[ R ]. 北京:原地矿部水文地质工程地质研究所, 1990, 36 - 38.

[16] 郭卫星,卢国平,等. MODFLOW模块化三维有限差分地下水流动模型. 美国地调局水资源调查技术丛书之六[ R ]. 南京:南京大学,

1999.

Simula tion of Groundwa ter Explo ita tion at Early Stage and Ana lysis onWa ter Budget in Ha ihe River Ba sin

ZHANG Xiao - ming1 , XUE Li - juan1 , ZHANG Qi2 , L I J ian - xin1

(1. Haihe River Commission, Ministry ofWater Resources, Tianjin　300170;

2. Nanjing Institute of Geography & Limnology, Chinese Academy of Sciences, Nanjing　210008, China)

Abstract Overexp loitation of groundwater in Haihe River basin has caused significant environmental p roblems. The decline of groundwater resource also affects the sustainable development of the local economics. The shortage of water resources in Haihe River basin has become a major attention of the government. As one of the management tools, a groundwater model was constructed for the Haihe River basin. The model was used to simulate the groundwater exp loitation at early stage in the area. The model indicated that the annual inflow to groundwater is 156. 6 ×108m3 , ofwhich 17. 5 ×108m3 for lateral flux from mountains, 136 ×108m3 for rainfall recharge, 3. 1 ×108m3 for river recharge. Groundwater discharges via rivers ( 88. 5 ×108m3 ) , evapotransp iration ( 50. 1 ×108m3 ) , sea (15. 2 ×108m3 ) and inner flooded area of Yellow river (1. 8 ×108m3 ). Groundwaterwas at a natural balance. Groundwater flows from the north, west and south boundaries to the Bohai Sea, which is consistent with the slop ing direction of the land surface. The modeling results p rovided data forwater resource management and remediation ofwater environment. The model also p rovided initial head condition, physical parameter values and boundary conditions for the simulation of the mass exp loited condition as the next step.

Key words: groundwatermodel; MODFLOW; p rimeval exp loited groundwater; Haihe River basin

作者简介:张晓明(1960 - ) ,男,高级工程师,主要从事水文地质及地下水模拟方面的研究。

来源：《干旱区资源与环境》2008.9