基于数值模拟及可开采系数法的广花盆地地下水可开采量研究

李 丹， 李志威， 张明珠， 庞 园

(广州市水务科学研究所， 广东 广州 510220)

1 研究背景

广花盆地是广州市重要的地下水应急备用水源地[1-2]，同时也是珠江三角洲岩溶塌陷及地面沉降地质灾害重点防治区[3]。盆地内自来水暂未全面供给，分散性开采岩溶水的工业及商业机井仍较多，由于不注意合理布局和科学开采，历史上广花盆地局部地区诱发了地面沉降、塌陷、地裂缝、房裂等地质问题[4-6]。如何让这一重要的应急备用水源地在水量和水质上满足人类社会和生态系统长期稳定发展的需要，且能有效避免地质-生态灾害，则需要对地下水资源可持续开发利用问题进行探讨。

国内学者对广花盆地开展过地面塌陷成因、特征、形成机理、风险评估、监测设计与防治方面[3-8]以及水质评价及趋势[1-2]的研究，然而尚未对广花盆地地下水可开采量进行计算及预测，不利于实现地下水资源可持续开发利用。

地下水系统的状态随着水文气象、开采强度、开采地段和开采层位等因素的变化而变化，因此确定地下水合理开采量，关键是要选择能够真实客观地反映近期地下水资源状况的开采量。国内诸多学者使用数值模拟的方法计算不同研究区的地下水可开采量，并通过模拟研究区按照设计开采量进行开采时地下水水位的变动规律、降落漏斗的发展，论证设计开采量的可行性[9-13]。本文使用数值模拟及开采系数法计算广花盆地的可开采量，两种方法相互验证，论证设计开采量的合理性。本文可为实施地表水-地下水联合开发提供基础技术支持，对实现地下水高效利用与可持续发展、保障广州市供水安全、改善环境地质问题具有重要意义。

2 研究区概况

广花盆地的地理范围为北纬23°00′～23°30′，东经 112°57′～113°26′，北高南低，面积约858 km2，以平原为主，台地和丘陵次之。地下水可分为松散岩类孔隙水、岩溶水、基岩裂隙水3大类。

大气降雨是研究区地下水补给的主要来源。平原区除主要接受大气降雨的入渗补给外，还接受地表水入渗、灌溉入渗、潮水顶托反补给、基岩山区及丘陵台地地下水的侧向补给。基岩山区及丘陵台地的地下水主要接受大气降雨的入渗补给，山塘、水库及渠道等对其入渗补给的规模较小。

总体上研究区地下水从北向南流动。外围丘陵山区一般地势较高，浅部基岩裂隙水获得补给后经过短暂的径流，便以泉或渗流形式排入附近的溪流；深层地下水则通过断层和裂隙向广花盆地汇流。

丘陵山区基岩裂隙水多以泉或渗流形式向沟谷排泄。江村-龙归一线以北由于地下水位高于河水位，孔隙水和岩溶水主要向河流排泄。江村-龙归一线以南地区地势平坦，河涌较发育，丰水期由于地下水位埋深浅，地下水受潮水顶托影响，使地表水与地下水之间呈一定互补关系；枯水期地下水位埋深较深，地下水接受河涌水的入渗补给。此外，地下水的排泄方式还有开采和地表蒸发等。

21世纪60至80年代期间，广花盆地城镇供水地下水所占比例很高，农村用水及灌溉用水大量使用地下水，采石场抽取地下水特别突出，建立了江村、肖岗、石井、新市、新华、雅瑶等多个水源地。

3 概念模型

3.1 模型边界

第四系含水层北部和东部边界为平原和山区分界线，山前第四系的孔隙介质富水性好，容易接受山区降水潜流形式的侧向补给，模型中将该边界设为给定流量边界；研究区西部边界是芦苞涌和白坭河，河流与潜水有较好的水力联系，模型中将该边界设为给定水头边界；地势北高南低，通过地下水流场判断南部为流量流出边界，其他实测资料缺乏的边界概化为通用水头边界。

岩溶含水层是模拟区主要的集中开采层，补给主要通过第一含水层越流补给，小块岩溶出露地区接受降水补给，南部边界为地下水流出边界，不同时期地下水等水位线表明，岩溶层在不同年份和季节，存在着地下水的流出，由于广花盆地三面被断裂包围，而且两侧多为基岩，或两侧大部分地层不衔接，设定为隔水边界。模拟范围及地下水水位观测井分布情况见图1。

潜水的自由水面为系统的顶部边界，地下水主要以垂直运动为主，第四系和浅层岩溶地下水水力联系比较密切。钻孔揭露岩溶层底板以下岩溶水流动缓慢，几乎没有地下水的开采，因此将钻孔打到的位置作为本次数值模型的底部边界。

3.2 模型结构

依据研究区内水文地质勘察成果(主要为开采井及观测井的钻探成果)将整个含水系统自上而下分为3层，第1层是潜水含水层，主要的岩性为素填土及细砂等，其底板标高在-30～35 m之间；第2层为相对隔水层，主要岩性为粉质黏土，其底板标高在-50～25 m之间；第3层为岩溶含水层，主要岩性为碳酸盐岩，其底板标高在-100～- 50 m之间。

图1 模拟范围及观测井分布图

3.3 水文地质参数

潜水层水文地质参数按照岩性的不同共划分为8个参数区，在河流及其支流的两侧以及冲洪积扇中上部由于岩性主要是粗砂和砾石，故渗透系数K取值范围为4～50 m/d，河道中下游逐渐过渡到细砂，中下部地区K取值较小，给水度取值范围为0.02～0.2。第四系弱透水层水文地质参数按照岩性的不同共划分为6个参数分区，渗透系数取值范围为0.014～3 m/d，弹性释水系数取值范围为1×10-5～3×10-4。岩溶含水层水文地质参数根据抽水试验的结果，用GIS建立泰森多边形，其渗透系数取值范围为0.03～123 m/d，弹性释水系数取值范围为1×10-5～3×10-4。

3.4 源汇项

研究区的补给项主要包括降水入渗补给、河流与水库入渗补给和灌溉入渗补给。降水入渗补给是主要的补给来源，河流补给以及灌溉入渗补给较小。排泄项主要有泉、地下水开采、蒸发、侧向排泄，其中地下水开采是主要的排泄方式，根据研究区内的工厂和农田布局情况，结合计算需要，做如下概化：工业区主要集中在盆地南部区域，主要开采岩溶承压水，因此均匀布设部分承压水开采井和一定数量的潜水开采量。

4 数学模型

4.1 数学模型的建立

岩溶发育的不连续性、非均匀性和各向异性，加之实际勘探资料的制约，使得岩溶地下水建模成为岩溶水研究的难点和热点。广花盆地的含水岩组主要是可溶性的碳酸盐岩，水流运动理论上是不符合达西定律的，但是盆地的岩溶发育相对比较均匀，因此仍应用等效多孔介质数值模拟模型分析广花盆地的地下水流运动。该方法基于达西流假设和典型单元体的概念，用等效水力参数表征含水介质渗透特性。大量的实例表明，以岩溶地下水系统水资源调查评价为目的，利用等效介质方法模拟岩溶裂隙含水介质是有效可行的。MODFLOW、FEFLOW、GMS等基于连续孔隙介质的模拟软件，已有很多成功应用案例。

从空间上看，地下水流整体上以浅层垂向运动与深层水平运动相交替为特征流动，为了用模型分析不同开采条件下含水层之间的互动，在水文地质概念模型基础上，建立潜水和承压水水流数值模型。本研究采用地下水模拟系统GMS。对研究区进行水平剖分，研究区被剖分为160行、160列，共计有4 462个有效单元。模拟期为2011年4月至2013年3月。根据观测井的水位资料，确定本次数值模拟的初始流场。由于第四系和下伏岩溶层有密切的水力联系，监测资料也表明潜水水位和岩溶层的水位在前汛期的4月基本一致，故主要观察岩溶地下水初始流场(图2)。

图2 岩溶含水层初始流场图

4.2 模型的校正与验证

采用等效多孔介质单一连续方法的关键是确定等效的水力参数。本文使用手工校准(正向校准)进行水流模型的校准，即手工调节各参数，直到模拟结果(地下水流场、地下水水位的变化趋势、边界条件、水文地质参数和含水层结构)与实地观测值基本一致[14-16]。反复调参后得到比较理想的模型校正结果。模拟时间总共730 d，分别计算第60、180、540和730 d的标准化残差均方根，分别为4.996%、4.043%、6.517%、6.787%。将第60、180、540和730 d的不同观测井的地下水水位计算值与观测值进行对比(图3)，计算值与观测值总体分布在y=x线两侧，表明观测井的观测值与计算值相差较小，地下水位的拟合程度较高。图4为地下水水位等值线拟合图，从图4可以看出，所建立模型总体反映了模拟区的实际情况，模型经校准得到的地下水流场与实际流场基本上一致，地下水水位模拟值与实测值的动态变化过程基本吻合。研究区地下水流场分布与地形趋势一致，自东北向西南，地下水位呈递减趋势。3口代表性观测孔(A1-12、A1-47和A1-23)的水位拟合过程线分别见图5～7，观测值总体分布在计算值附近，拟合效果较好。因此，所建立模型总体反映了研究区的实际情况，可用于地下水资源计算与预测。

图3地下水水位计算值与观测值对比图4地下水水位等值线拟合图(单位：m)

图5 A1-12观测井拟合曲线图6 A1- 47观测井拟合曲线图7 A1-23观测井拟合曲线

5 可开采量计算结果

5.1 数值模型法

盆地内岩溶含水层大部分为埋藏型，上覆厚度约5～25 m的第四系松散沉积物，属于岩溶塌陷易发区。地下水的运动会使覆盖层土体的力学性质发生改变，一般认为抽取地下水引起塌陷的水位约束值与抽水井所在地区的土层厚度相对，即水位降深以不超过岩溶含水层顶板为宜。广花盆地地下水水位至少高于第一含水层底板5 m，因此取降深5 m作为岩溶含水层抽取水量的约束条件。《广东省地下水保护与利用规划》表明，在岩溶覆盖区由于大幅度抽取地下水易引起地面塌陷，而控制开采水位降深则是避免大面积、大幅度塌陷的主要措施，开采水位降深宜控制在5～8 m 以内。因此本文确定岩溶开采水位降深最大值为5 m。

可开采量是在不大于该开采强度条件下，预测时段内的地下水平均水位能保持在安全开采水位降深之内。以此为基础，以降深为5m的水位约束条件对模拟区未来5年(2015年9月至2020年9月)地下水可开采量进行模拟预测。

在规划的开采条件下模拟末期研究区含水层流场结果(图8)显示，2020年广花盆地中心城区出现明显的降落漏斗区域，降落漏斗中心的水位均控制在给定的降深内，并且漏斗的中心是广花盆地开采井最密集的区域，由此可见供水水源地集中开采地下水是引起地下水水位降落漏斗空间位置变化的主要影响因素。图9为模拟期3个观测井水位过程线，由图9可知，中心城区附近的观测井(A1-22和A1-25)的水位随着时间逐渐降低，但是降低的幅度在减小，最大降深约为4 m，偏远地区的监测井(A1-40)的水位变化幅度不大。

图8 研究区2020年9月模拟流场图

中心城区降深最大的开采井的水位预测如图10所示。随着开采井运行，水位开始下降，在3月份水位达到最小值，这时控制水位降深最小，为0.3 m，从4月开始地下水水位开始回升，控制水位降深也逐渐增大，直到9月份达到峰值。10月至次年3月，地下水水位又持续下降，之后又上升，如此循环，2016年9月水位降深达到最大值2.36 m。通过不断的调节抽水井的抽水强度，使预测时段内的地下水平均降深能保持在安全开采水位降深之内，模型计算得到的地下水允许开采量为12 997.78×104m3。

图9 模拟期3个观测井水位过程线

图10 中心城区观测井水位降深变化过程线

5.2 可开采系数法

广花盆地地下水资源量采用补给量法计算。

降水入渗补给量的计算公式为：

Q降=F·α·P

(1)

式中：F为可入渗面面积，m2，将研究区离散为若干1 km2的单元格，参与计算的各单元格面积是指该单元格扣除遥感解译得到的地表水系和硬化地面部分后的剩余面积；P为降水量，m，将降水量插值到研究区每一个单元格内；α为降水入渗补给系数，将每一个自动监测井的年内数据，绘制成表，选取次降水强度较大的若干场次降水，划定对应次降水时段的地下水水位增幅，根据公式(2)求取土壤的降水入渗补给系数。

α=μ· (∑ΔH)/P次

(2)

式中：μ为水位变动带含水层给水度；ΔH为次降水造成的水位增幅，m；P次为引起该水位增幅对应时段的次降水量，m。

灌溉入渗补给量的计算公式为：

Q灌=15a0·m·fg

(3)

式中：a0为灌溉入渗补给系数，农作物a0值为0.157，园林作物a0值为0.12；m为广州市灌溉定额量，m3/hm2；fg为灌溉面积，hm2。

平原区的山前侧向补给量计算公式为：

Q侧补=K·I·L·B

(4)

式中：Q侧补为山前侧向流出量，m3/d；K为剖面位置的渗透系数，m/d；I为垂直于剖面的水力坡度；L为山丘平原界线的长度，m；B为剖面含水层厚度，m。

补给量采用达西公式计算：

Q河补=K·I·L·B

(5)

式中：Q河补为单侧河道侧向补给量，m3/d；K为剖面位置的渗透系数，m/d，通过钻孔抽水试验成果插值确定；L为河段长度，m；B为计算剖面厚度，m，通过钻孔抽水试验成果插值确定；I为垂直于剖面的水力坡度。

经计算，盆地内白坭河、新街河、天马河的侧向补给量较大，分别为1 075.31×104、579.13×104、279.43×104m3。此外，盆地内流溪河段、珠江(西航道)段侧向补给量分别为136.82×104、6.76×104m3。区内河道侧向补给量共计2 077.45×104m3。

利用补给系数法计算湖、库、渠的入渗补给不仅需要确定入渗面积，还需判别其水位与地下水水位之间的大小关系。计算方法为：

Q湖库渠=F·I·ψ，I>0

(6)

式中：Q湖库渠为湖、库、渠及微小支流对地下水的补给，m3/a；F为计算单元格内纳入计算的湖、库、渠面积，m2；I为计算单元格与周边含水层之间的平均水力坡度；ψ为单位水力坡度的渗漏补给系数，m/a，本次计算中将区内水系划分为主要河道和湖、库、渠及微小支流，并分别计算其与地下水之间的补排量，事实上这种划分并没有严格的界定方法，人为主观因素明显。因此，对于河道侧向补给量，依据公式(5)与公式(6)的计算结果应较为一致。通过试算，当ψ=1.1 m/a时，两者结果基本一致。由此确定单位水力坡度的渗漏补给系数为1.1 m/a。

经计算，广花盆地各补给项结果见表1。表1中地下水资源量的计算公式为：

Q平=Q降+Q灌+Q河+Q湖库渠+Q山前

(7)

式中：Q平为平原区的地下水资源量，104m3；Q降为降水入渗补给量，104m3；Q灌为灌溉入渗补给量，104m3；Q河为河道侧向补给量，104m3；Q湖库渠为湖、库、渠渗漏补给量，104m3；Q山前为山前侧向补给量，104m3。则广花盆地地下水资源量为18 972.99×104m3。

表1 广花盆地地下水资源量计算结果 104 m3

本文根据1∶100000比例尺区域水文地质图资料确定广花盆地各计算单元的富水性，按各级富水性相对应的可开采系数范围，结合开采条件、抽水试验成果及实际开采量，确定可开采系数如表2所示。经计算，广花盆地地下水可开采量为14015.09×104m3。

表2 广花盆地可开采系数表

6 结论与建议

6.1 结 论

(1)本文运用水文地质知识分析数据资料，按照规范化的建模流程(GMS)构建了广州市广花盆地的地下水水流模型。结果表明，该模型可为广花盆地提供合理的地下水水流模拟结果，能揭示地下水流动规律，可用于地下水资源计算与预测。

(2)采用补给量法计算广花盆地的地下水资源量为18 972.99×104m3，结合广州市地下水基础资料条件，考虑开采条件、实际开采现状和地下水资源状况等因素，采用可开采系数法计算广花盆地地下水可开采量为14 015.10×104m3；以可开采水位最大降深5 m作为岩溶含水层抽取水量的约束条件，采用数值模型计算的地下水可开采量为12 997.78×104m3，两种方法计算的结果相对误差为7.3%，两种方法相互验证，表明本文可开采量计算结果较为可靠。

6.2 建 议

(1)结合已有的地下水工作成果及广花盆地水文地质条件，进行地下水超采区划定及评价。

(2)提出广花盆地应急水源地的开采控制目标及措施：

广花盆地岩溶地下水区一般情况下严禁开采、严格保护，但在启动应急预案或专供水明显不足时，并经论证不会产生环境地质问题的前提条件下，可适量开采。该目标与《广东省地下水功能区划》关于应急水源地管理目标相一致。

除应急备用水、矿泉水、地热水等对水温、水质有特殊要求的地下水外，不再批准抽取地下水取水许可申请。水行政主管部门对批准的地热水、矿泉水取水工程应核定地下水开采量和年度用水计划；对已批准的地下水为自备水源的取水户，要实行计划用水，采取相关节水措施，逐年核减地下水开采量。

对于在城市供水管网覆盖区域内，新建、改建或扩建工程项目自备取水设施取用地下水的申请不予批准。市区各自备水源单位原持有取水许可证到期后，城市公共供水管网能满足用水需求的，不再办理续证手续，不得再开采地下水。在公共供水管网未覆盖区域内，建设单位应经建设项目水资源论证、取水许可获得批准后，方可兴建地下水取水工程或者设施。

(3)加强地下水动态观测，逐步构建地下水水位、水温和水质的自动监测、信息共享及预警预报体系。