GMS在滦河三角洲浅层地下水数值模拟中的应用

张佳宁

（中国市政工程东北设计研究总院有限公司，长春 130021）

海水入侵是指由于自然因素或人为作用，破坏了沿海地区地下水中的淡水与海水中的咸水之间的平衡状态，海水（含盐量较大或矿化度较高）沿着与其具有水力联系的含有淡水的地下水含水层向内陆方向扩散，使得原淡水含盐量升高的过程[1]。导致海水入侵现象发生应具备三个重要条件，前提条件是区域处于干旱或半干旱地区，降水较少，水资源量不足，存在潮汐与潮流；基础条件是淡水含水层与海水之间存在水力联系，决定性因素是人为导致的地下水超采或截断含水层补给源，三个因素综合起来，使得滨海地区海水入侵范围逐步扩大[2]。

1 研究区地质概况

研究区位于燕山山前自北向南倾斜的冲洪积平原之上，地形起伏较小。地面绝对高程在10～50 m。环渤海海岸带，滦河—洋河冲洪积平原及滨海冲海积平原组成了研究区的主要地貌，由海积和海积-冲积共同堆积共同作用形成了滦河—洋河冲洪积平原，地面高程在2～5 m。在河相沉积及海相沉积共同作用下，大多低洼地区排水不畅，加上海潮托顶，该区域盐渍化较为普遍，涝灾频发。海积平原绝对高程在1～3 m，地形起伏较小，坡比总体小于0.1‰。海积冲积平原地处上述两种地貌单元之间，地面绝对高程3～5 m，主要不良地质现象为土壤盐渍化和沼泽化。根据不同岩性，海岸带可以分为淤泥质、砂质及基岩海岸等类型。海岸带地表出露岩性自东向西逐渐由细粒过度到粗粒。涧河口—大清河口多为淤泥质海岸，地势极低，岸线外潮间带发育，并残留一些砂岛，盐田、养虾池遍布[3]。

本区基底为太古界和下元古界变质岩系。第四纪沉积物主要为滦河水系建造，由一套陆相碎屑沉积物组成，滨海地区夹有2～6层海相层。Q1～Q4地层连续沉积，没有明显的沉积间断。

评价区地下水主要由潜水及岩溶水组成，其中潜水主要赋存于第四系松散孔隙之中，岩溶水含水介质为奥陶系灰岩，由于研究区内主要工业、农业及生活用水来自于第四系松散地下水的开采，因此第四系含水层为此次研究的目标含水层。研究区主要由滦河补给的地下水系统，东至渤海西岸，西部主要通过黏性土组成的弱透水层与其下白河-蓟运河作为主要补给源的含水层在空间上进行了相接。根据地下水含水层内水质不同，又可分为冲洪积淡水系统与冲海积咸水系统。

通常情况下，相对不透水层将上部孔隙潜水与下部岩溶水在空间上进行了分割，二者之间几乎不存在水力联系。但是，随着人类活动的日益增多，研究区地下水得到广泛利用（广泛采用通天滤水管），致使天然隔水层作用减弱，发生兼之复杂的岩相变化，导致上下两含水层水力联系较为密切。总体上，大气降水为研究区浅层地下水的主要补给来源，深层地下水主要为浅层地下水越流补给。

研究区地下水主要补给方式有大气降水补给、河流侧渗补给、地下水侧向径流补给及人工灌溉回渗补给。天然条件下地下水流场主要受研究区地形地貌影响，但研究区长期大面积人工开采地下水使得流场形态变化较大，第四系松散孔隙潜水主要渗流方向为自北向南，从北部山区指向南部滨海平原，渗流方向于区内滦河方向近似平行，水力梯度不大于2‰。在地形坡度与地下水赋存条件双重影响下，南北两侧地下水水质及其开发利用程度有很大差异。模拟区北部属低山丘陵区，含水介质以圆角砾为主，颗粒较粗，地形坡度起伏较大，径流条件较好，水质较好，人为开发利用程度相对较高；南部属于滨海平原区，含水层介质以粉细砂为主，地势起伏较小，径流滞缓条件极差，水质较差，由于咸水原因，开发利用程度较低。

根据地下水含盐量，可将浅层地下水划分为全淡水区和咸水区。全淡水区内孔隙潜水为研究区内生活饮用水主要来源，也是农业灌溉的主要开采层，开采量较其他含水层大。随着经济发展，深层岩溶水开采程度逐渐加大，第四系松散孔隙潜水水位逐渐降低，开采强度也日益增加，导致研究区内局部出现小范围的降落漏斗。在全咸水地区内，咸水、卤水在地下水中所占比例较大，开发利用程度较差。受地形地貌及含水介质渗透性的影响，南北部浅层地下水侧向径流也存在较大差异。

2 水文地质模型

2.1 文地质概念模型

根据岩性及透水性不同，将研究区内松散岩层划分四层，自上而下依次为四个相对独立的含水层。各含水层组之间存在相对不透水的黏土层，其分布较为分散，不连续，局部存在天窗，因此其无法构成完整的隔水层。本次数值模拟采用了准三维非稳定地下水流模型。模型北边界是根据项目工作区范围划定的边界，垂直于区域多年平均等水位线，与边界流线相重合，设定为二类零流量边界。模型东边界定在滦河西岸。模型南部边界划定为渤海北部海岸，海水与内陆地下水含水层之间存在水量交换微弱的黏土层，将其概化为一类水头边界。根据此次研究的范围，人为划定了模型的西部边界，由于该区域资料较为详细，因此概化为二类零流量边界。滦河位于模拟区的西北侧，因此将其概化为一类水头边界。研究区下部为中风化基岩，其透水性极差，构成了此次模拟的下部边界，将其视为隔水底板，模型向下流量为零。选取潜水顶面为上部边界，其为开放型边界，参与研究区内水文循环。

模拟区内含水层的主要补给源主要有：以面状形势补给的大气降水入渗、以线状形式补给的河渠侧向渗漏、以面状形势补给的农田灌溉回渗等，由于模拟区内内河渠密布，纵横交错，无法将其清晰区分，因此此次模拟将其视为面状补给。

地下水流的数学模型如下：

式中：W为流入汇或源项的水量；Ss为储水率（承压含水层），（给水度）潜水含水层；Kh为水平透系数；Kv为垂直透系数。

对于该模型来说，式（2）为初始条件，式（3）为第一类边界条件，式（4）为第二类边界条件，式（5）为第三类边界条件。

2.2 地下水流数值模型及识别

此次模拟利用GMS的Grid frame模块采用六面体的矩形网格剖分法，将模拟区剖分为150行、150列、2层。根据评价区的主要岩性及地下水赋存条件，将浅层潜水含水层分为6个分区进行水文地质参数赋值，将深层潜水含水层分为3个分区进行水文地质参数赋值，这些经验数据主要来自资料收集，基本可以反映研究区的实际情况。滦河三角洲数值模拟采用的时间单位是天，拟合时间为4年。采用这种时段划分，模型计算水位是天平均值，与观测水位吻合。根据搜集相关区域的入渗系数，将模拟区按4个分区进行降水入渗补给系数的赋值。区内小型河流及渠系成网状密布，在区域性地下水流场中表现为面状渗漏补给，只能以面状入渗的方式进行赋值。此次模拟将灌溉较为集中的4-6月的总入渗量平均分配到各个时段，将其作为模拟的初始流量，本次模拟是在保持深层开采总量不变的条件下，对浅层及深层地下水开采量的空间位置进行了重新分配。

2.3 地下水流数值模型校正及水位预测

此次模拟选取多年水位观测资料中的一个较为完整水文年的观测数据，同时绘制时间-水位变化曲线图，通过分析可以发现研究区内部分水位观测孔内水位动态变化曲线与区域地下水补给及排泄变化曲线相关性较好，整体趋于一致，较好地反映了模拟区地下水的年内变化特征。但是，部分观测井受人为地下水开采影响，其地下水位动态变化曲线与区域地下水补给及排泄变化曲线相关性较差，局部出现矛盾，不能用以反映地下水年内变化特征，对于受附近开采井影响的水位观测孔，模型拟合时在观测孔附近设置了抽水井。

结合研究区近期、中期、远期的地下水水位预测结果，笔者发现，研究区地下水流场基本稳定，由于该区域地下水补排状态呈负均衡，将导致降落漏斗的中心水位会进一步加大降深，降落漏斗面积亦呈小幅度增长，将会导致研究区由于区域降落漏斗引发的环境地质问题进一步加剧，因此需针对这一问题，在水资源利用方面做出合理分配。

3 控制海水入侵问题的取水方案优化

据调查，研究区内主要河流等均发育在断裂之上，各条河流直接入海，在其上游均修有水库，从而减少了入海河水水量和泥砂量。另外，河流的下游入海河段，河床低而平缓，在涨潮时海水沿河道上溯，尤其是遇有天文大潮和风暴潮时，其上溯距离更远（一般上溯3～4 km，最远出现在洋河可达10 km），充满河床的海水与河水之混合水沿河床周边渗漏补给地下水，使沿河两岸地下水受到污染。因此，近海松散的砂砾质含水层及构造破碎带与通海河流的砂砾质河床是本区海水入侵的主要通道[4]。咸水分布区主要分布于研究区东南及南部，面积约3511 km2，从区域上看，自东北向西南，地下水的底板埋深同咸水体的厚度均呈加大的趋势，而沿海地区分布最广，其厚度自北向南、自西北向东南逐渐增大。将2015年年末水位与底面沉降速率等值线图及地面沉降等值线图进行配准，得出海水入侵范围与地下水位降落漏斗分布特征。据分析，评价区东南部地区海水入侵情况相对较为严重，咸水体的分布位置与区内的浅层地下水集中开采中心基本吻合。据此可推测，研究区南部地下水降落漏斗为形成该区海水入侵的主要诱因，通过调整地下水开采方案可以达到控制海水入侵范围的目的。

本次地下水开采方案的优化以海水入侵范围缩小或保持现有海水入侵范围不再进一步扩大为目的，通过对海水入侵区域的空间演化进行分析，以控制研究区南部浅层地下水降落漏斗降深。在减少现有地下水集中开采区地下水开采量的同时，在研究区中部及西部地下水未超采地区建立分散型地下水开采水源地。

从以上优化方案模拟的结果来看，2020年即优化方案实施5年间原地下水位降落漏斗面积和中心水位降深变小，地下水流场开始变化。至2030年原区域降落漏斗基本消失，形成新的面积及中心水位降深较小的多个小型降落漏斗，基本形成内陆地区地下水位高于海水位，地下水流向由原来的海水向浅层降落漏斗流动，基本改为内陆地区向海洋排泄地下水。所以，该方案在改变现阶段由于地下水超采造成的区域降落漏斗问题上是可行的。

图1 研究区2020年、2030年浅层地下水水位预测图

4 结论

本文通过数值模拟方法，应用GMS地下水流数值模拟软件，针对冀东地区目前存在的与水资源开发利用有关的海水入侵这一环境地质问题，通过模拟不同的水资源利用方案，得到了有利于改善当前地下水超采问题的水资源调度方案。但是，本次模拟仍然存在一些问题。

本次模拟中监测孔内地下水水位变幅计算值整体小于地下水水位变化的实际值，未能准确反映研究区黏性土释水、储水系数大的特点。此次模拟将模型西侧概化为一类水头边界，属于静态边界，不能很好地反映模拟区内水位动态变化特征；本次地下水数值模拟采用准三维地下水流模型，计算值可以近似代替各含水层之间的越流补给量，但不能有效刻画相对不透水层内的水头分布特征。准三维模型虽然可以很好地反映该区域浅层含水层的水位变化，但是浅层对深层的越流补给、越流及压缩释水过程都不能得到较为准确的刻画；当前收集的水文地质资料，主要针对各主要含水层，对相对不透水的黏土层研究较少。由于目前资料有限，难以满足高精度的模拟要求，在模型反演校正过程中，模型对水平渗透系数的识别精度不高。