浅析采用三维数值分析方法计算水库地下水壅高

裴 颖

（长安大学环境科学与工程学院 西安 710021）

1 概述

水库蓄水后，库区周围地下水相应壅高而接近或高于地面，会产生一些不良现象，如农田盐渍化、沼泽化及建筑物条件恶化、矿坑涌水等。水库蓄水浸没问题是水库工程勘察与评价的重要内容之一，是确定库区征地和移民迁建水位的重要依据。严重的水库浸没问题影响到水库正常蓄水位的选择，甚至影响到坝址的选择。浸没评价结论是水库设计、建设非常重要的因素和先决条件，如果水库蓄水后的地下水埋深值小于浸没地下水埋深临界值，就会产生浸没。

江巷水库位于江淮分水岭北部、滁州市定远县连江镇境内，其主要任务是供水和灌溉，兼顾防洪等综合利用。水库总库容为1.30亿m3，水库正常蓄水位43.0m。工程规模为大（2）型。洪水标准为100年一遇设计，2000年一遇校核，设计、校核洪水位分别为44.42m和45.50m。

江巷水库为典型的平原地区水库，库周地势较为低缓，易产生浸没问题，可研设计阶段需开展浸没问题研究，通过一系列勘探、试验、观测等方法进行浸没评价，判定水库蓄水后是否会产生浸没问题。作为重要判定指标的库岸地下水位壅高值可通过地下水动力学方法计算，主要方法有解析法、三维数值分析法等。本文主要通过三维数值分析法来计算水库地下水位壅高值，得出库岸地下水埋深，为水库的浸没评价提供依据。

2 三维渗流模型计算地下水壅高方法

2.1 数学模型

假定库区为非均质各向同性，三维地下水流非稳定运动的数学模型可表示为：

式中：Ω为渗流区域，H为含水层的水位标高，K为渗透系数，μ为给水度；ε为含水层的源汇项，H0为含水层的初始水位分布，Γ1为渗流区域的二类边界，包括承压含水层底部隔水边界和渗流区域的侧向流量或隔水边界，n为边界面的法向方向，Kn为边界面法向方向的渗透系数，q（x，y，z，t)定义为二类边界的单位面积流量，流入为正、流出为负，隔水边界为0。

2.2 有限元离散

对式（1）进行有限元离散，得到整个渗流区域的代数方程：

式中：[G]为总体渗透矩阵（传导矩阵），[P]为贮水矩阵，｛H｝为未知节点水头列阵，｛F｝为已知的右端项。

采用“隐式差分格式”，有：

整理后，得：

2.3 水文地质参数确定

库区地层较为简单，表层均为第四系上更新统（Q3a）l地层覆盖，以重粉质壤土、粉质粘土为主，下伏基岩为白垩系红层，地表无出露。在库区郭小圩、孟严、东方曹、齐岗、小庙陈布置了5个地质剖面（图1），

图1 库区地质剖面布置图

库区地下水类型主要为基岩裂隙水和孔隙水两大类。含水层位主要为下伏厚层红层碎屑岩层顶部风化带裂隙水和上覆粘性土孔隙潜水。地下水赋存条件很差，水量较为贫乏。

（1）渗透系数

根据地质资料表明，库区均为第四系土层覆盖，层厚15～30m，以重粉质壤土、粉质粘土为主，该土层渗透系数一般都在n×10-6～n×10-7cm/s之间。下伏基岩为泥质、钙质粉砂岩和砂质泥岩，亦不易形成渗透通道，不具备形成地下大规模强渗流的地质条件，可以当做隔水边界处理。另外，通过对库区的地质调查发现，库区内所有的冲沟都是由第四系土层组成，为微透水或极微透水土层，构成库盆的地层都是第四系上更新统（Q3al）的粘性土，地下不存在强渗透含水层。渗透系数取值参考注水试验结果取值（表1）。

表1 渗透系数及压力传导系数取值表

（2）给水度

库区的岩性主要为重粉质壤土、粉质粘土，给水度取值为5%。

（3）降水入渗补给系数

降水入渗补给系数α是指降水渗入量与降水总量的比值，α值的大小取决于地表土层的岩性和土层结构、地形坡度、植被覆盖以及降水量的大小和降水形式等，它是一个无量纲系数，其值变化于0～1之间，库区的年均降雨量约为934mm，主要岩性为重粉质壤土、粉质粘土，降水入渗补给系数取值为0.1。

（4）潜水蒸发系数

潜水蒸发系数主要年水面蒸发量、岩性和地下水位埋深有关。库区的多年年水面蒸发量为798mm，地下水位埋深为0.5～1.5m，主要岩性为重粉质壤土、粉质粘土，蒸发系数取值为0.18。潜水蒸发量主要与潜水位的埋深、包气带岩性、地表植被和气候等因素相关。根据前人的研究成果，一般认为水位埋深大于5m的地区潜水蒸发很小。库区地下水位埋深一般小于5m，因此要考虑蒸发影响。

2.4 水文地质概念模型

水文地质概念模型是在综合分析地下水系统的基础上，对评价区地质、含水层实际的边界条件、内部结构、渗透性质、水力特征和补给排泄等水文地质条件进行科学的综合、归纳和加工，从而对一个复杂的水文地质实体进行概化，便于进行数学或者物理模拟。因此，建立水文地质概念模型主要应该考虑如下几个方面：概化后的模型应该具备反映库区水文地质原型的功能，概化后的各类边界条件应符合库区地下水流场特征，概化后的模型边界应该尽量利用自然边界，人为边界性质的确定应从不利因素考虑等。本次模型采用46m等高线作为模型边界，按照水文地质勘察钻孔所监测的地下水位插值获得边界水位（表2和图2）。

表2 库区观测孔地下水位统计表

图2 通过观测孔的地下水位插值得到的库区地下水位等值线图

模型四周为水头边界，在43m以下高程的也为水头边界，顶部为降雨入渗和蒸发边界，模型底部基岩为隔水边界，建立相应水文地质概念模型。

2.5 区域离散

计算区域以项目所在地中心位置为坐标原点，正北方向为y轴正向，正东方向为x轴正向，垂直向上为z轴正向，垂向上考虑1层，将库区域离散为51001个节点，57860个单元。

2.6 边界和初始条件

边界条件：库区为一个相对独立的水文地质单元，地下水由周边向库区流动，外边界为水头边界，通过实测地下水位插值获得；内边界分别为高程低于常水位42.5m和正常蓄水位43m的区域，上部边界为降雨入渗和蒸发边界，底部为隔水边界。初始条件：考虑稳定流计算。

3 库岸剖面地下水壅高计算

在库区选取了5个地质剖面计算水库蓄水后的地下水位变化，本文列出其中的两个地质剖面的计算结果（图3，其他剖面情况类似）。从图中可以看出，不同剖面地下水位存在一定差别，但总体上地下水位浸润线比较平缓。

图3 库岸剖面地下水位等值线图

4 计算结果分析

根据蓄水位为43m时的地下水位的分布情况，库区周边地下水位略高于库水位，地下水对库水位进行补给，其水位差不大，水力坡度较小。但大坝下游地下水等值线密集，水位降幅较大，水力梯度较大，表明大坝防渗性能较好。本次同时也计算了水库水位为常水位42.5m时的地下水位的分布情况，其地下水位分布规律与水库蓄水位为43m时基本一致。

三维数值分析法计算结果表明，与初始地下水位相比，蓄水后地下水壅高为0.09～1.53m。当水库蓄水位为42.5m时，库区周边区域地下水位高于常水位0.32～0.56m，平均值为0.41m。当水库蓄水位为43m时，库区周边区域地下水位高于常水位0.48～0.66m，平均值为0.52m。这是在库水位长期稳定在42.5m和43m时的地下水位变化，实际库水位一直在变化，因此，稳定库水位条件下计算的地下水位是最终的理想结果，比实际地下水位偏高。

从计算结果可以看出库岸地下水位的变化规律和趋势：水库运行多年后，从水库向四周，地下水位逐渐升高。当库水位上升时，与库水边距离越近，地下水位上升越快；距离越远，地下水位上升越慢。地下水位上升（壅高）较高的位置一般在库水边，越远离库水边上升越少。虽然不同区域土层参数有所差异，影响范围和水位上升高度也有差异，但地下水位上升趋势总体一致，地下水位浸润线比较平缓，库岸地下水位随库水位上升需要很长时间才能达到稳定。

5 结语

库岸地下水位壅高值计算主要采用地下水动力学方法，较常用的有解析法和三维数值分析法。两种方法各有其适用条件和特点。与解析法相比，三维数值分析法更适用于渗流场较复杂的浸没区，通过对可能产生浸没的区域内地下水位情况进行模拟，可以有效解决预测的浸没范围与实际浸没范围相比偏大的问题。对于规模等别较高的水库，可分别采用解析法和三维数值分析法进行浸没计算，对比、验证、分析计算结果，对更准确、合理地提出浸没评价结论是大有裨益的■