1. 引言

地下水库是利用天然地下储水空间兴建的有拦蓄、调节和利用地下水流作用的特殊蓄水工程 [1] 。从20世纪二三十年代起，世界各国就开始了运用地下水含水层来调蓄淡水的工程实践，建立了为数众多的地下水库，并取得了显著的经济、社会和生态效益。瑞典、荷兰和德国的地下储水库分别可提供该国总供水量的20%、15%和10%以上 [2] 。美国从20世纪50年代就开始在地下咸水层内进行贮存淡水的试验，经过20几年理论和实践的发展，1973年后形成了含水层储存与回采(Aquifer Storage and Recovery, ASR)技术 [3] 。80年代以来，开始在佛罗里达州等干旱和半干旱地区推广这种利用含水层调蓄水资源的做法。到2002年，美国正在运行的ASR系统共有56个，而在建的系统则在100个以上 [4] [5] [6] [7] 。近年来我国兴建了八里沙河水库、黄水河水库、山东王河等多个地下水库 [8] [9] ，60年代，在上海市开展了地下水人工补给试验研究工作。70年代，采用地下水夏灌冬用和冬灌夏用以节省能源 [10] 。80年代，多省开展了地下水人工补给的试验工作 [11] [12] 。林学钰论述了地下水库的定义、设计论证及运行管理 [13] 。赵天石探讨了地下水库的概念、结构、勘察要点等 [14] 。李旺林等 [15] [16] [17] 对反滤回灌井进行室内试验研究和工程应用并推导出单井回灌量的计算公式。康华等 [18] 对秦岭山前的地下水库进行调蓄功能分析并建议在汛前加大开采强度以腾空库容。刘记来等 [19] 对北京西郊地下水库进行数值模拟和人工补给能力研究。魏晓燕等 [20] 对肥城盆地岩溶地下水进行数值模拟，预测出肥城盆地未来六年的地下水流场等。

近些年来，随着电子计算技术的快速发展，地下水数值模拟研究已经发展到空间三维流运动，也从单一的含水层发展到多个含水层，从稳定流发展到非稳定流，从饱和态发展到非饱和态，目前的地下水数值模拟技术的应用领域大大增加，且已经成为了当前研究地下水动态变化的主要手段与技术。目前常用的求解地下水流方程定解问题的方法有有限差分法和有限单元法。其中，美国地质调查局开发的地下水流模拟软件MODFLOW是基于有限差分法编制的一套用于孔隙介质中地下水流动数值模拟的软件。不需要对源程序进行任何修改，MODFLOW就可以直接用来解决大多数地下水模拟问题，已经被世界上许多官方和司法机构所认可。美国Brigham Young大学环境模拟研究实验室(Environmental modeling research-Laboratory)开发的GMS (Groundwater Modeling System)是基于概念模型的地下水模拟软件，具有良好的使用界面、强大的前处理、后处理功能及优良的三维可视效果，目前已成为国际上最受欢迎的地下水模拟软件，该软件中Modflow是地下水流模拟最常用的模块。本次工作即选用GMS中的MODFLOW模块作为建立地下水流模型的模拟软件。

乳山市位于山东半岛东南端，是山东省工农业发展水平较高的地区之一，用水需求不断增加，亟需开辟新的水源工程以增加供水量。本文分析了乳山市的水文地质结构，研究了地下水库建库条件，利用FEFLOW等大型地下水计算软件对乳山河地下水库进行地下水动态分析，建立水文地质概念模型，对地下水库进行数值模拟和调蓄分析。研究成果可为乳山市水资源开发和综合利用规划参考使用，为区域水资源合理配置、促进水资源开发利用与生态环境的和谐发展以及水利工程综合规划提供科学依据。

2. 乳山河地下水库建库条件分析

乳山河为乳山市第一大河，发源于马石山南麓的垛鱼顶，属半岛独流入海河流，于乳山口注入黄海，为季风性河流，河流全长65 km，流域面积1043 km²，其中威海乳山市境内流域面积为875.1 km²。乳山河流域属典型的暖温带海洋性季风气候，据降水资料统计分析，乳山河流域多年平均降水量为753.2 mm。乳山河流域呈阔叶型，在崔家河处河道缩窄，具备建设地下坝(地下截渗墙)的地形条件。

乳山市第二、第三水厂靠近乳山河附近，地下水库水源利用条件较好。乳山河上游龙角山水库、院里水库等水库，可利用其放水洞、溢洪道放水补源，地下水库有可靠的补给水源。乳山河流域多年平均年径流量为22,145万m³，多年平均径流深为271.3 mm，大部分径流流入大海，因此乳山河下游拦蓄利用地表水的潜力很大。

乳山流域中下游主要含水层为第四系孔隙潜水含水层，含水层岩性主要为砾质粗砂、中细砂、中粗砂，渗透系数为1.35 × 10⁻²~4.76 × 10⁻² cm/s，具有强透水性。砂层分布面积广，厚度为2.0~13.1 m，透水性和蓄水条件较好。含水层下伏花岗岩，且无穿越地表分水岭的断裂输水构造，基本不存在渗漏的问题，因此具备良好的地下水库建设条件。乳山河地下水库水文地质图详见图1。

3. 地下水库地下水流场数值模拟与调蓄分析

本次模拟区主要为乳山河地下水水库库区，包括乳山河自宅子夼村北至西耿家村河段，崔家河自官庄村东至西耿家村，司马庄河自楼村至凤台顶村东河段，总面积约23.8 km²。选用MODFLOW进行地下水流场数值模拟。

3.1. 水文地质模型

1) 含水层结构概化

乳山河地下水库模拟区含水层结构按照地下水埋藏条件确定为潜水含水层。综合考虑各地层的水力联系程度及含水层渗透性大小，以及岩性、特征，库区含水层岩性自上而下可概化为4层。

1、地表含砂壤土、壤土层 分布在乳山河地下水库库区地表，多位于地下水位以上，渗透性相近。

2、中上部细砂层 主要沿乳山河分布，多分布于两侧漫滩表层，部分位于地下水位以下。

3、中部淤泥质壤土层位于第四系覆盖层中部，渗透性较弱，富水性弱，分布不连续。

4、砾质粗砂层 大面积分布于库底，少量分布主河槽内，以砾质粗砂为主，渗透性强，富水性好。

2) 边界条件

通过对乳山河地下水库地下水系统的分析，认为乳山河地下水库为河道型地下水库，其储水地层均位于乳山河及支流河床上，是一个相对独立水文地质单元。其边界条件如下：

1、库区基底透水性较弱，隔水性能良好，基岩富水性差，为零通量边界。

2、库区地表接受大气降水、夜间凝结水，也存在蒸发排泄作用，当作一个源汇项赋值。

3、东部边界，南北山村东北基岩裸露处–邓家村–官庄村–西耿家村一线，视为相对隔水边界。

4、西部边界，宅子夼村东–车村–乳山寨–凤台顶村一线,视为相对隔水边界。

5、河道上游边界，宅子夼村北、崔家村南、圈港村南、楼村西，为主要补给边界。

6、河道下游边界，乳山河下游官庄村–乳山寨橡胶坝–凤台顶村北一线，是排泄边界，通过建设地下防渗墙，可形成不透水边界。

3.2. 水文地质参数

采用的渗透系数和给水度见表1。

3.3. 地下水主要补排项

模拟区含水层主要接受降雨入渗、乳山河河水侧渗和地下径流等补给，消耗于潜水蒸发、下游地下水径流排泄和人工开采。各均衡要素处理如下：

1) 补给量

1、大气降水入渗补给量

降雨入渗补给量按下式计算。

$Q_{降}=\alpha \times F\times P$ (1)

式中：α——降雨入渗系数；F——计算区面积(m²)；P——多年平均降雨量(m)。

其中，降水入渗系数壤土为0.1~0.2，细砂为0.27~0.32，砾质粗砂为0.35~0.45。

2、地下水侧向径流量

上游径流补给量按下式计算。

$Q_{径流补给}=K\times M\times I\times L\times T$ (2)

式中：K——计算断面的加权平均渗透系数(m/d)；M——计算断面的平均含水层厚度(m)；I——计算断面的平均水力坡度(‰)；L——计算断面长度(m)；T——补给时间(d)。

3、农业灌溉回渗量

由于目前农业生产技术条件有限，回渗量较小或者没有，故忽略不计。

4、乳山河侧渗补给量

乳山河侧渗补给量按下式计算。

$Q_{河流补给}=K\times M\times I\times L\times T$ (3)

5、沟谷洪流入渗补给量

沟谷洪流入渗补给量按下式计算。

$Q_{沟渗}=M_{地表}\times F\times \alpha$ (4)

式中：M_地表——地表径流模数(M = α径流系数 × P/T)；F——山区汇水面积(m²)。

2) 排泄量

1、农业区开采量，包括农业灌溉量和农村生产生活用水量

2、潜水蒸发量

潜水蒸发量按下式计算。

$Q_{蒸发}=C\times F\times \epsilon$ (5)

式中：Q_蒸发——浅层地下水蒸发量；C——潜水蒸发系数；ε₀——水面实际蒸发强度(m/a)，各气象站水面蒸发观测数据乘以折算系数0.62，转换成E601型蒸发皿水面蒸发值。

3、地下水下游径流排泄量

地下水下游径流排泄量按下式计算。

$Q_{径流排泄}=K\times M\times I\times L\times T$ (6)

4、水源地开采量

水库区周边分布多个村庄，部分村庄居民及牲畜用水通过开采地下水解决。水库区及边界处还分布了威海龙峰硅胶、富康生物肥料公司、乳山市华冶铸造有限公司等工业生产单位进行地下水的开采。

3.4. 地下水流数学模型

1) 地下水流数学模型的建立

根据模拟区地下水系统水文地质概念模型，建立了地下水流数学模型。当不考虑弱透水层贮水能力时，浅层潜水和承压水系统地下水运动的数学模型概化为三维不稳定流动系统，模型见式(7)。

$\{\begin{array}{l}\frac{\partial }{\partial x}\left(K_x\frac{\partial H}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_y\frac{\partial H}{\partial z}\right)+\frac{\partial }{\partial z}\left(K_z\frac{\partial H}{\partial z}\right)+\epsilon =S_s\frac{\partial H}{\partial t}\\ {H\left(x,y,z\right)|}_{t=0}=H_0\left(x,y,z\right)\\ {K\frac{\partial H}{\partial n}|}_{\left(x,y,z\right)\in \Gamma }=q\left(x,y,z\right)\\ H—地下水水头\left(\text{m}\right)\\ K_x,K_y,K_z—x,y,z方向渗透系数\left(\text{m}/\text{d}\right)\end{array}$ (7)

式中：K——边界法线方向渗透系数(m/d)；S_s——含水层比储水系数；Г——模拟区域第二类边界；H₀——含水层初始水头(m)；q——含水层二类边界单位面积过水断面补给流量(m²/d)；ε——源汇项强度(包括开采强度、入渗强度、蒸发强度等) (1/d)；Ω——渗流区域；n——渗流区边界的单位外法线方向，矢量。

2) 模型的设计与实现

1、网格剖分

采用矩形网格对渗流区进行离散化(剖分)。

2、参数初始化

本次数值模拟模型的模拟期为2014年9月到2015年9月，整个模拟期划分为12个应力期，每个应力期为一个相应的自然月，计算的时间步长为一天。每个应力期中，所有外部源汇项的强度保持不变。

边界条件：将所计算出来的各边界流入流出量输入到模型之中，适当调整边界流入流出量。

3、模型的识别与验证

由于受工作区地下水水文系列资料的限制，本次研究采用的方法为试估—校正法。根据己有的地下水动态观测资料，选取2014年9月1日至2016年9月1日为本次模型研究的识别阶段，两个月为一个时间段；选取2016年9月1日至2018年9月1日为本模型的验证阶段。选择2014年9月1日乳山河地下水库地下水位观测值为本次模型识别的初始流场，通过2014年10月1日和2016年9月1日的具有代表性的8眼观测井的的模拟地下水位和观测水位进行拟合，拟合结果显示模型在识别期和验证期内水位观测值与模型计算值拟合误差小于0.5 m的观测井数占到总井数的85%，而小于1 m的观测井数占总井数的93%。由此可以表明，该模型在识别和验证阶段均能够真实反映乳山河地下水流场，这也能够说明对研究区的边界条件及含水层结构等的概化基本符合实际情况。

3.5. 地下水库调蓄能力分析

乳山河地下水库调蓄方案：乳山河在库区内由西北向东南穿过，在冬春季节灌溉时上游来水量对地下水的补给以及雨季的弃水作为乳山河地下水库的主要补给源，区内的大气降水入渗补给及经过处理后的污水也是可利用的补给水源。在地下水流下游入海处建设地下截渗墙，减少库区地下水的天然排泄量，阻挡海咸水入侵。在乳山河不同位置建设橡胶坝，可以增大乳山河河道的蓄水量，增加地下水的入渗补给量，进而来增加储蓄量。

借助前面所建立的地下水流场模拟模型，利用美国地质调查局开发的地下水流模拟软件MODFLOW进行地下水库地下水流调蓄分析和计算，可分别求得各特征水文年的调蓄水量，计算结果见表2。

根据表2，乳山河地下水库建成后，有较大开发潜力。丰水年(25%)时水库的调蓄水量为1790万m³；平水年(50%)时水库的调蓄水量为1664万m³；枯水年(75%)时水库的调蓄水量为782万m³；特枯水年(95%)时水库的调蓄水量为404万m³。但从含水层水力条件及其富水性特征上分析，从调蓄空间上看，潜水含水层可调蓄空间适用性较大，但其水量在枯水年相对较小，合理的调蓄，优化储水用水方案是合理利用地下水资源的关键。

4. 结语

通过上述分析，得出以下结论：

1) 在乳山河流域中下游可以建设地下水库，地下水库实际总库容约2197万m³。

2) 从补给水源条件上看，拦蓄地表径流是主要补给水源，乳山河下游拦蓄潜力很大。

3) 乳山河地下水库建成后，有较大开发潜力。丰水年(25%)时水库的调蓄水量为1790万m³；平水年(50%)时水库的调蓄水量为1664万m³；枯水年(75%)时水库的调蓄水量为782万m³；特枯水年(95%)时水库的调蓄水量为404万m³。

4) 建设地下水库，不占用过多耕地、无大量移民搬迁问题且水量蒸发损失少，是提高乳山地区水资源利用率的有效手段，同时可有效防止海水入侵，具有良好的发展前景。

参考文献