基于Modflow和LAK3的三泉水库渗漏量研究

王亚琴 杨军耀 刘庆 任柳妹

摘要：水库渗漏补给是泉域复流的重要形式，确定水库渗漏量及影响因素对于预测泉域复流很有必要。以三泉水库为例，运用地下水模拟软件CMS中的Modflow和LAK3程序包，建立水库和地下水含水层耦合模型，模拟预测不同关井压采方案对水库渗漏量及古堆泉口水位的影响程度，并结合正交试验选取渗漏区域的垂向渗透系数、平均渗透性能、给水度和贮水系数进行参数敏感性分析。结果表明：实施关井压采后，水库渗漏量增加，泉口水位恢复加快，对渗漏量的影响大小顺序为平均渗透性能>给水度>垂向渗透系数>贮水系数。

关键词：水库渗漏量;GMS;参数敏感性;三泉水库;古堆泉

中图分类号：TV697.2; P641.2

文献标志码：A

doi：10.3969/j.issn.1000- 1379.2019.02.013

1 研究区概况

古堆泉域地处山西省南部，地理坐标为东经111°5′-112°5′，北纬35°28′-36°0 ′。古堆泉出露于整个泉域的西侧，位于运城市新绛县古堆村，主要含水层为寒武、奥陶系灰岩及其上覆第四系松散沉积物层。受区域构造影响，古堆泉成为山西省重要的岩溶大泉之一。随着社会经济的发展，受井群开采和矿坑排水等影响，古堆泉水出现干枯断流。

三泉水库也位于泉域的西侧，坝址临近古堆泉。该水库在试通水期间，水位下降，水库渗漏量较大，库尾渗漏区出现了18处大小不等的渗水塌陷坑，塌陷坑和主要泉眼高程为442.0～ 443.6 m。经渗漏分析，在覆盖层较厚的区域，库水沿着原有的泉水出溢通道回灌补给岩溶水。主渗漏区位于水库库尾，即古堆泉口，水库渗漏量的大小直接影响着泉域的岩溶地下水量及古堆泉口水位变化。

水库渗漏作为泉域复流的一种方式，研究水库渗漏量对于预测泉域复流具有重要的现实意义。水库渗漏量的计算方法有很多，较为成熟方便的是利用软件模拟计算渗漏量[1-4]。笔者运用地下水模拟软件GMS中的Modflow和IA K3程序包，建立库水和地下水耦合模型，预测水库在不同运行水位并施行关井压采情况下的渗漏量变化及对泉域岩溶水位的影响，通过监测泉口水位计算古堆泉复流年限，并结合正交试验来分析水库渗漏量的主要影响因素。

2 地下水流数值模拟

目前应用广泛的地下水数值模拟软件GMS不仅拥有良好的可视化操作界面，而且运算效率高，处理功能强大[5]。本研究选用GMS中的Modflow和LAK3程序包进行三泉水库渗漏过程模拟。首先根据泉域的水文地质条件建立水文地质概念模型，然后运用Modflow程序建立地下水流模型，将水库条件通过LAK3模块添加进去，最后将库水与地下水耦合起来计算。

2.1 水文地质概念模型

（1）计算区域及边界概化。模拟区地下水由北东、东南方向向古堆泉方向流动，为保证东部边界水位变化不受三泉水库渗漏的影响，设定475 m岩溶地下水等水位线作为东部边界，概化为定水头边界：北部以汾阳岭北缘断裂为边界，该断层以北岩溶含水层埋深在800 m以上，岩溶地下水处于循环缓慢的滞留性状态，故将该边界概化为滞留性隔水边界：西部九原山一汾阳岭西缘断裂概化为隔水边界：南部侯马盆地南缘断裂概化为隔水边界;模拟区计算面积约1 026 km。

（2）含水层概化。第一层为上部广泛分布的第四系孔隙水含水层和局部分布的岩溶水、风化裂隙水含水层，底板埋深30～350 m;第二层为第三系弱透水层和岩浆岩极弱透水层，底板埋深为175—495 m;第三层为下部的奥陶系岩溶水含水层。

（3）源汇项。研究区内的源汇项主要由5部分组成：①降水人渗补给;②灌溉回渗补给（主要为引黄灌溉水）;③三泉水库渗漏补给（渗漏主要发生于水库库尾）;④蒸发（只考虑汾浍河谷范围蒸发）;⑤人工开采（模拟区内的农业灌溉水井和农村生活用水井均概化为面井，开采量采用估值并在模型识别时调整）。

2.2 数值模型

（1）数学模型建立。研究区内含水层水文地质参数不均一，位于泉口附近的水库渗漏明显，故按非均质、各向异性三维非稳定流进行数值模拟，数学模型为

（2）水库（LAK3）模块。水库通过LAK3模块概化，利用Triangle Interplation方法进行插值生成水库TIN，从而确定水库的形状和位置。需要设置的参数主要有初始水位、水库平均渗透性能、降水量、径流量、供水量，其中：水库平均渗透性能K’由水库库底淤积物渗透系数和淤积物厚度比值决定，水库水位根据水库实际运行方式赋值。

（3）参数分区。根据水文地质特征、地形地貌单元及地层岩性将模拟区划分为汾浍河谷阶地区（I）、碳酸盐岩深埋区（Ⅱ）、碳酸盐岩浅埋区（Ⅲ）、碳酸盐岩裸露区（Ⅳ）、塔儿山一九原山中低山区cv）s个参数分区，其中：中低山区根据岩性差异性又可分为碳酸盐岩一岩浆岩浅覆盖区（ vi）和碳酸盐岩裸露区（v2）两个亚区（见图1）。

（4）模型求解。根据模拟区水文地质条件、三泉水库位置和模拟区岩溶含水层埋深状况，将模型剖分为109x105x3（列×行×层）共18 138个计算单元（活动单元）。模拟期为2016年6月-2056年1月，模拟总时长为39.5 a。

（5）模型校正和验证。选取2016年6-9月作为校正期，对模型中的参数进行校正。选取2016年10月-2017年1月作为模型验证期，运用该阶段观测孔资料来验证模型。经过计算水位和实测水位的多次拟合，观测孔的拟合误差均小于10%，符合《地下水資源管理模型工作要求》[6]对模型误差的要求，说明该模型参数取值合理。校正末期岩溶水位等值线见图2。

3 渗漏量与水位预测

3.1 方案预设

为研究水库不同蓄水位状态，结合不同关井压采方案的水库渗漏量以及古堆泉口岩溶水位变化特征，设定4种方案：①现状开采条件下三泉水库低水位运行，运行水位为443.0 m;②现状开采条件下水库正常调节库水位运行，运行水位为447.5 m，③水库正常蓄水运行下实施一期关井压采方案：④水库正常蓄水运行下同步实施一二期关井压采方案。其中：一期关井压采方案为2017-2018年压采岩溶水342.9万m：二期关井压采方案为2020-2021年压采岩溶水706.9万m。

3.2 预测分析

（1）模拟期内，水库渗漏影响最为显著的古堆泉口水位变化在4种方案下均经历3个阶段，即短期快速抬升一多起伏波动抬升一等幅抬升（见图3）。古堆泉口水位由379.7 m恢复至443.0 m（即复流），4种方案下分别需要39、39、37、35 a。

（2）水库低水位和正常蓄水位条件下，末期渗漏量相差不大，均约为1.105 0万m/d，初期的水库渗漏量在正常蓄水位下为2.183 8万m/d，较低水位下大0.368 2万m/d。水库正常蓄水位条件下实施一期关井压采和同步实施一二期关井压采方案后的末期渗漏量相差不大，较正常蓄水位条件下分别大0.002 0万、0.002 1万m/d。4种方案下渗漏量变化趋势相似，以水库低水位运行条件下为例，水库渗漏量变化和泉口水位变化历时曲线见图4，随泉口水位上升水库渗漏量迅速变小，随后波动上升并逐渐趋于稳定。

（3）水库蓄水之后，库水通过泉口渗漏通道人渗，渗漏速度快，渗漏量大;地下水漏斗具有中间小、四周大的形态特征，古堆泉水位在开始时恢复速度相对较快;随着水库淤积加深，淤积厚度增大，渗漏速度变慢，渗漏量减少：库水渗漏补给地下水漏斗，随着泉口水位抬高，相应的水头压力变大，漏斗面积增大，水位回升速度相对变慢：随着水位持续回升，在水库渗漏补给古堆泉的过程中，含水层结构发生变化，包气带变为饱水带，水库相当于顶托式地往漏斗内灌水，渗漏量相对增大，渗漏速度也变快;最后水位恢复到泉口时，库水与岩溶水达到相同水头，渗漏量和水位回升达到平衡。

4 参数敏感性分析

4.1 介质参数及灵敏度指标选取

在地下水流模型建立过程中，选取与岩溶水流场相关性较强的因素作为敏感性分析因子，模拟在水库正常蓄水条件下同步实施一二期关井压采方案的水库渗漏情况。选定水库所在分区的垂向渗透系数、水库平均渗透性能、给水度、贮水系数4个因素作为敏感性分析因子，同时选出5组水平值得到四因子五水平的正交设计表L25（5）[7]。在设计正交试验表时，将以上选取因素作同等幅度的变化，以便在相同条件下比对参数的敏感度。选取模型最终识别值为标准值，其余水平值在该值基础上分别上下调整20%、40%。选取水库的初期渗漏量作为评价指标。各参数水平值见表1。

4.2 试验结果分析

将这25组试验值代人模型中进行GMS计算，利用在地下水数值模拟研究中创建专用图形的GW_Chart程序的运算器和水文提取器功能输出LakePolts.从而得到渗漏量随时间的变化值。正交试验结果及初期渗漏量见表2。

最常用的正交试验数理统计分析方法是极差分析法，即通过极差的大小来反映评价指标对各个因素的敏感程度，极差越大说明指标对该因素的敏感性越高，反之则敏感性越低[8]。对正交试验结果采用极差法进行分析，结果见表3。影响水库渗漏量的因素极差大小顺序为平均渗透性能>给水度>垂向渗透系数>贮水系数。

通过计算4个因子在不同水平下的极差可知，水库初期渗漏量对平均渗透性能最为敏感，其次为给水度，而对渗漏区域的垂向渗透系数和贮水系数敏感性较低，因此在数值模拟过程中对水库的平均渗透性能准确性要求最高。

为研究水库平均渗透性能K’与水库初期渗漏量Q的关系，在其他条件不变的情况下只改变K’的大小，得到K’与水库初期渗漏量Q的关系曲线（见图5）。可知，水库平均渗透性能K’与初期渗漏量Q基本成线性正相关关系，表达式近似为Q=1 030.7K'+ 0.086 7，即水库平均渗透性能K’越大，初期渗漏量Q也越大。

5 结论

（1）数值模拟结果表明，水库低水位和正常蓄水位条件下，水库末期渗漏量相差不大，均约为1.105 0万m/d;实施关井压采后，方案③和方案④末期渗漏量较正常蓄水位条件下分别大0.002 0万、0.0021万m/d。

（2）在方案③和方案④下，泉口水位恢复至443 m（即复流）分别需要37、35 a，较未实施关井压采方案分别提前2、4a，说明关井压采力度加大后可提高泉口水位的恢复速度。

（3）参数敏感性分析表明，影响水库渗漏量各因素极差大小顺序为平均渗透性能>给水度>垂向渗透系数>贮水系数，其中：水库的平均渗透性能K’值与初期渗漏量Q成线性正相关关系，表达式近似为Q=1030.7K’+0.086 7。因此在地下水流随机模拟过程中，应对水库平均渗透性能的准确性提出更高要求。

参考文献：

[1]程春龙，束龙仓，鲁程鹏，等，应用MODFLOW和LAK3计算水库渗漏量[J].工程勘察，2013，41（7）：31-34.

[2] 吕路，杨军耀，秦嘉楠，汾河三期蓄水工程人工湖渗漏对太原市浅层地下水影响的数值模拟[J].水电能源科学，2017，35（1）：89 - 92.

[3]介飞龙，朱志新，封丽华，等.LAKE在水库地下水影响预测中的应用[J].中国水运，2016，16（1）：121-123.

[4]刘涵宇，夏强，张世殊，等，基于水文地質参数敏感度分析的水库渗漏量不确定性评价[J].科学技术与工程，2017，17（14）：32-38.

[5]祝晓彬，地下水模拟系统（GMS）软件[J].水文地质工程地质，2003（5）：53-55.

[6] 王兆馨，林学钰，余国光，等，地下水资源管理模型工作要求：GB/T 14497-1993[S].北京：中国标准出版社，1994：8.

[7]刘瑞江，张业旺，闻崇炜，等，正交试验设计和分析方法研究[J].实验技术与管理，2010，27（9）：52-55.

[8] 郝静，贾仰文，张永祥，等，应用正交试验法分析地下水流模型参数灵敏度[J].人民黄河，2015，37（9）：66-68.