黑龙洞泉域岩溶地下水水质运移模拟分析

马志敬， 边 凯， 庞 宇， 刘 博

(1.河北工程大学水利水电学院, 邯郸 056038； 2. 中国煤炭地质总局水文地质局， 邯郸 056004; 3.河北工程大学地球科学与工程学院， 邯郸 056038)

黑龙洞泉域是邯郸市的主要饮用水源地之一[1-2]，由于长期的不合理开发利用，造成该泉域的地下水水位下降，形成了地下水位降落漏斗，水质变差，出现了许多地质问题[3]。地下水污染不同于地表水污染，当发现地下水被污染时其污染程度已经很高，很难进行恢复,由于人类活动对区域水资源可持续利用的影响会随着人类社会的发展和进步继续扩大[4]，所以合理开发利用和有效保护地下水资源就显得尤为重要。

建立地下水水源地保护区是保护地下水水质安全的有力手段[5]，虽然国内对地下水水源地保护区的研究主要借鉴国外的3级划分法[6]。但是随着研究的不断深入，中国水源地保护区的划分方法主要形成了经验值法、公式计算法、解析模型法、数值模拟法等[7]。近年来，岩溶水越来越受到学者的关注，常用的岩溶水数值模拟的方法主要有等效多孔介质法、双重连续介质法和三重介质法[8-9]。地下水数值模拟方法正逐渐取代传统技术,已成为现代研究地下水系统各种问题的重要技术手段之一[10-11],地下水数值模拟方法的各种模型及研究范围已经基本满足了国民经济发展建设的需要[12]。

利用地下水流数值模型分析评价地下水资源量的研究有许多，但针对地下水系统动态演变规律的区域研究就相对匮乏。白晓等[13]采用Visual MODFLOW建立峰峰矿区岩溶地下水流数值模型，在目前开采条件下，对研究区地下水资源量及均衡状态进行分析，并结合ARIMA降水量预测模型，预测了2018—2025年地下水位动态变化，为峰峰矿区地下水资源合理开发利用提供理论依据和参考；秦欢欢[14]通过建立北京市平原区三维地下水流稳态及非稳态数值模型，采用情景分析法模拟了北京平原2015—2030年的地下水流动及水均衡状态，对北京平原地下水储量、水均衡及可持续利用状况进行一定分析，制订出最优的地下水开采情景，以保证北京地下水资源的合理、可持续开发利用，控制地下水环境的继续恶化，为北京地下水的可持续利用和发展提供科学的指导。

史启朋等[15]和杨钊等[16]都利用Visual MODFLOW软件建立了污染物运移的数值模型，以模拟和预测研究区污染物运移过程与趋势特征，为地下水环境评价和治理提供了依据；张文喆等[17]将FEFLOW应用于徐州市奎河地下水流数值模拟，分析了河水与地下水的补排关系，为研究地下水污染提供了依据。

肖杰等[18]通过应用数值模拟的方法对崇州市城区饮用水源地进行了各级水源地保护区的划分；吴乐等[19]通过对北京西山地区进行区域含水层系统数值模拟，分析了在不同的开采条件下含水层系统响应特征，为地下水资源持续利用提供依据；孟茜等[20]通过采用地下水数值模拟的方法，对准格尔地区地下岩溶水资源量变化进行了预测评价，为本区的地下水合理开采提供了依据，对类似条件的地下水数值模型建立提供了一定的借鉴。研究表明，地下水数值模拟法在实际应用中结果准确可靠，具有广泛的实用性。但在地下水数值模拟中结合遥感影像来确定水文地质参数分区的研究相对匮乏。

现以黑龙洞泉域岩溶地下水为例，在充分了解地下水数值模拟法的理论和解决地下水问题的步骤后[21]，对研究区建立概念模型和数学模型，在已有资料的基础上结合遥感影像来确定含水层参数，对其数学模型进行识别和验证,通过模拟获取地下水流场和以开采井为中心溶质质点向外迁移的运动轨迹，按照《饮用水水源保护区划分技术规范》(HJ 338—2018)[22]要求划分各级水源地保护区，进而为有关部门提供合理的划分结果。

1 区域概况

黑龙洞泉域在河北省邯郸市峰峰矿区，西边紧邻太行山，华北平原在泉域的东面，总体地势趋于西北、西部高，东南、东部低。

泉域内从老到新发育的地层依次为：震旦系、寒武系、奥陶系、石炭系、二叠系、三叠系、古近系、新近系和第四系。其中震旦系和三叠系地层仅在局部出漏分布。泉域内以新华夏系为主要构造体系，构造方向主要是北东、北西及东西向，构造形态主要有褶皱和断裂。

2 地下水流数值模型构建

2.1 水文地质概念模型

泉域内主要含水岩系有寒武系、奥陶系碳酸盐岩裂隙岩溶含水岩系，石炭系、二叠系、三叠系碎屑岩夹碳酸盐岩裂隙含水岩系，第三系、第四系松散岩类孔隙含水岩系三种，其中寒武-奥陶系灰岩作为研究区的主要含水岩组。含水层中的岩溶裂隙非常发育，且都为非均质各向异性，由西北向东南逐渐变薄，厚度在109～230 m。上覆由石炭系本溪组的铝土质泥岩构成了相对稳定的隔水顶板，使地下水处于承压状态，但由于局部地区的人工开采，造成水位下降，地下水已由承压状态转化为无压状态。下伏的寒武系灰岩由于埋藏位置较深，节理裂隙极不发育且多被充填，因此形成了相对稳定的隔水底板。虽然它们之间都有相对隔水层,但是还会通过各种构造来发生水力联系,因此在区域上可以当成是统一的含水岩体。在研究区内除了有集中式开采的羊角铺水源地外，还有工业用水井、矿井疏干排水井、分散式开采的乡镇工农业及生活用水水源井。

根据前述并结合羊角铺水源地周边水文地质条件进行合理概化，东部边界：以岳城、新坡、中史村一线为边界，灰岩的埋深大多在标高-800 m以下，岩溶已经极不发育，地下水循环变得滞缓，所以设为阻水边界；西部边界：以由许多相对阻水及导水断层组成的贾壁断裂带为边界，该边界的西面有大面积的碳酸盐岩裸露区，在雨季时接受大气降水补给后，以侧向径流的形式通过导水断层流入该系统，因此设置为补给边界;南部边界：以老爷山背斜及陡贡-杨花庄断裂隆起带为边界，是一个地下水分水岭，阻水性能良好，将此边界设为隔水边界；北部边界：以八特至康二城西一线为边界，其中西段八特至杨二庄一线，由于苑城地堑使奥陶系灰岩埋藏加深，岩溶极不发育，为阻水边界，中段由于鼓山断层将断层两侧含水层错开，同样为阻水边界，东段奥陶系埋深较深，且北部有紫山阻水岩体存在，也可视为阻水边界，总体来说北部边界为阻水边界。

将研究区的地下水流概化为三维非均质各向异性的承压-无压非稳定流，地下水溶质通过对流-弥散随地下水流一起运移。水文地质概念模型如图1所示。

图1 水文地质概念模型图Fig.1 Hydrogeological conceptual model diagram

2.2 数学模型

选取北东18°主构造线方向-主渗透方向为x轴，根据上述水文地质概念模型，建立下列与之相适应的数学模型。

2.2.1 地下水流数学模型

根据已经建立的水文地质概念模型，用以下的数学模型进行描述：

(1)

H(x,y,z)t=0=H0(x,y,z,t0), (x,y,z)∈Ω

(2)

(3)

(4)

式中：kxx、kyy、kzz为含水层各向异性主方向渗透系数，m/d;H为点(x,y,z)在t时刻的水头值，m;H0(x,y,z,t0)为点(x,y,z)处的初始水位，m;Ss为贮水率，m-1;W为源汇项，d-1;t为时间，d；Ω为计算区；qw为自由面单位面积上大气降水入渗补给量与地下水蒸发量之和，m/d;cos(n,x)、cos(n,y)、cos(n,z)分别为流量边界外法线方向与坐标轴方向夹角的余弦；μ为饱和差(当自由面上升时)或给水度(当自由面下降时)；q(x,y,z,t)为第二类边界上单位面积的补给量，m/d。

2.2.2 溶质运移数学模型

(5)

C(x,y,z,t)t=0=C0(x,y,z,t0)， (x,y,z)∈Ω,t≥0

(6)

(7)

式中：R为迟滞系数；θ为介质孔隙度；C为组分的浓度，ML-3;T为时间，d;x,y,z为空间位置坐标，L；Dij为水动力弥散系数张量，L2/T;Vi为地下水渗流速度张量，L/T;qs为源和汇，T-1；Cs为源或汇水流中组分的浓度，ML-3;C0(x,y,z,t0)为已知浓度分布；Ω为模型模拟区域；Γ2为通量边界；fi(x,y,z,t)为边界上Γ2已知的弥散通量函数。

2.3 数值模型

2.3.1 研究区空间网格剖分

研究区空间离散采用等距或不等距的正交长方体对研究区进行网格剖分，从而得到了研究区的空间离散单元，对于需要提高模拟精度的区域进行局部加密剖分，最终在研究区平面上剖分成110×70的矩形网格单元，共计7 700个网格单元，其中有效单元为4 962个，无效单元格2 738个，平面剖分网格和垂向剖分网格如图2～图4所示。

图2 研究区平面剖分图Fig.2 Plane section map of the study area

图3 第32列剖面图Fig.3 Column 32 profile

图4 第51行剖面图Fig.4 Line 51 profile

2.3.2 源汇项处理

根据研究区地下水系统的补给、径流及排泄条件，最终确定大气降水入渗补给和侧向径流补给为研究区的地下水补给来源。大气降水补给主要分布在中部鼓山和西部广大的碳酸盐岩裸露区，入渗系数取0.5(根据河北省邯郸市水资源研究)。侧向径流补给主要在研究区西部边界通过导水断层流入该系统。

地下水排泄以人工开采井和黑龙洞泉群自然排泄为主要途径，排泄项采用抽水井的方式表示。根据自来水公司及峰峰矿区水利局、邯郸市水利局收集2018年地下水排泄量数据，如表1所示。

表1 地下水排泄量统计表

2.3.3 水文地质参数分区

已建立的数值模型把渗透系数、储水系数、给水度等作为主要的水文地质参数。早在20世纪70年代，专家学者就已对黑龙洞泉城研究做了大量的工作，如1976—2000年，煤炭部水文勘探指挥部、河北省电力勘测设计研究院等对峰峰矿区奥灰地下水疏放降压进行了水文地质勘探，并在羊角铺和王凤等地进行了多次抽水试验。通过以前多次的水文地质勘探和抽水试验，对黑龙洞泉域内的水文地质条件有了充分了解。故奥陶系灰岩含水层参数主要依据20世纪70年代以来多次水文地质勘探成果资料进行分区并赋值。根据黑龙洞泉域内水文地质条件差异，对该区域进行参数分区划分，如图5所示，共分为了14个参数分区，并赋予每个参数分区不同的给水度初值、储水系数Ss和渗透系数k，来作为模型模拟计算的基础。

黑龙洞泉域岩溶水系统补给区位于西部山区和中部鼓山灰岩裸露区，面积约1 262 km2。为减少因模型参数的不准确性带来的误差，在上述参数分区的基础上结合遥感影像(图6)，依据灰岩裸露区覆盖程度对其再进行分区，赋予不同的入渗系数，如表2所示，最终模拟计算时采用最准确的参数分区。

2.3.4 模型的识别与验证

模型识别与验证也称为“参数反演”，通过野外现场水文地质试验获得的水文地质参数是模型的解析解，与实际值相差较大，但水位观测数据是比较可靠、准确的，因此水位观测数据可用于反求水文地质参数。当计算水位与观测水位趋势相同时，说明修正后的水文地质参数可反映实际的水文地质条件。

图5 水文地质参数的分区图Fig.5 Zonal map of hydrogeological parameters

Ⅰ区为灰岩表面覆盖程度较大区域； Ⅱ区为灰岩表面覆盖程度较小区域图6 系统边界的遥感图Fig.6 Remote sensing diagram of system boundary

模型选择2017年1月1日—12月31日为模型的校正、识别期，分为10个制度期，每个制度期又分为10个计算时间步长。在全区共建立了20口水位观测井，其分布如图7所示。该模型选取了9眼地下水位观测井，因为9眼井大致分布均匀，具有长期的水位观测资料。

通过该模型的运行，根据预测水头值和观测井实测值的拟合情况，对参数不断进行试算检验，直到两者拟合较好为止，得到各分区的参数值。经校正、识别，各分区的参数识别结果如表3和表4所示。各地下水位观测井的水位拟合精度如图8所示，平均水位误差在1 m左右。

从上述拟合图中可以看出，在模型识别的验证期内，每层观测孔的观测水位基本与计算出的水位一致，满足了模型识别和验证的要求。结果表明，水文地质概念模型是合理的，所建立的数值模型基本反映了模拟区域地下水运动规律，可用于水源保护区的划分。

表2 泉域补给区入渗系数分区统计表

图7 地下水位观测井分布示意图Fig.7 Schematic diagram of groundwater level observation well distribution

表3 各分区参数表

表4 各分区参数数据来源表

图8 不同时间的地下水位拟合图Fig.8 Groundwater level fitting map of different time

2.4 溶质运移模型(流线追踪)

利用上述识别、验证后的数值模型，以2018年1月1日作为预测的初始时刻，通过运行Visual Modflow模拟，并结合《中华人民共和国国家环境保护标准》关于地下水源地一、二级保护区划分的原则，模拟得到运移100、1 000 d的运动轨迹。

受地下水径流和渗透系数的影响，质点在羊角铺水源地西部和北部运移迹线较长，南部和东部相对较短。在羊角铺水源地北部，粒子运移100 d时北部最大运移半径约为1 120 m，西部约770 m，南部约960 m，东部约618 m；粒子运移1 000 d时北部最大运移半径约为3 230 m，西部约3 980 m，南部约2 215 m，东部约950 m。

3 保护区范围的初步确定

保护区范围初步确定其拐点坐标系采用2000国家大地坐标系(CGCS2000)。

3.1 一级保护区范围的确定

根据《饮用水水源保护区划分技术规范》(HJ 338—2018)有关规定，利用上述识别、验证后的数学模型，获取以羊角铺地下水开采井为中心溶质质点向外迁移的运动轨迹,根据溶质质点运移100 d的距离确定地下水水源地一级保护区范围。是由下列地理坐标点所连曲线包围的范围，面积2.52 km2，如表5和图9所示。即羊角铺水源地地下水开采井为中心距粒子运移100 d时北部最大运移半径约为1 120 m，西部约770 m，南部约960 m，东部约618 m，包括羊角铺村和三河底村辖区范围。

表5 一级保护区边界拐点坐标

图9 一级保护区示踪粒子运移轨迹及定界图Fig.9 Tracer particle migration trajectory and boundary diagram in the first class reserve

3.2 二级保护区范围的确定

根据数值模型预测，在划分完一级保护区后，根据溶质质点运移1 000 d的距离确定地下水水源地二级保护区范围。粒子运移1 000 d时北部最大运移半径约为3 230 m，西部约3 980 m，南部约2 215 m，东部约950 m。是由下列地理坐标点所连曲线包围的范围，面积为26.28 km2，北至山底村北300 m，南至豆腐沟村南220 m，西至宿凤村西560 m，东至S222省道，如表6和图10所示。

表6 二级保护区边界拐点坐标

图10 二级保护区示踪粒子运移轨迹及定界图Fig.10 Tracer particle migration trajectory and boundary diagram in the second class reserve

4 结论

(1)通过运用数值模拟的方法，建立研究区的水文地质概念模型、数学模型和数值模型，为减少因模型参数的不准确性带来的误差，当在模型水文地质参数分区时结合遥感影像来确定含水层参数，在借助Visual MODFLOW软件模拟地下水流场和溶质质点的运动轨迹，并根据在不同的时间标准内所迁移的距离来界定地下水水源地各级保护区的范围。

(2)以羊角铺水源地为实例，进行了数值模拟法的应用研究。在地下水水流模型的基础上，通过流线追踪技术，以溶质质点迁移100 d和1 000 d的距离分别作为一级保护区和二级保护区的边界，划定了水源地保护区。

(3)在最终确定保护区的边界时，应充分考虑水源地的实际地形情况，以便确定明显的标志物，使相关地下水管理部门更容易操作和制订管理措施。