平原水库渗流场模拟及渗漏量计算——以斗门水库试验段为例

李文溢， 杨阿敏， 周维博

(1.长安大学 环境科学与工程学院， 陕西 西安 710054；2.长安大学 旱区地下水文与生态效应教育部重点实验室， 陕西 西安 710054)

1 研究背景

平原水库是在平原区天然湖泊、洼地、河道等地形的基础上，通过圈筑围堤而建成的水库[1]。在我国西北干旱及半干旱地区，由于受到人类活动和经济发展等多方面的影响，区域性的水环境问题、生态环境问题日益突出[2]，严重影响了地区社会发展及“一带一路”建设，而平原水库在解决这些问题和促进区域经济发展等方面起着重要作用。但平原水库的运行也往往会对区域地下水、地表水、土壤环境等带来一定的负面影响，产生库水富营养化、泥沙淤积、浸没、土壤盐渍化等一系列问题[3]，这对水库的运行效率和城市供水安全造成了严重的隐患。

针对平原水库运行引起的环境问题，国内外学者展开了大量研究。李荣荣等[4]采用遥感影像法分析了水库下游土壤盐分指数变化，结果表明水库蓄水使得下游土壤盐渍化程度不断加重；Winter[5]在充分考虑介质各向异性、地下水流系统、渗透系数等因素的情况下，运用垂向二维稳定流模型研究了水库和地下水之间的相互作用；韩菲等[6]系统分析了国内外湖泊及水库富营养化模型的优缺点，结果表明富营养化模型应当着重考虑水动力、化学、生物3方面因素的影响。

目前关于平原水库的研究多集中在库水水质评价及周边土壤生态环境等方面，对水库渗漏及周边地下水流场演变等方面的研究较少，而这正是威胁平原水库安全运行的常见问题[7-9]。在天然条件下，地下水系统会保持一定的均衡状态，而当水库蓄水之后，库水的渗漏补给作用改变了区域地下水系统的收支特征，打破了原有的地下水均衡动态[10]；另外，库底沉积物的分布、渗透性特征、地下水位变化会对水库的渗漏速率产生影响[11]。因此在充分考虑水库沉积物渗透性、地下水位特征的基础上，分析计算水库渗漏量，探究水库蓄水后区域地下水流场的变化就显得尤为重要。

本文以陕西省西咸新区斗门水库试验段为例，采用数值法和解析法相结合，分析库区渗漏量及水库蓄水对区域地下水流场的影响，以期为库底防渗、水库调蓄、地下水环境保护及斗门水库的进一步建设提供相关科学依据。

2 研究区概况

斗门水库在陕西省原汉代昆明池遗址的基础上修建而成，位于西咸新区沣东新城境内(原西安市长安区)，北接鱼斗路，西靠新韦斗路，距离沣河约2.6 km。斗门水库是引汉济渭工程的重要一环，在输配水工程中起着重要作用，是一座集调蓄、防洪、改善生态环境等多项功能于一体的水资源综合利用平原水库[12]。斗门水库试验段是斗门水库工程建设的一部分，于2017年3月正式蓄水运行，库区面积0.53 km2，其中水面面积约0.47 km2，库底高程396.8 m，正常蓄水位400.47 m，蓄水深度3.67 m，蓄水量155×104m3，地理位置如图1所示。

图1 研究区地理位置图

2.1 地形地貌

斗门水库试验段所处地貌单元为秦岭北麓近山前冲洪积平原，库区内渭河一级阶地大面积分布，阶面高程为396～410 m，整体上东南高、西北低，地形开阔，高程起伏较为平缓。

2.2 地层岩性

库区地表50 m以内地层主要形成于第四系，可分为以下几种类型：(1)上更新统(Qp3)风积层，其岩性主要为风积黄土(Q32eol)，层厚约5～12 m，外表呈淡灰黄色，垂直节理及大孔隙发育，质地较为疏松，分布于库区西南部。(2)全新统(Qh)冲积层，分为上下两部。下部冲积层(Q41al)分上下两层：上层厚约10～40 m，主要为壤土、粉土，夹少量中砂，于库区南部、西部(渭河一、二级阶地)广泛分布；下层岩性主要为粉质黏土，局部夹薄砂层，夹层厚度0.3～0.5 m。而上部冲积层(Q42al)岩性主要为壤土，层厚约5～10 m，受到人类活动影响，顶部人工堆积零星分布(约0.5 m，主要为水泥块、砖块等建筑垃圾)，于库区东北部(渭河一级阶地)分布。

2.3 水文地质条件

库区地下水类型为孔隙潜水，埋深在7.5～14.2 m之间，主要赋存于第四系黏性土层(夹部分少量砂层)中，大致由东南流向西北，含水层埋藏于山前洪积扇、阶地及漫滩。受到斗门水库建设影响，周边农田面积及人为灌溉有所减少，因此降水入渗补给和人工开采分别为斗门水库区域地下水的主要补给和排泄途径。

3 库区渗漏量计算

3.1 水库渗漏模型

在水库的建设过程中，库坝体内侧面做了较好的防渗处理，因此可以忽略内侧面的渗漏，仅考虑库底区域，根据库底沉积物岩性将水库划分为不同的渗漏分区进行分析计算。库水对地下水入渗补给的渗漏模型如图2所示，可概化为：(1)库水在水头差的作用下通过库底沉积物垂向入渗到含水层中；(2)在地下水流场的作用下，发生侧向渗流。根据达西定律[13-14]，采用公式(1)所示的离散化模型计算库区渗漏量：

(1)

式中：Q为库区总渗漏量，即库水补给地下水的总量,m3/d；Qi为分区渗漏量，即库水在各个分区对地下水的补给量,m3/d；Ai为各分区库底的水平面积,m2；ΔHi为库水位与所对应的地下水监测井水位之差,m；Li为监测井与所对应的分区形心的距离，m；Ki为各分区沉积物的垂向渗透系数，m/d；n为分区个数。

图2 水库渗漏模型示意图

3.2 水库渗漏分区

根据斗门水库试验段工程地质勘察资料，按照库底沉积物岩性，将库区划分为3个区域，渗漏分区及监测井位示意图如图3所示。其中Ⅰ-Ⅲ为不同的渗漏分区，GW1-GW6为沿库北岸自西向东布设的6口地下水位监测井。

图3 库区渗漏分区图

考虑到水库较深，不适宜进行野外原位试验测定沉积物渗透系数。根据土工试验标准[15]，在库区内共取样12组，采用室内变水头渗透试验法[16]测定了库底沉积物的渗透系数。变水头渗透试验法原理如图4所示：将野外所取的原状土样作为渗流通道，将变水头管的供水作为水源来模拟地下水的一维渗流过程；通过多次调节(5～6次)不同的水头高度进行重复试验；记录多次试验的时间和水头来计算分析试验结果。

图4 变水头渗透试验法示意图

根据试验数据，采用公式(2)计算库底沉积物的垂向渗透系数：

(2)

图4及公式(2)中：a为变水头管垂直于水流方向的横断面面积,cm2；L为变水头管两端过水断面之间的距离(渗径)，在本试验装置中等于土样高度,cm；KV为沉积物的垂向渗透系数,m/d；t1、t2分别为同一次试验中初始、最后读取变水头管水头的时间,s；H1、H2分别为前后两次读取的变水头管的水头,m。

通过野外取样及室内试验求得了各个分区库底沉积物的垂向渗透系数。各渗漏分区面积、库底沉积物岩性及垂向渗透系数如表1所示。

由表1可以看出：库底沉积物渗透系数量级为10-3m/d，其中库区东北部(Ⅱ区)沉积物渗透系数相对较大。

表1 渗漏分区相关参数

3.3 水库渗漏计算

结合表1中的相关参数，根据库区2018年3-7月的地下水位监测资料(由于库区施工限制，4月份监测资料缺失)，求得了库区渗漏量及日渗漏强度[17]，结果如表2所示。

表2 库区渗漏量及渗漏强度

由表2可以看出：库区日平均渗漏量为19.37 m3/d，由此可得年渗漏量约为7 070.05 m3/a，约占水库总库容(155×104m3)的4.6‰，库区平均日渗漏强度约为0.041 mm/d，年渗漏强度约为14.97 mm/a。从时间上看：3-7月库区渗漏量变化呈下降趋势，分析原因是由于3-7月随着降雨量的增大，库区地下水位有所抬升，而由于水库的运行调节，库水位基本无变化，因此使得区域水力坡度减小，库水对地下水的渗漏补给减小；从空间上看：库区东北部(II区)的渗漏强度较大，分析其原因是II区库底沉积物岩性以淤泥质土为主，且夹有部分中砂，因此造成了较大的渗漏损失，存在较大的渗漏风险。在水库运行管理中，应当着重加强库区东北侧的地下水位监测，做好防渗措施。

4 地下水流数值模拟

4.1 水文地质概念模型的建立

4.1.1 模拟范围 由于斗门水库试验段面积较小，大尺度的区域水流模拟不能精确反映斗门水库试验段蓄水对周边地下水流场的影响。因此采用等流量交换法[18]，概化模拟区范围。如图5所示：模拟区北至常七队，南到杨家庄村，距离约1.42 km，西接西安龙门中学，东至袁其寨一带，距离约1.85 km，总面积约1.83 km2。模拟计算的目标含水层为第四纪孔隙潜水含水层，厚度约55～60 m，含水层上边界为潜水面，下边界为隔水黏土层。

4.1.2 水力特征及边界条件概化 模拟区内含水层地下水主要表现为水平运动，仅在河道、库岸周边等区域存在垂向运动，整体地下水流速较小，符合层流运动的特征，适用于Darcy定律；渗透系数、给水度等水文地质参数在含水层不同点存在一定差异，为非均质性；地下水系统的补给、径流、排泄等要素随时间和空间有所变化，为非稳定流。综合分析模拟区的水文地质条件及地下水流的补、径、排条件，将模型边界条件概述如下：将黏性土隔水底板概化为潜水底部隔水边界；潜水面接受大气降水补给、库水入渗补给、蒸发排泄等垂向入渗补给，为不稳定水面；将东北及西南侧边界概化为零流量边界，东南侧及西北侧为二类流量边界。

4.2 地下水流系统数学模型

根据斗门水库试验段模拟区的水文地质概念模型，可建立如式(3)所示的数学模型来描述模拟区含水层地下水的二维水流运动：

(3)

式中：x，y为坐标变量；D为模拟区范围；h为水位,m；h0为初始水位,m；K为渗透系数,m/d；W为垂向补给强度,m/d；μ为给水度；Zb为含水层底板标高，m;q为二类边界上单位面积流入(流出)水量，m3/(d·m2);Γ1、Γ2分别为一、二类边界；n为二类边界外法线方向。

4.3 模型设计

4.3.1 模型网格剖分 根据实际水文地质条件和计算需求，采用不等距剖分的方法对所建立的模型进行了网格剖分，将模型剖分成了若干行和列，在库岸边进行了一定程度的网格加密。如图6所示：共剖分132行、160列，21 120个单元格，其中活动单元格15 789个，非活动单元格5 331个。

4.3.2 模型顶板和底板高程 通过DEM数据提取获取模拟区顶板高程数据，通过收集模拟区钻孔资料及工程地质勘察资料获取模拟区底板高程数据，采用克里金法进行插值，生成模拟区顶、底板高程如图7所示。

4.3.3 模型源汇项的处理 根据相关研究[19-20]，潜水埋深超过5 m，其蒸发可忽略不计，而模拟区潜水埋深在9～14 m之间，因此在模型源汇项的处理中不考虑潜水蒸发，主要考虑以下源汇项：(1)降雨入渗补给及蒸发。根据模拟区气象站实测降雨蒸发资料，模拟区降雨量多年平均值为605 mm，水面蒸发量多年平均值为592.9 mm。由于模拟区位于渭河一级阶地，根据西安市地貌图及模拟区相关地质资料[21-22]，确定模拟区降雨入渗补给系数为0.37，在模型中通过PRECHARGE模块进行设置。(2)对于地下水侧向径流补给，运用MODFLOW的Wells模块，在模型中设置抽注水井，通过侧向补给量确定井的流量。

4.3.4 模拟时长及初始流场 模拟期为2017年3月-2022年3月，在设置好模型各项参数后进行地下水稳定流场的模拟，得到模型的初始流场，如图8所示。

4.4 模型识别和校正

选择地下水流数值模拟中常用的“试错法”[19]对模型参数进行校正，即不断地对比观测井的模拟计算水位和实测水位，根据对比结果针对性地调整水文地质参数，直至观测井的模拟计算水位和实测水位误差在合理范围之内，则认为模型可用于区域流场的预测。在建立确定区域水文地质参数分区的过程中，参考斗门水库建设可行性研究结果，充分收集斗门水库试验段建设前期工程勘测资料[22-23]，划分出模型的水文地质参数分区。在此基础上对模型进行水文地质参数校正：选取2017年3月-2018年3月为模型的识别和校正期，根据2018年3月库岸边监测井的实测地下水位绘制校正对比流场，从而得出校正后的各分区及其对应的参数值，如图9、表3所示。

表3 模拟区各分区水文地质参数表

图5 模拟区范围示意图 图6 模拟区网格剖分示意图

图7 模拟区顶、底板高程(单位：m)

图8 模拟区初始流场(单位：m) 图9 模拟区水文地质参数分区图

4.5 结果与分析

利用模型识别和校正后的水文地质参数和边界条件，采用有限差分软件Visual MODflow 模拟了斗门水库蓄水5年后的地下水流场变化，模拟结果如图10所示。

图10 斗门水库蓄水5年后地下水流场(单位：m)

通过蓄水前后(图8及图10)的地下水流场模拟结果对比可以看出：蓄水前，受地形因素控制，模拟区地下水流向为东南向西北，整体地下水位介于387～396 m之间，库区地下水位介于391～394 m之间；蓄水5年后，模拟区地下水位抬升约0.1～0.5 m，整体地下水位介于388.5～396 m之间，其中库区东北部水位抬升较为明显，等水位线变得更为密集，水力梯度变大，地下水流向西发生偏转，结合表1、表2分析其原因是由于库区东北部渗漏强度较大，对原有的地下水流场产生了较大影响；库区南部地下水流场变化较小，地下水流向基本保持不变。综合以上分析可知，水库蓄水会对区域地下水流场产生一定影响，使得区域地下水位尤其是库区东北部地下水位产生了较为明显的抬升，地下水流向发生偏转，应当加强库区东北部的的地下水位动态监测，做好防渗措施。

5 结 论

(1)斗门水库试验段库底沉积物岩性以淤泥质土和壤土为主，东北部地层中夹有少量中砂；库底各分区渗透系数差异较小，量级均为10-3m/d，库区平均渗漏强度约为0.041 mm/d，即14.97 mm/a；库区东北部沉积物渗透系数及渗漏强度均较大，应当加强库区东北侧的地下水位动态监测，做好库底防渗措施。

(2)库区平均渗漏量约为19.37 m3/d，即7070.05 m3/a，年渗漏损失约占试验段总库容的4.6‰；库区渗漏量受降水量影响较为明显，丰水季节的渗漏量较小。

(3)水库蓄水对周边区域原有的地下水流场产生了一定影响，蓄水5年后使得模拟区地下水位整体抬升0.1～0.5 m左右，库区东北部水力梯度变大，地下水流向发生偏转。