基于有限元法的傍河水源地地下水资源评价及预测

雷 磊，龙子华，闫晓君，徐晓青

(1.中国电力工程顾问集团华北电力设计院有限公司，北京 100120；2.国网北京市电力公司，北京 100031)

模型计算指的是进行地下水流或溶质运移正反演计算，常用的方法主要是有限差分法、有限元法、边界元法等。FEFLOW软件是一套专门用于孔隙介质中的地下水流动的三维有限元数值模拟软件，近年来，该软件在地下水水资源评价得到了广泛应用。经多年的开采实践证明，傍河取水是保证长期稳定供水的有效途径，本文研究的某傍河电厂水源地已运行了约30年，周边地下水位下降情况比较严峻，已影响到电厂正常的生产用水。在对水源地水文地质条件现状再分析的基础上建立水文地质概念模型和数学模型，在此基础上运用地下水模拟软件FEFLOW建立研究区地下水流模型，并基于所建模型评价了当前地下水的水资源和进行了趋势分析。

1 研究区概况

1.1 自然地理

研究区位于宣化县城西北，为一条带状山间盆地，见图1。洋河自西向东从盆地流过(现基本已干涸)，将盆地地形形态进一步分为三部分，即沿洋河为河床、河漫滩及一、二级冲积阶地；洋河以北为冲洪积倾斜平原；以南主要为坡洪积台地。洋河阶地为上迭阶地，沿河断续分布，阶面较窄，高出河床5～10 m，冲洪积倾斜平原(冲洪积裙)和坡洪积台地前缘高出河床约5～20 m左右。电厂位于洋河北岸二级阶地，其配套水源地分二期，一期水源地位置在沙岭子镇至样台一带洋河北岸的河漫滩地段，二期水源地位于洋河南岸的西前所村至古树营村一带洋河河漫滩地中。

1.2 水文地质条件

研究区地貌成因属于构造侵蚀堆积盆地区，含水层岩性为中细砂、圆砾和卵石与灰黄色和褐黄色含砾粘性土互层，最厚可达160 m，为地下水赋存了良好的条件，地下水类型为弱承压水～潜水。沿洋河两岸地区富水程度极强，单位涌水量在65～90(m3/h·m)；沿洋河向南、向北富水性逐渐变弱，基岩山区由于基岩岩性及裂隙不发育，富水程度较差；地下水溶解性总固体含量基本上小于1 g/l，地下水水化学类型一般为HCO3－Ca－Mg型水，个别为HCO3－Ca－Mg－Na型水。水源地运行初期，研究区地下水的主要补给来源是洋河地表水，河床一个水文年平均渗漏量为1.154 m3/s，区域内潜水水位埋深一般0～3 m，经过十几年地下水的开采，研究区水文地质条件已发生了变化。现洋河基本断流，洋河已改为景观河，河底和两侧都铺设水泥石板，阻断了地表水渗透补给地下水，现区域内潜水水位埋深一般大于20 m，地下水补给现在主要为地下水侧向径流，其次为田间灌溉及大气降水入渗补给；地下水主要以侧向径流、人工开采等方式排泄。地下水流向自西北流向东南。

2 水文地质概念模型

2.1 计算区范围的确定

从地下水流动理论出发，第四系孔隙含水系统渗流场数值模拟的范围应取至流动系统的自然边界，即基岩与第四系的交界处。本次地下水资源评价的模型范围东北与西南边界为地表自然分水岭，西北边界为清水河和洋河的交汇处，东南边界为盘长河与洋河的交汇处，模型区计算范围与研究区面积一致，见图1。

图1 模拟区及计算区边界

2.2 边界性质的概化

根据研究区地形地貌、水文地质及工程地质条件等，可将研究区的边界概化见图2。

BC、DA边界：概化为零通量边界，此两处边界为地表自然分水岭，与地下水流场的流线平行，为隔水边界。

CD边界：概化为流量边界(第二类排泄边界)，即在水源地的影响范围之外，地下水的流动仍保持原来的流动状态，向外流出，但由于开采量增大，会袭夺一部分水量流入井内，从而使侧向流出减少。

AB边界：概化为流量边界(第二类补给边界)，补给来源主要为上游地下水的侧向流入。

在垂向上模型的上边界为潜水水面，接受大气降水、田间灌溉渗漏入渗补给，以人工开采、大气蒸发来排泄地下水。下边界取第四系单一结构潜水含水层的下界面为相对隔水层。

图2 水文地质概念模型示意图

2.3 含水层的概化

研究区内有供水意义的含水层主要为第四系全新统和上更新统松散岩类孔隙水，含水层岩性主要为卵砾石和中细砂，含水层之间并无连续稳定的隔水层，相互间水力联系密切，并且区域内地下水皆为统一开采，故可将计算区含水层概化为非均质各向异性的统一的潜水含水层。本次模拟模型的底部边界根据岩性和第四系含水层的厚度以及下部渗透系数较差的粘土层厚度确定，模型的顶部根据实测地面标高插值得出，三维建模见图3。

图3 含水层顶底板标高示意图

3 地下水数值模拟

3.1 数学模型

根据计算区水文地质概念模型，对应的数学模型选用一层非均质各向异性三维非稳定流数学模型，所建立的数学模型可表示见式(1)。

式中：Ω 为模拟区域；L2；Kxx、Kyy和Kzz分别为X、Y和Z方向的渗透系数(L/T)，Kxx=Kyy；H为水头值(L)；H0为初始地下水位；ε为源汇项(L/T)；S为给水度；n为边界面的外法线方向；Г为AB、CD侧边界；B为底边界和BC、DA侧边界。

3.2 网格剖分

对计算区在空间上的离散包括平面上的网格剖分和垂向上的分层。本次研究区的面积为230 km2，研究区剖分的单元格共有57005个，节点57240，剖分单元格为三角形，平均边长约120 m，面积约6000m2。在电厂水源地和自来水厂水源地区域进行了网格加密，加密三角形平均单元格边长约0.2 m，如图4 所示，垂向上概化为一层。

图4 模型网格剖分图

3.3 边界条件

将底边界处理为隔水边界；上边界作为开放边界，考虑入渗及蒸发；将侧边界条件概化为流量边界和零通量边界。

3.4 参数性质的概化

模型中每个单元格的含水层的厚度、地形高程及底板高程由ARCGIS软件自动插值获得。根据渗流场的特征，将研究区划分为5个参数区，每个参数区根据利用已有抽水试验资料计算渗透系数和给水度，结合不同岩性的经验值，确定不同岩性的水平、垂向渗透系数和给水度，模型拟合、核正后参数见表1。模型参数分区及渗透系数赋值见图2(a)。降雨入渗系数依据地形地貌和出露岩性确定。

表1 模拟模型核正参数

3.5 源汇项的处理

田间入渗补给量则是分不同灌区，根据各灌区的地下水位埋深、包气带岩性和灌溉方式及时间等因素影响，确定回渗系数，然后根据收集到各乡镇的相关农业用水量计算各单元的农田灌溉回渗量，按各月的用水比例进行分摊，将其概化为单位面状补给量。

地下水开采用well程序包模拟。区内现状地下水开采以电厂水源地生产用水、农业开采用水和自来水厂开采用水为主，本次模拟根据实际调查的结果进行分配。其中，电厂水源地生产用水和自来水厂开采用水均按井边界条件赋值。农业开采用水由于开采井数量大，且开采具有季节性，在每年4～5月、7～8月为灌溉期，为简化农业开采赋值，又不失准确性，在本次模型中，我们将农业开采量总额比农业开采井分布面积值，在模型中以面源来进行赋值，其值为已扣去灌溉入渗补给量的大小。

潜水蒸发量根据多年平均蒸发量以及模拟区内的水位埋深分布情况将其分配到各个单元格上。

洋河地表水对地下水的补给量在本模型中可不考虑。由于评价区内地表水体(洋河)河道已改为景观河道，河槽下及两侧均铺设水泥石板，造成地表水与地下水联系不密切，且洋河上游建设有橡皮坝，造成洋河常年断流，地表水入渗量极小，可忽略。

3.6 初始水位

采用2014年1月实测地下水等水位线为初始流场，对模拟期后期2016年10月的流场进行拟合。模型按非稳定流运行，通过初始流场的拟合，在调整参数的同时，可以判断边界处理的合理性，经反复调试计算，流场拟合结果见图5。

图5 流场拟合图

4 模型识别与校正

4.1 水位拟合检验

为了验证对地下水系统模拟是正确的，研究区内设有3个水位观测点(图5)，分别位于洋河北岸和南岸，分布较均匀，具有较好的代表性。通过不断地改变水文地质参数，重复计算，直到整个模拟期内观测水位与计算水位的水位差在精度要求的范围内。

从图5可以看出，观测水位和计算水位在计算区范围内变化趋势总体保持一致，即上升和下降趋势一致。依据计算区2014年1月至2016年10月地下水观测水位资料，时间步长以天为单位。研究区内三个观测孔的拟合水位过程如图6。从图6可以看出，整体曲线拟合很好，通过以上模型的识别和校验，所建立的地下水数学模型能够较为真实地反映研究区地下水系统的实际情况。

图6 水位拟合图

4.2 基于模型的地下水均衡

FEFLOW软件提供了记录计算过程的输出文件，本次模拟的总补给量为1.35×104m3/d，总排泄量为1.85×104m3/d，均衡差为－0.50×104m3/d，呈现负均衡，这与该地区地下水位逐年下降的趋势吻合。

4.3 研究区地下水水资源的趋势分析

为了评价研究区水资源今后趋势变化，我们选取当地1996～2005年这10年降雨量作为本次模型降水量的输入值。因为这10年不仅包括了丰水期、平水期及枯水期不同年份的组合，而且年平均降水量偏小，较能体现今后降水量减少的趋势。本文从两种情况来对比分析，一种是保持现开采规模不变的情况，另一种是电厂水源地停采地下水采用城市中水作为生产用水或部分采用中水部分采用地下水的情况。10年后，其水位预测见图7、图8。

图7 开采规模不变情况10年后水位预测图

图8 电厂水源地停止开采情况10年后水位预测图

从水位预测图上可以看出，如果研究区开采规模不变，10年后，研究区水位下降了约12 m，相当于水位平均下降1.2 m/a。如果电厂停采地下水，采用城市中水作为生产用水，10年后，水位下降了约8 m，相当于水位平均下降0.8 m/a。说明，电厂采用中水，缓解了水位下降的速率，但并没有改变研究区水位持续下降。这是因为研究区内人工开采为地下水的主要排泄方式，其中，电厂生产用水只占用了其中之一，约17%；农业灌溉用水比例较大，灌溉持续用水依然会导致区域水位下降。但电厂采用中水后，评价区内水源地周边漏斗有所恢复，水位下降速率变缓，对整个研究区的地下水环境也有所改善。

在评价区内，以AS－56观测孔为监测点，对比两种情况下，水位变幅情况，见图9。从图中，我们可以看出电厂采用城市中水作为生产用水后，水位下降较正常情况下开采，下降幅度小很多。在0～900天范围内，年水位变幅基本持平，主要是因为这两年为丰水年，年降雨量在466.7 mm，给当地水资源有了很大补给，之后年份为平水年和枯水年，降雨量较丰水年份有所减少，水位又开始下降。

图9 两种情况下AS-56观测孔内水位变幅图

5 结论

通过对电厂傍河水源地的水文地质条件进行再分析，建立概念模型，依据水文地质概念模型对应的数学模型选用一层非均质各向异性三维非稳定流数学模型。经过模型的识别和校验，所建立的地下水数学模型能够较为真实地反映研究区地下水系统的实际情况。

本次模拟的总补给量为1.35×104m3/d，总排泄量为1.85×104m3/d，均衡差为－0.50×104m3/d，呈现负均衡。通过对两种开采情况地下水水资源趋势变化可以看出，如果保持现开采量不变，10年后地下水水位持续下降，由于将来研究区内农业的持续发展，水位的下降速率可能比预测更加严峻，电厂生产用水更加得不到保障。因此，从电厂安全用水考虑，电厂将来可采用城市中水，满足生产用水。并且，根据电厂采用城市中水预测的结果，电厂减少开采地下水，水位下降的速率变缓，遇到丰水期时，地下水水位还能恢复，可有效地对区域内地下水环境进行改善。

[1]丁继红，周德亮，马生忠.国外地下水模拟软件的发展现状与趋势[J].勘查科学技术，2002，(1).

[2]刘国义，等.傍河取水水源地地下水资源评价及实例研究[J].安全与环境工程，2008，15(2).

[3]郇心善，等.沙岭子电厂水源地供水水文地质勘测报告(详勘阶段)[R].北京：水利电力部电力规划设计院，1987.

[4]马小波.基于有限差分法的康滩水源地地下水资源评价[J].地下水，2016，38(6).