作者：金芳义 郑秀清 陈军锋 马瑞群 宁显林
论文 关键词：地下水 数值模拟 modflow 水源地
论文摘要：本文首先建立了汾河二库及其下游地区水文地质概念模型，然后根据地下水水头长观资料识别模型的水文地质参数，运用可视化标准软件visual modflow，建立了地下水三维不稳定流模型来刻画汾河二库地下水下游区域流动的基本 规律 。运用识别后的模型预测了未来2010年和2015年水库不同蓄水位条件下的地下水动态变化。模拟结果表明，汾河二库的建成使周围地下水流动发生了变化，水库蓄水造成流域地下水位的上升，但对兰村水源地的影响有限。
汾河二库位于兰村泉域中西部太原市区西北约32km的汾河峡谷区，坝址下游8km为兰村水源地。兰村饮用水水源地是太原市规模最大、产水量最高的水厂。汾河二库建成以后，由于水库渗漏严重，对其周边地下水环境产生了一定的影响，研究汾河二库蓄水以后对兰村水源地地下水水位的影响具有一定的现实意义，可为当地政府合理、可持续开发利用、保护地下水资源提供决策依据，对汾河二库调蓄运行具有一定的指导意义。WWw.11665.cOM笔者将从研究区地下水的形成、循环规律入手，建立地下水三维流数值模型，进而利用该模型对区域地下水的水动力条件进行分析，预测分析水库在不同蓄水水位方案下对兰村水源地地下水动力场的影响。
1 研究区概况
研究区面积约162km2，多年平均降雨量为484.8mm，多年平均蒸发量为1806mm，属典型的北温带大陆性干旱、半干旱季风气候。研究区地势总趋势是北高南低，依地貌特征分成两部分：其北、西、西南部为侵蚀、溶蚀为主的石灰岩剥蚀构造地貌，出露地层主要为奥陶系和寒武系灰岩，海拔高度为1000m～2000m；东南部为盆地区洪积倾斜平原分布于西山区边缘，出露地层主要为第三、第四系松散堆积物，海拔高度为800m～830m。
研究区地下水补给、径流、排泄条件分成山区和盆地区。山区裂隙岩溶水主要接受大气降水和河流的入渗补给，地表径流渗漏主要是汾河渗漏，沿汾河罗家曲－镇城底、古交（寨上）－小塔为两个主要漏失段。山区岩溶裂隙水接受大气降水和河流的入渗后，沿东、西、北三个方向向盆地方向径流。从补给区至径流排泄区，地下水由无压水向承压水逐步过渡，排泄带岩溶水承压水头由于受地形地貌及地质条件的控制，承压水头从几米到几百米不等。天然状态下排泄主要是以泉的形式排泄，现状条件下排泄以人工开采为主。盆地孔隙水主要接受大气降水的入渗补给，汾河干流段地表水的渗漏补给，灌溉渠系及田间灌溉的入渗补给，山区岩溶裂隙水的侧向径流补给。由平原北部向南部，由东西两侧向中部径流运移。但由于平原地形相对较平坦，其流动缓慢。浅层孔隙水以潜水蒸发和向下游径流，现状条件下，由于长期大量超采，形成地下水位下降漏斗，地下水由四周向漏斗中心径流，以人工开采的形式集中排泄。
2数值模型的构建
2.1 水文地质概念模型
根据区内水文地质钻孔资料及地下水储存条件，以边山断裂带为分界线将研究区分为两个区。汾河二库所在区为奥陶系及部分寒武系岩层出露，属灰岩裸露区，由于灰岩厚度较大，故将其概化为一层。兰村水源地所在区地表为第四系覆盖物的覆盖，由于其透水性明显弱于其下面的奥陶系灰岩，故将该区在垂向上概化为两层，第一层为弱透水层，第二层为承压含水层。含水层均为奥陶系岩溶介质，概化为非均质各向异性介质，同一参数分区内含水层可视为均质，水流服从达西定律。根据区内的地下水动力场特征，研究区所有边界均为第二类边界。受 自然 条件及人为开采地下水的影响，研究区地下水流为非稳定流。根据实际长观系列资料，以2003年1月1日的水头分布作为研究区的初始水位。
2.2 数学模型
根据模拟区地下水系统水文地质概念模型，建立了 计算 区内水位分布数学模型。潜水和中深层承压水系统地下水运动的数学模型分别为:
潜水系统：
2.3 模型剖分
研究区的总面积为162km2，在平面上将其剖分为209行，184列，在垂向上分为两层，其网格为矩形，共计76912个单元格。其中，有效单元格45465个研究区的总面积为162km2，在平面上将其剖分为209行，184列，在垂向上分为两层，其网格为矩形，共计76912个单元格。其中，有效单元格45465个，无效单元格31447个，每个单元格的平面面积约为0.007km2。
2.4 均衡要素的确定
地下水系统的均衡要素是指其补给项和排泄项，均衡区为整个研究区。在均衡区内地下水补给项主要包括有大气降水入渗补给、地下水侧向径流补给、河道渗漏补给和二库渗漏补给；地下水排泄项主要包括有人工开采和径流排泄。均衡要素计算的目的在于确定地下水的各个补给排泄项随时间和空间的变化规律，为建立地下水流数值模拟模型准备数据。各均衡要素计算结果如表1。
表1 水文地质参数分区图
各均衡要素
大气降水入渗补给
河道渗漏补给
汾河二库渗漏补给量
人工开采量
侧向径流
补给量
（万m 3 /a）
1487.2
110
138
4666.8
排泄量
（万m 3 /a）
5107.6
1399.2
3 模型模拟结果及讨论
3.1 模型调较与验证
根据研究区地下水位观测资料的实际情况，选取2003年1月1日至2003年6月30日作为模型识别时段，历时半年，计算时段按应力期进行划分。在每个应力期内的开采量和补给量不随时间变化，分10个步长来计算研究区内的水位。
研究区水平方向分成四个参数分区：ⅰ区为西部和北部岩溶裸露地区，是本研究区的主要补给区；ⅱ区为沿河及河谷地区，渗透性较强；ⅲ区为薄层黄土覆盖区，受山前断裂带阻隔，水位雍高，是本区主要排泄区；ⅳ区为石炭二叠系出露，与奥陶系相比渗透性能稍弱。研究区水文地质参数分区见图1。
图1 水文地质参数分区图
利用试错法对模型参数进行率定，经过反复调参，得到较为理想的模型识别结果和合理的参数，各分区水文地质参数见表2。
表2 含水层水文地质参数分区表
地层序号
参数分区
渗透系数（m/d）
给水度
贮水率
x
y
z
上层含水层
ⅰ
42.4
42.4
8
0.253
ⅱ
62.3
62.3
8
0.225
ⅲ
16.3
16.3
8
0.137
ⅳ
36.3
36.3
8
0.137
下层含水层
ⅰ
42.4
42.4
8
0.000012
ⅱ
62.3
62.3
8
0.000013
ⅲ
46.6
46.6
8
0.000020
ⅳ
36.7
36.7
8
0.000013
根据所建立的数学模型，计算观测孔所在单元的水位，并和实际观测的水位进行对比，从而反求相关水文地质参数。由于研究区内的观测点较少，本次计算选用兰村水源地代表性水位观测孔s1进行水位拟合，水位拟合曲线见图2。由图可知，观测孔的实测水位和计算水位的变化趋势基本一致。
图2 s1观测孔模型识别阶段水位拟合曲线图
图3 模型验证s1观测孔水位拟合曲线图
本次识别各时段长观孔的水位 计算 值与实测值的拟合误差（绝对误差）小于0.1m的达到80%，符合识别要求。上述模拟结果表明，正负误差比较均匀，系统稳定性较好，模拟的精度和效果是理想的。
为进一步验证识别后的模型和水文地质参数的可靠性，用校正后的模型及参数组合计算出s1观测点的地下水位，将计算结果与实测水位相比较，对模型进行检验。根据计算区地下水位观测资料的实际情况，选择2003年7月1日至2003年12月31日期间的水位观测资料对模型进行验证。模型验证阶段的水位拟合结果见图3。
本次检验各时段水位拟合误差小于0.1m达85%，符合验证要求。结果表明通过水文地质条件的概化、边界条件的确定、水文地质参数的选取以及源汇项的处理后，所建的数学模型较好地反映了研究区的水文地质特征，该模型可用于地下水位动态预报。
3.2预测方案与结果分析
（1）模拟方案
为了能模拟二库蓄水以后对兰村水源地的影响，采取了两种模拟方案，即在现状年蓄水位和正常年蓄水位下，预测分析研究区岩溶地下水2010年和2015年的流场变化情况。
表3 模拟方案
不同方案
蓄水位
方案一
现状年水位(879m)
方案二
正常蓄水位（905.7m）
（2）结果分析
在现状年蓄水位下，各预测年地下水位降落漏斗以水源地为中心向四周扩展。2010年研究区内水位升幅主要变化在0m～30m之间，在兰村水源地地下水位降落漏斗一带水位降幅最大，可达30m，以770m等水位线所围成的面积为6.9km2；2015年与2010年相比，水位在水源地周围稍有减低，以770m等水位线所围成的面积为7.2km2。
在方案二下，与方案一相比，各预测年地下水位在二库库区明显上升，兰村水源地附近水位降幅减小。2010年研究区内水位降幅主要变化较方案一，二趋缓，在兰村水源地地下水位降落漏斗内以770m等水位线所围成的面积为6.8km2；2015年与2010年相比，水位在水源地周围稍有减低，以770m等水位线所围成的面积维持不变，二库库区地下水位明显上升。
图4 现状年水位下2010年等水位线预测图 图6 现状年水位下2015年等水位线预测图
图5 高水位下2010年等水位线预测图 图7 高水位下2015年等水位线预测图
参考 文献
[1] 冯玉明.泉枣沟水源地地下水保护区划分报告[r].太原，2009.
[2] 李俊亭.地下水流数值模拟[m].北京地质出版社，1989.
[3] 周念清，朱蓉，朱学愚.modflow在宿迁市地下水资源评价中的应用[j].水文地质工程地质，2000（6）:9-13.
[4] 刘建忠.太原市汾河蓄水工程对兰村泉域水资源的影响[j].太原科技，2003(4):47-49.
[5] 赵兰霞.左海凤.兰村泉域地下水位预测[j].研究人民黄河，2007(3):49-50.
[6] 王秀杰，练继建，杨弘.石家庄淖沱河地区地下水库研究[j].水利水电技术，2004 (9):5-7.
