基岩裂隙水数值模拟探讨——以七台河市特大型应急水源地为例

祁福利，李永利，张峰龙，鲁守刚，王宏存

(1.中国地质大学(武汉)，湖北武汉 430047;2.黑龙江省九〇四水文地质工程地质勘察院，黑龙江哈尔滨 150027;3.黑龙江省地下水文学省级领军人才梯队，黑龙江哈尔滨 150027)

近年来由于降水量普遍偏小，作为七台河市三区一县城区唯一供水水源地的桃山水库蓄水水位、蓄水量持续下降。为确保城市居民生产生活用水问题，在桃山水库上游倭肯河河谷发现一处地下水类型为白垩系砂岩裂隙水的大型水源地。为确保该水源地的可持续开发利用和极端干旱年份的应急供水需求，准确可靠的确定该大型水源地的可开采量并对地下水水位的变化情况进行预测预报十分必要。因此，建立可靠的地下水流数值模拟模型对当地地下水资源及其环境问题进行综合评价，是十分必要和迫切的，具有重要的理论意义和实用价值。

目前，国内外专家学者开发了许多功能多样的地下水系统数值模拟软件，以其模块化、可视化、交互性、求解方法多样化等特点得到广泛的使用，其中以MODFLOW最为突出。不同模型各有优势，而MODFLOW的数据结构更易于实现与GIS的整合［1］;可从多角度校准模型，目前已经通过补给、排泄和水均衡的评估对模型进行了合理的校准，并达到广泛应用和推广［2］;应用范围广，从传统的模拟已经推广到区域地下水系统的水化学和同位素变化情况模拟，并结合实例进一步强调了做细致准确的稳态流分析对于瞬时流分析的重要性［3］。

近几年，我国学者在开发新模块上也做了大量的工作，主要针对以往模拟模型中暴露出的种种问题，在新理论和新方法的指导下，不断创新，开发新的模块，以提高模拟结果的可靠性和软件的实用性和快捷性。

本文以七台河市特大型应急水源地为例采用MODFLOW对区内的基岩裂隙水进行数值模拟探讨，进而为七台河市特大型应急水源地提供合理的地下水可开采量以及极端干旱年份时的地下水应急可开采量，并对不同地下水开采量的地下水位变化情况进行预测预报，进而实现地下水资源的合理利用与当地生态环境的可持续发展。

1 地质概况

七台河市属大陆性气候，多年平均气温为4.8℃，多年平均降水量为553.1 mm/a，主要降水月份集中在7—9月。

1.1 地质特征

第四系广泛分布于倭肯河河谷平原内，包括上更新统顾乡屯组，岩性上部黑灰色、黄灰色粉质黏土，下部含砾细—粗砂，砂砾石为主，全组地层厚度10～28 m;全新统温泉河组，岩性上部灰黄、黄色粉质黏土，含少量砂、砾，厚度1～4 m;下部为浅黄色含砾中细沙，含砾中粗砂，砂砾石互层，厚度小于10 m。现代河床冲积层，厚度小于15 m，岩性下部为青灰色、灰黄色含砾中粗砂、黑白杂色砂砾石互层为主，上部为灰黄色含砂粉质黏土。

在区内地表出露及河谷平原下伏被揭露的地层主要是白垩系下统，自老到新分为城子河组、猴石沟组。城子河组下段为海陆交互相含煤地层，中段为灰黄色中、细砂岩、粉砂岩与多层煤交互组成，上段为粉细、中砂岩和多层煤。底部在多数情况下都以或厚或薄的砾岩、含砾砂岩、中粗粒砂岩为主。猴石沟组岩性主要由黄、灰绿杂色凝灰质粉砂岩、中粗砂岩、中粗砾砂岩，泥质粉砂岩夹煤线等组成;中段夹流纹质凝灰岩，底部为底砾岩。全组明显分三段，构成一个从粗至细的沉积韵律，为河湖相沉积。

由于中生代以来的构造运动对沉积地层进行了多次的水平挤压与拉张作用，而且还有不同幅度、间歇性的垂向上的断陷与抬升作用，形成背斜轴部断裂构造带。根据物探结果，沿倭肯河河谷解译出10条断裂破碎带，为地下水的赋存提供空间，为地下水的运移提供有利通道(图1)。

1.2 水文地质特征

图1 断裂带分布示意图Fig．1 Distribution of fault zones

研究区地貌类型主要为泥沙质河谷平原、砂岩丘陵和玄武岩熔岩台地，赋存地下水类型为第四系松散岩类孔隙水、白垩系砂岩裂隙水和新近系玄武岩孔洞裂隙水三种类型，水文地质界线与地貌界线、岩性界线基本统一，值得注意的是在河谷平原地区存在双层含水层结构，上部为第四系松散岩类孔隙水含水岩组，下部为白垩系砂岩裂隙水含水岩组。

研究区地下水主要以大气降水入渗补给为主，直接补给裸露的第四系松散岩类孔隙水含水岩组、白垩系砂岩风化裂隙水含水岩组和新近系玄武岩孔洞含水岩组，埋藏型的白垩系砂岩构造裂隙水不直接接受大气降水入渗补给，依靠上部含水层和风化带沿垂向构造裂隙入渗补给。

进入到各含水岩组的地下水，总体上受地表分水岭控制，由分水岭向边界径流，相邻含水岩组间存在越层径流，其中在地表分水岭及其两侧地下水由上伏含水岩组向下覆含水岩组径流。

研究区内地下水的排泄方式有，潜水蒸发排泄，各含水岩组的人工开采，局部边界侧向径流排泄，被河谷切割处各含水岩组地下水向地表水体的排泄等。

(1)第四系松散岩类孔隙水

根据施工的钻孔抽水试验资料，可将第四系松散岩类孔隙水分为两个富水性分区:水量较丰富区，涌水量100～1 000 m3/d;水量较贫乏区，涌水量＜100 m3/d。

(2)新近系玄武岩孔洞裂隙水

主要分布在西北部熔岩台地区，船底山组玄武岩呈致密块状、隐晶质结构，具气孔状或杏仁状构造，风化裂隙和柱状节理较发育，含水不富水，透水性极好，在局部地段的低洼处地下水沿下伏基岩顶面以泉的形式流出地表，泉水流量较小，一般均小于10 m3/d。在有隔水层与阻水岩层的接触部分可赋存且水量较大。

(3)白垩系砂岩裂隙水

主要含水岩组为白垩系猴石沟组(K1～2h)、城子河组(K1c)砂岩、砂砾岩。因燕山期岩浆作用、压应力挤压成岩较好，孔隙不发育;断裂构造与裂隙发育形成破碎带，因而赋存基岩裂隙水，其中破碎带中的裂隙带总体呈层状分布，大概为2～3层含水带，累计厚度可达30 m，其间夹有泥岩或砂质泥岩而具有承压性;断裂带呈北东向展布，发育深度一般在170 m。

结合钻孔资料，倭肯河干流河谷平原地区造裂隙水非常发育，水量丰富区，单井涌水量1 000～3 000 m3/d(图2)。

图2 水文地质剖面图Fig．2 Section map of hydrogeology

2 七台河市特大水源地数值模拟

2.1 水文地质概念模型

本次模拟区将桃山水库也包含在内，模拟面积共计120 km2，模拟区地貌类型主要为第四系冲洪积河谷平原(Ⅰ区)、砂岩丘陵(Ⅱ区)和玄武岩熔岩台地(Ⅲ区)(图2)，各分区基本信息详见表1。

表1 渗透系数和给水度分区初值表Table 1 Initial-value table of zoning of permeability coefficient and specific yield

图3 模拟区地下水系统划分图Fig．3 Division of groundwater system in the model area

(1)含水层概化

第四系松散岩类孔隙水分布在倭肯河及其支流的河谷地带;新近系玄武岩风化裂隙水分布于模拟区的东北部;白垩系砂岩风化裂隙水分布于河谷两侧。根据物探资料，模拟区内倭肯河第四系下伏白垩系中张性构造断裂发育，形成了沿河谷呈北东向分布的构造裂隙富水带。这一赋水带从上至下呈现出3个裂隙发育富水带，第一层破碎带为风化裂隙及构造裂隙非常发育层位，埋深40～54.3 m，厚度4.7～10.2 m。第二层和第三层为构造裂隙发育层位，其中第二层埋深122～130 m，厚度3.7～11.6 m;第三层埋深147～162.3 m，厚度3.8～16.5 m。3个富水层位累计厚度约14.8～33.0 m。模拟区白垩系中除了河谷地带的张性构造断裂发育带之外的地区富水性一般。因此本次模拟将每个构造裂隙发育的富水带概化为渗透性良好的孔隙介质处理。对于和每个构造裂隙发育层位同一层的模拟区内的其他地区概化为透水性较差的孔隙介质处理。

(2)边界条件概化

在模拟区，各河流流入流出口为第四系地下水的流入流出边界，可以作为流量边界，全区边界上白垩系流入流出边界，同样概化为流量边界。

模拟区东北侧边界定义为给定流量边界，地下水均由外面流向模拟区;模拟区西南侧边界定义为给定水头边界，地下水均由该处流入水库。将东南、西北侧边界定义为给定流量边界，方向均为流进模拟区。倭肯河常年有水，与孔隙潜水水力联系密切，多数地段成为潜水的排泄边界，可概化为已知水头的一类边界。取潜水面为模型的上边界。地下水通过该边界接受河流、大气降水、田间灌溉回归水等入渗补给，并以蒸发的方式在该边界上产生排泄。将底边界确定为埋深180 m左右的位置，由于该深度处以下地层岩性为泥岩，裂隙发育不良，将该边界定义为零通量边界(图4)。

图4 模拟区水文地质概念模型图Fig．4 Conceptual model of hydrogeology in the model area

2.2 地下水流数学模型的建立

根据建立的水文地质概念模型，将计算区地下水流模拟的数值模型为非均质各向同性的孔隙潜水、裂隙承压水由越流联系的非稳定流混合模型，可以分别建立四层地下水的数学模型，用越流量将其耦合起来。由此所建立的模型可以看作是准三维的模型。

计算区面积120 km2，采用GMS进行100×100自动矩形剖分，剖分成7层，共剖分活动单元格22 246个，每个单元格200 m×200 m，面积0.04 km2。

模型识别期自2013年4月1日—9月30日，以一个月为一个应力期，每个应力期包括6个时间步长。根据统测资料确定的2013年4月1日含水层的初始流场(图5)。

图5 模拟区地下水等水位线图(2013年4月)Fig．5 Contour of groundwater table in the model area(April 2013)

2.3 模型的识别与检验

2.3.1 模型的识别

选择2013年4月—9月(共6个时段)进行模型的识别，以一个月为一个时间段，每个时间段包括6个时间步长。该时段为丰水期，源汇项较多，地下水均衡要素包括降水入渗补给、地表水灌溉渗漏补给、侧向流入补给、河流渗漏补给、潜水蒸发排泄等，此期间地下水流场变化快，且变化幅度较大，模型识别后的流场特征能较好地反映出含水层结构、水文地质参数和含水层边界性质的变化。对计算出的地下水水位与实测水位拟合误差进行统计，计算水位与实测水位在河谷地区的整体拟合程度良好。

2.3.2 模型的检验

本次采用多年平均状态下的地下水源汇项，输入模型运行20年，分析地下水流场的变化。分别选取上、中、下游段的代表性钻孔分析;而以运行20年后地下水流场作为模型验证的末期流场，与2013年同期水位进行对比，在20年的模型验证期内，地下水位呈周期性变化，总体保持不变，说明含水层结构、边界条件的概化、水文地质参数的选取是合理的，所建立的数学模型能较真实的刻画出模拟区地下水系统特征，可以利用该模型对水位的变化进行预测预报。

2.4 地下水位的预报

2.4.1 预报方案

本次对模拟区地下水位的预报根据水源地的使用方式，设计2个方案:

方案1:作为常规水源地开采，水源地可开采资源量5.2×104m3/d，模型运行20年;

方案2:作为应急供水水源地，(只在极枯年(P=97%)，降雨量为303.87 mm，桃山水库仅死库容1 900×104m3)确定某一应急取水量，在该取水量连续运行1年，地下水漏斗中心水位下降50 m。

2.4.2 预报结果分析

为了分析应急水源地在该开采条件下地下水流场的变化，将水源地所在河段分为上、中、下游三段分析。

2.4.2.1 方案1预报结果分析

此方案下，按照资源评价结果确定开采量5.2×104m3/d进行开采，将有关数据输入经过识别和验证的地下水流数值模型进行水位预报，预测到2014年、2018年、2023年、2028年和2033年地下水位的变化情况。模型结果显示在开采5.2×104m3/d开采基本没有出现含水层疏干现象，因此采用5.2×104m3/d可作为可持续开采量的参考依据。

为了分析新建水源地在该开采条件下地下水流场的变化，将水源地所在河段分为上、中、下游三段分析，上游段选择代表钻孔ZK100，中游段选择代表钻孔ZK29，下游段选择代表性钻孔ZK2。

(1)上游段钻孔ZK100:随着水源地地下水的开采，在前期下降速度较快，第一年水位下降0.92 m，到第5年水位下降1.49 m，第10年水位下降1.67 m，15年水位下降1.74 m，20年水位下降1.77 m。从水位下降量随时间的变化关系可以看出，水位在前期下降迅速，第10年以后基本保持不变，说明由于水源地的开采，地下水的降落漏斗基本稳定，不再继续下降，该处的水位在该开采量下的水位下降量维持在1.80 m(图6)。

(2)中游段钻孔ZK29:随着水源地地下水的开采，同样在前期下降较快，第一年水位下降1.84 m，到第5年水位下降3.26 m，第10年水位下降3.64 m，15年水位下降3.73 m，20年水位下降3.77 m。从水位下降量随时间的变化关系可以看出，水位在前期下降速度较快，到了第10年以后基本保持不变，说明由于水源地的开采，地下水的降落漏斗基本稳定，不再继续下降，该处的水位在该开采量下的水位下降量维持在3.80 m(图7)。

图6 预测ZK100地下水位历时曲线Fig．6 Forecasting of duration curve of groundwater table in ZK100

图7 预测ZK29地下水位历时曲线Fig．7 Forecasting of duration curve of groundwater table in ZK29

(3)下游段钻孔ZK2:随着水源地地下水的开采，地下水位总体变化量不大，主要原因是由于受到了下游桃山水库的影响。第一年水位下降1.00 m，到第5年水位下降1.35 m，第10年水位下降1.37 m，15年水位下降1.39 m，20年水位下降1.40 m。从水位下降量随时间的变化关系可以看出，水位下降量很小，到了第5年以后基本保持不变，说明由于水源地的开采，地下水的降落漏斗基本稳定，不再继续下降，该处的水位在该开采量下的水位下降量维持在1.50 m(图8)。

图8 预测ZK2地下水位历时曲线Fig．8 Forecasting of duration curve of groundwater table in ZK2

此方案下进行地下水开发利用，含水层没有出现疏干现象，河谷地区地下水位全部出现下降，但是下降深度不大。水源地地下水位下降量中游大，上游和下游逐渐减少。靠近桃山水库处地下水位基本不变(图9)。

2.4.2.2 方案2预报结果分析

根据方案2，若该水源地作为应急供水水源地，一般情况下不会启用，只有极端干旱年份——水库水位降至171.57 m时(死库容对应水位)，由于水库库底崎岖不平，水面难以连续，呈现多个分散的“水泡子”现象，但仍有1 900×104m3死库容难以抽取。这时启用应急供水水源地，按理论推算，如果地质条件可行的情况下，1 900×104m3死库容完全补给地下水，按照15×104m3/d的供水计算，死库容仍可维持一年，因此利用地下水数值模型根据方案2研究水库死库容水与裂隙地下水的水力联系极为重要。

为了分析新建水源地在该开采条件下地下水流场的变化，确定降落漏斗的中心位置，上游段选择代表钻孔ZK100，中游段选择代表钻孔ZK29，下游段选择代表性钻孔ZK2地下水位随时间的变化趋势。

(1)降落漏斗中心:通过调整开采井的开采量，最终确定当地下水的开采量达到18.2×104m3/d时，地下水的降落漏斗中心位置位于中游钻孔ZK50附近，水位下降量49.40 m(图10～图11)。

图9 地下水水位预测图(a)1年后、(b)5年后、(c)10年后、(d)20年后Fig．7 Forecasting of groundwater table(a)1 year later(b)5 years later(c)10 years later(d)20 years later

图10 预测ZK50地下水位历时曲线Fig．10 Forecasting of duration curve of groundwater table in ZK50

图11 2014年地下水水位预测图(12个月后)Fig．11 Forecasting of groundwater table in 2014(12 months later)

(2)上游段钻孔ZK100:随着水源地地下水的开采，在前8个月水位基本呈均匀下降，8个月以后下降趋势增大4倍，3个月后水位下降3.3 m，6个月后水位下降到6.45 m，9个月后水位下降18.45 m，1年后水位下降46.11 m。从图12可以看出，前8个月水位降深为0.82 m/月，而之后的5个月水位降深为9.15 m/月，表明前8个月水库补给充分，大量抽取地下水得到水库的及时补给，当连续抽取8个月，水库死库容逐渐消耗，水头差减小，补给能力衰减，逐渐开始抽取第一层破碎带中的净储量，导致水位急剧下降，模型继续运行至1年，水位急剧下降约50 m，达到第一层破碎带水的底板，第一层破碎带疏干。因此，采用18.2×104m3/d进行开采，可以保障应急1年的用水量(图12)。

(3)中游段钻孔ZK29:随着水源地地下水的开采，在前期6个月水位基本呈均匀下降，到6个月后下降趋势增大4倍，3个月后水位下降2.60 m，6个月后水位下降到15.27 m，9个月后水位下降至21.06 m，1年后水位下降40.67 m。分析6个月后地下水位下降深度增大原因同ZK100。从水位下降量随时间的变化关系可以看出，水位下降量基本与ZK100相同。因此，采用18.2×104m3/d进行开采，可以保障应急1年的用水量(图13)。

图12 预测ZK100地下水位历时曲线Fig．12 Forecasting of duration curve of groundwater table in ZK100

图13 预测ZK29地下水位历时曲线Fig．13 Forecasting of duration curve of groundwater table in ZK29

(4)下游段钻孔ZK2:随着水源地地下水的开采，3个月后水位下降1.39 m，6个月后水位下降1.76 m，9个月后水位下降至2.6 m，1年后水位下降15.49 m。从水位下降量随时间的变化关系可以看出，水位在前期较慢，8个月后水位下降明显加快，但水位下降量较中游和上游要小很多，主要原因是下游靠近桃山水库，在水位下降后更容易接受桃山水库的地表水补给。1年后桃山水库补给量平均为12.20×104m3/d(图14)。

图14 预测ZK2地下水位历时曲线Fig．14 Forecasting of duration curve of groundwater table in ZK2

此方案下进行地下水开发利用，上覆潜水含水层全部疏干，地下水位在后期下降明显，水源地地下水降落漏斗在中游地段，若应急开采1年后桃山水库补给地下水的平均量为12.20×104m3/d。若采用该方案进行开采，对当地地下水资源和桃山水库的库容具有很大的影响，因此建议只在应急时开采，应急开采时间不宜大于12个月。

该方案情况下，9个月后桃山水库补给地下水量平均为8.5×104m3/d，一年后桃山水库补给地下水量平均为12.2×104m3/d;前9个月累计水库补给量1 213.6×104m3，理论上前12个月累计水库补给量2 736.2×104m3，由于水库初始库容 1 930 ×104m3，在11月末水库死库容疏干，完全失去补给能力。因此，采用当地地下水资源和桃山水库的库容联合运营模式开采18.2×104m3/d维持一年是有保障的。

3 结论

本次研究在深入分析七台河倭肯河河谷桃山水库上游区地质及水文地质的基础上，利用MODFLOW软件建立了该研究区的地下水流模型，通过识别和验证后，提出两种不同方案对未来长期开采及应急开采条件下，地下水的变化趋势进行了预测。得出以下结论:

(1)建立了研究区多层富水带结构的地下水流数值模拟模型，对研究区白垩系砂岩构造裂隙水、第四系松散岩类孔隙潜水进行了数值模拟，利用实测资料对模型进行识别和验证，识别和验证效果良好，建立的水流模型达到了较高的精度，因此可以作为地下水流预报模型。

(2)研究区的地下水资源主要富集在白垩系砂岩构造破碎带中，构造裂隙富水带与下游桃山水库具有强烈的水力联系，在上游地下水集中开采的驱动下，水库能及时响应补给地下水。应急开采初期水位下降趋势相对较缓，而在大量集中开采的条件下后期水位下降迅速。反映出研究区开采地下水初期来自构造裂隙富水带赋存的水量，随着开采时间的增加，地下水开始接受下游水库的补给和上游激发条件下的侧向补给量。而且应急开采量越大，接受水库补给发生的时间越早，补给量越大。

(3)根据不同预测方案结果得出了长期可持续利用前提下的地下水开采量及不同应急开采量下的应急开采时间。

(4)从该地区地下水的长期可持续利用要求出发，该区日开采量以不超过5.2×104m3为宜。在该开采条件下，地下水位的下降深度较小，全区水位降深均维持在4.0 m以内，对研究区及周边环境不会造成不良影响。

(5)研究区作为应急水源地时应采用陶山水库与地下水联合运营模式，全力开采应急水源地，必要时可以动用桃山水库的死库容，开采量按18.2×104m3/d算保障开采一年。

研究区第四系松散岩类孔隙水和白垩系砂岩裂隙水具备较大的地下水开发利用潜力，地下水开发利用前景良好，适宜作为应急供水水源地。研究区地下水资源为正均衡，但属于气象型水源地，水源地补给来源主要为大气降水，虽可以短时间内大规模开采但是不适宜长时间大规模开采，以丰补歉、多年调节，方可作为应急水源地用以保障完达山七台河地区供水安全。

［1］ JUANC S， KOLMK E. Conceptualization，characterization and numerical modeling of the Jackson Hole alluvial aquifer using ARC/INFO and MODFLOW［J］.Engineering Geology，1996，42:119-137.

［2］ OLSTHOORNTN.A comparative review of analytic and finite difference models used at the Amsterdam Water Supply［J］.Journal of Hydrology，1999，226:139-143.

［3］ HARRINGTONG A，WALKERG R.A compartmental mixing-cell approach for the quantitative assessment of groundwater dynamics in the Otway Basin［J］.Journal of Hydrology，1999，214:49-63.

［4］ 卢文喜.地下水运动数值模拟过程中边界条件问题探讨［J］.水利学报，2003(3):33-36.［ LU W X.Discussion of the boundary conditions in numerical simulation of groundwater movement［J］.Journal of Hydraulic Engineering，2003(3):33-36.(in Chinese)］

［5］ 寇文杰.格尔木河流域地下水数值模拟［J］.水文地质程地质.2013.40(1):34-40.［KOU W J.Numerical simulation of the groundwater system in the Geermu river basin［J］.Hydrogeology ＆ Engineering Geology，2013，40(1):34-40.(in Chinese)］

［6］ 杨青春，卢文喜，马洪云.Visual Modflow在吉林省西部地下水数值模拟中的应用［J］.水文地质工程地质，2005，32(3):67-69［YANG Q C，LU W X，MA H Y.Application of Visual Modflow in Western Jilin Province numerical simulation of groundwater［J］.Hydrogeology ＆ Engineering Geology，2005，32(3):67-69.(in Chinese)］

［7］ 魏亚强，乔小娟，李国敏.MODFLOW不同算法及参数设定对计算精度的影响［J］.水文地质工程地质，2015，42(1):14-21. ［WEI Y Q，QIAO X J，LI G M.Influence on calculation accuracy caused by different algorithms and parameters in MODFLOW［J］.Hydrogeology ＆ Engineering Geology，2015，42(1):14-21.(in Chinese)］