基于GMS的长春齐家地下水水源地二级保护区划分

鲍新华, 于瀚博, 计 量

1.吉林大学新能源与环境学院，长春 130021 2.中国电力工程顾问集团东北电力设计院有限公司，长春 130021

0 引言

水源地是指提供城镇居民生活及公共服务用水取水工程的水源地域，分为地表水水源地和地下水水源地，包括河流、湖泊、水库、地下水等。为保护我国居民正常生产生活用水安全，建立水源地保护区是一种成本低廉、效果显著的措施[1-4 ]。常用的水源地保护区划分方法有经验值法、经验公式法及数值法等[5-8]。经验值法和经验公式法总体而言成本较低、操作便捷，但由于这两种方法一般采用固定的参数值模糊化处理水文地质特征，导致保护区划分得不够精确；而数值法拥有结果精确、能够对含水层进行细致刻画的优势，因此本文主要采用数值法，通过建立地下水流及水质模型进行保护区划分。现有的采用数值法对地下水保护区进行划分的研究主要包括数值法建立地下水水流模型、利用粒子示踪法来划分保护区范围[9-11]等，如韩京龙等[12]采用数值法对磐石市下水水源地保护区进行了划分，其模型能够较为准确地描述研究区含水层，提出的保护区划分方案实用性较高[13-14]。

齐家水源地是长春市重要的后备水源地和当地农业生活用水的主要来源。由于其临近石油开采区，随着油田前期勘探试采工作的进行，原有的水源地一级保护区已经无法解决油田发展与用水安全之间的矛盾。为了保护周边居民用水环境不受石油污染[1-4]的同时能兼顾油田开采工作的持续有效进行，在一级保护区的基础上划分范围更大的二级保护区就显得尤为重要。由于研究区临近石油开采区，传统的利用粒子示踪法建立的保护区并不适用；且由于污染质受分子扩散作用影响，在其随地下水渗流迁移过程中迁移速度要快于地下水流速，因此仅采用水流数值模拟并不足以准确描述污染物的运移范围。有别于前人以水源地为中心进行溶质运移模拟确定保护区边界会带来的不够细化等问题，本文在地下水水流模型的基础上，在本研究区采取假定污染源位置反向模拟溶质到达水源井的时间，以期更精准地确定保护区边界[15-21]。

1 研究区概况

齐家饮用水水源地研究区面积约260 km2，位于长春市双阳区东北部，北临石头口门水库，南距双阳城区3 km，西北距长春市区45 km，双阳河自南向北流经研究区，如图1所示。

研究区北部平原区地下水资源较为丰富，南部低山丘陵区水资源较贫乏。根据地下水赋存条件和水力特点，区内地下水分为5种类型：第四系松散岩类孔隙潜水，主要沿饮马河、双阳河分布；新近系碎屑岩类裂隙孔隙承压水，主要分布在双阳区北部平原区；白垩系碎屑岩类裂隙孔隙承压水，主要分布在饮马河、双阳河河谷地带；碳酸盐溶隙水，主要分布在低山丘陵地区；基岩裂隙水，主要分布在双阳盆地东西两侧及南侧。

研究区共有集中式饮用水水源地供水井30眼，井深200～220 m，一般单井出水量约40 m3/h，总计服务人口约10万人。作为长春市第二大地下水水源地，齐家水源地主要为周围居民提供生活用水及作为长春市后备应急供水水源地，可保证长春及周边地区供水安全[22]。为应对缓解我国能源危机，1981年对研究区周边进行全区勘察工作，探明齐家镇周边地区石油地质储量1 017×104t[23]，具有一定的开采价值。然而石油开采井的建设会对周围地下水水质安全形成一定威胁，保障居民用水安全与石油开采之间的矛盾十分突出，因此科学建立地下水水源地保护区的需求日益迫切[24-26]。2017年5月，长春市环保局双阳分局根据《饮用水水源保护区划分技术规范》(HJ/T 338—2007)[22]中的相关技术指标对研究区进行了水源地一级保护区的划分工作，保护区范围如图2所示。新规范(HJ/T 338—2018)[6]施行后，关于一级保护区的规定与(HJ/T 338—2007)相同。

图1 研究区地理位置图

图2 研究区水源地一级保护区范围

一级保护区的确定，为齐家水源地水质保护提供了基本的技术支撑，但由于其范围过小，不足以更全面地保护水源地。本文拟在此基础上，进一步对二级保护区进行研究划分。

2 水文地质条件概化

通过研究区的地层结构、水动力场、水化学场分析，对研究区进行水文地质条件概化。从采油区到水源地附近的区域为第四系潜水和新近系承压水的三维流场(图3)。新近系地层中有多套含水层与相对弱透水层互层分布，由于地层情况较为复杂，已有水文地质资料不足，将其概化为一套含水层。模拟中渗透系数采用了砂砾石参数作为概化后的参数，概化后局部地区会出现溶质沿渗透性较差的介质运移距离较近的情况，在划分保护区时应将砂砾石参数考虑在内进行适当调整；但宏观层面上对整体保护区划分影响不大，因此可将新近系地层中含水层概化为一套含水层。

2.1 确定模型模拟区

模拟区应尽量包含油田开发可能影响到的区域。根据研究区地下水的流向，确定可能影响的区域为下游齐家水源地的准保护区。

本次模拟中，区内河流较多，容易通过地表高程寻找相对独立的地下水分水岭。现有水位观测点绘制的地下水流场分布图表明了地表与地下水分水岭具有一致性；因此北部和西北部边界扩大到伊通河与双阳河分水岭，西南部以双阳河为界，东部以双阳河与饮马河的分水岭为界。如图4所示。

2.2 含水层概化

研究区地下水类型主要为第四系松散岩类孔隙潜水和新近系碎屑岩类裂隙孔隙承压水。

第四系松散岩类孔隙潜水含水层包括河谷一级和二级阶地含水岩组，地势低平。含水层岩性具有典型的二元结构特征，上部颗粒较细，主要为亚砂土、亚黏土等，下部颗粒较粗，主要为含砂砾卵石等，厚度总体由南向北渐增，颗粒由粗渐细，水位埋深由深渐浅。其含水层总体埋藏浅，厚度大，透水性较好，有利于邻区地下水的径流补给。地表岩性多为亚黏土，局部为亚砂土。区内水稻田多，灌溉回渗水可直接补给地下水。河流直接与第四系潜水含水层相接，地下水与地表水水力联系密切，水交替频繁，河水常年或间歇性补给地下水[27]。

第四系之下为新近系碎屑岩类裂隙孔隙承压水。含水层岩性为粉细砂岩、中粗砂岩、砂砾岩等，顶板埋深20～50 m。承压水接收上部孔隙潜水垂向和相邻地下水的侧向径流补给。

通过对研究区地质和水文地质条件分析，将研究区概化为非均质各向同性含水系统，渗流按三维达西流研究。

3 地下水流数值模拟

3.1 地下水流数学模型

根据研究区地质条件与水文地质条件，地下水流数学模型为：

(1)

根据文献[25]修改。

图4 模拟区范围及剖分图

式中：Kx，Ky，Kz分别为x,y,z方向渗透系数(m/d)；h为地下水水头(m)；W为源汇项(1/d)；μ为贮水率(1/ m)；t为时间(d)；h0为初始水头(m)；h1为一类边界水头(m)；Σ1为一类边界；D为研究区。

3.2 研究区剖分

根据研究区水文地质条件和渗流场特征，将模拟区垂向上剖分为3层，分别为第四系潜水含水层、粉质黏土相对弱透水层和新近系承压含水层。水平方向的网格大小为300 m×300 m，在油井区和水源地等主要模拟区附近加密到50 m×50 m。如图4所示。

3.3 边界处理

根据前期收集的水文地质资料，将研究区西北部伊通河与双阳河分水岭[27]概化为零流量边界，将西南部双阳河概化为水头边界，将东部双阳河与饮马河分水岭概化为零流量边界。

3.4 源汇项处理

模拟区潜水含水层主要补给为降水入渗补给、灌溉回渗水补给和河流侧向补给，主要排泄有蒸发排泄、人工开采和侧向排泄。承压水的补给主要为上部孔隙潜水的垂向补给，排泄主要为人工开采。降水入渗系数取0.17～0.23，稻田灌溉回渗强度取0.000 8 mm/d，最大蒸发强度埋深为1.9 m，蒸发极限埋深为4.25 m。

3.5 水文地质参数

模拟用到的参数主要有渗透系数、潜水含水层的给水度、承压含水层的贮水系数和降雨入渗系数等。各参数初值先通过抽水实验、土工试验、水文地质手册和经验值综合给定，再通过模型参数调整，最终获得模拟采用的水文地质参数。

3.6 模型的识别和验证

使用数值法软件GMS(groundwater model system)中的MODFLOW模块建立地下水流数值模型后，利用给定水文地质参数、边界和各均衡项数值初值进行模拟，得到概化后的水文地质模型地下水流场空间分布。通过2017—2018年丰水期、枯水期的流场，识别验证水文地质参数、边界值和各均衡项，以便建立的模型更好地拟合研究区的水文地质条件。模型识别拟合主要原则为：识别后的水文地质参数与其他资料获取的水文地质参数基本一致；模拟的地下水流场与流场基本一致；边界条件、各水均衡项与实际情况基本一致。

模拟采用2017—2018年丰水期的观测水位进行模型识别，取枯水期的观测水位进行验证。通过反复识别和验证，最终确定上部第四系潜水含水层的渗透系数为6.0～18.0 m/d、给水度为0.1～0.2，中间层粉质黏土相对弱透水层的渗透系数为0.1～0.2 m/d，下部新近系承压含水层的渗透系数为2.0～14.0 m/d、弹性释水系数为0.003～0.008，其中渗透系数(K)分区见图5。图5中参数区范围与图4不完全一致(西北侧)，主要是参数区为水源地所在区域，资料较多；而西北侧资料有限，且考虑采用伊通河与双阳河分水岭边界较好，故西北侧所用参数参照邻区给出。图6为水源地所在区域上模型计算丰水期流场与实测流场的对比，拟合情况较好，水位相对误差为5%～10%，说明模型建立基本合理。

采用2017—2018年枯水期的观测水位进行模型验证，图7为模型计算的水源地所在区域上枯水期流场与实测流场的对比。计算结果表明：相对误差在1%～5%之间，拟合情况较好；地下水流动方向基本与河流方向一致，并且在主要的饮用水开采区形成明显的降落漏斗，地下水流速明显加快，可能会提高污染物迁移速率，使其迅速进入饮水井，从而对周边居民饮用水安全造成威胁。选取2018年1月—2019年12月，对代表性观测井W18井和齐家观测井潜水水位变化计算值与实测值进行拟合观察，上述流场拟合检验结果如图8所示，结果表明水位拟合相对误差均不超过5%。经过模型验证，表明本次模拟具有较高的精度和可信性。

图5 模拟区潜水含水层(a)和承压含水层(b)渗透系数分区

4 溶质运移模拟

以上述地下水流数值模型为基础，通过使用GMS中MT3DMS模块建立溶质运移模型来进行饮用水水源地保护区圈定。

4.1 地下水溶质运移模型

4.1.1 地下水溶质运移数学模型

地下水溶质运移的对流-弥散方程为

(2)

4.1.2 评价目的层

油田开采过程中发生泄漏，发生泄漏的位置可能是油井的各个层位。研究区的主要地下水类型是第四系孔隙潜水和新近系碎屑岩类裂隙孔隙承压水。如果泄露事故发生在第四系亚黏土层或者新近系泥岩层，因为这两个层位的渗透性较差，污水的泄漏量将会很小，不会对齐家水源地构成威胁；如果泄露发生在第四系潜水含水层，采油井在一定的压力下可向潜水含水层释放一定的污水，污水有可能通过由南向北方向的水平运移，将污染物运移到齐家水源地所开采的位置，再通过垂向渗透对水源地的承压含水层水源造成影响；如果泄露发生在新近系承压含水层，采油井在一定的压力下会向新近系含水层释放较多的污水，污水有可能通过由南向北方向的水平运移，将污染物运移到齐家水源地所开采的位置，对水源地的新近系承压水水源造成污染。

综上，评价目的层拟选择发生泄漏事故造成污染较为严重的第四系潜水含水层和新近系承压含水层，以齐家水源地开采的新近系承压含水层为主。

4.1.3 确定污染运移模型参数

根据笔者前期弥散试验结果得知，研究区纵向弥散度为0.560 m，横向弥散度为0.056 m ，垂向弥散度为0.056 m。模拟层主要为中砂，取有效孔隙度为0.3。

图6 模拟区丰水期潜水水位(a)和承压水位(b)实测值与计算值拟合图

4.2 模拟过程

石油类中不易挥发的大分子组分与砂的吸附关系多为线性吸附关系；但考虑到研究区主要含水层为砂砾层，为保守起见，模拟预测中并没有考虑到这种吸附的阻挡作用。石油类液体本身的黏滞系数较大，在含水层介质中不易流动，但考虑到评价的保守性，在模型中可不加以计算。故模型可采用保守型污染质计算方法进行泄漏事故发生后的污染预测。模拟1 000 d后污染物在承压含水层中分布见图9。

由图9可见，如果位于水源地内的采油井发生泄漏，对水源地的影响将是明显的。被污染的地下水将很快进入水源井污染饮用水水质，开采井附近受开采影响，具有较大的水力坡度，使得地下水流速加快，污染质运移速度也将加快。根据数值模拟结果，一旦发生泄漏事故，石油类污染质在地下水中经过 1 000 d的运移距离在400～800 m之间。

根据图9，在水源地周边距水源井400～800 m范围内布置多个假想的污染点源，污染点源为假想石油开采井，石油开采井穿过承压含水层，污染物通过套管部分渗漏入含水层。根据研究区石油开采井特点，采用5%作为采出液渗漏比，利用已有的开采井年开采量及采出液含水率资料可计算出模拟年份污染源向地下水渗漏石油污染物速率为0.25 m3/d、质量浓度为6.2×104mg/L，以此为源强进行预测。

图7 模拟区枯水期潜水水位(a)和承压水位(b)实测值与计算值拟合图

图8 研究区W18井(a)和齐家观测井(b)水位拟合图

图9 石油开采井泄漏污染物1 000 d后在承压含水层分布

利用先前建立的溶质运移模型进行模拟，预测污染发生1 000 d后的最大污染边界。根据污染晕最大边界与水源井之间的位置关系，在先前划定的范围内移动污染源的位置，再次进行模拟，直至污染晕在1 000 d的最大边界刚好到达水源井。由于含水层污染难以治理，石油类污染物将长期被含水层介质吸附并缓慢释放，因此上述模拟得出的位于水源井上游的污染源位置应再向水源地以外移动100 m左右留作缓冲区，此时污染源的位置即为二级保护区的边界。污染源最终位置1 000 d污染晕范围如图10所示。

从模拟结果来看，1 000 d时明显沿地下水流方向污染晕扩散范围更大，证明石油类污染物在含水层中运移主要受地下水流的影响。此外污染物同时也向其他方向运移，说明污染物受水动力条件影响外，本身分子扩散也起作用。

最终规划齐家水源地二级保护区范围如图11所示，面积为33.5 km2(含一级保护区面积8.5 km2)。

由最终划定的水源地二级保护区可见，齐古1、昌25、昌31、昌37、昌88、昌90、昌91等部分规划石油开采井位于二级保护区之内，会对水源地水质安全构成直接威胁，建议取消规划、停止建设；昌36-4、昌39、昌84、昌87、昌105-11等规划石油开采井位于二级保护区边界附近，存在对水源地水质安全造成破坏的风险，建议将开采井规划位置向二级保护区外迁移500 m以上距离以降低石油勘探、开发过程中对地下水破坏的风险。

图10 1 000 d污染晕范围

5 讨论

文中采用的研究方法对划分一级保护区同样适用，文中涉及到的一级保护区范围采用长春市环保局相关资料，仅做模拟计算验证。本文主要对研究区地下水水源地二级保护区划分进行研究，研究中有关问题进一步说明或讨论如下：

1)主要模拟层为新近系，为多套含水层与相对弱透水层互层，模拟中概化为一套含水层。在分层细化资料缺失时一般采用在模拟中概化为一套含水层[28]。这种概化与实际情况有一定出入。污染物运移时，一般在透水性强的地层中运移快[29]。模拟中渗透系数采用了砂砾石参数作为概化后的参数，计算的模拟污染晕范围总体较实际结果可能偏大，这考虑到了污染物迁移的最大可能。这种概化在资料有限及解决主要污染物最大迁移问题研究中，精度应该是有保证的[29]。

2)水源地西北侧非自然边界，边界条件不好给定，研究将模拟区边界扩展到伊通河与双阳河分水岭天然边界(模拟区边界见图4)。但模拟区西北部资料较为有限，因此进行模型运算时相关水文地质参数等资料参照水源地邻近区给出，这与实际相比会有一定误差。而模拟的重点区在水源地区域，模型调参与校正主要看水源地所在区域。这样的处理对水源地区域污染问题研究影响较小。

图11 水源地二级保护区范围

根据图6—图7流场图，地下水流方向基本沿主径流方向由西南向东北，位于水源地上游的石油开采井对水源地影响更大，因此在划分水源地二级保护区时水源地上游划分范围应较下游更大，通过溶质运移模型确定的二级保护区范围也证实了这一点。本文最终划定水源地二级保护区符合这一基本判断。

6 结论与建议

1)在地下水水流模型的基础上使用MT3DMS模块建立研究区溶质运移模型，在距水源地400～800 m范围设定若干污染源，1 000 d污染晕最大边界距水源井30～80 m的污染源位置定为二级保护区边界，划定二级保护区面积为33.5 km2(含一级保护区)。

2)建议齐古1井等规划建于二级保护区范围内的7口油田开采井取消建设，昌36-4井等规划建于二级保护区边界的5口油田开采井向保护区边界外迁移500 m以上距离；并建立完善的地下水污染监控系统体系和应急响应措施，最大限度地降低石油勘探、开发过程中风险事故对地下水环境的影响。

3)建立在地下水水流模型基础上的溶质运移数值模型与解析法相比，能够更加精细地刻画溶质运移细节，进行更精细的水源地保护区划分。本文的方法对其他类似问题的研究有借鉴意义。但同时应注意到，模型精细的刻画程度主要取决于建模依据的地质-水文地质基础的正确认识。没有大量的基础数据支撑，仅靠模型剖分的细化无助于精度的提高。油田运行过程中，还应注重实际污染物运移监测工作，一方面检验本文的计算结果，另一方面在资料进一步积累后，可对模型再进一步修改完善。