汾河兰村—柴村段河流水位抬升对潜水流场影响的数值模拟研究

王瑞芬 季晨龙 李海明

(1.山西省水利水电科学研究院，山西 太原 030002；2.天津科技大学滨海地下水利用与保护研究室，天津 300457；3.天津市海洋资源与化学重点实验室，天津 300457)

由于经济社会的快速发展和人口的急剧增长，汾河流域生态环境受到严重破坏[1]。山西省先后组织了4次大规模治理，基本实现了汾河干流全年不断流，流域地下水位持续回升，但流域自身的产流功能仍未恢复[2]。汾河流域历史上拥有丰富的湿地资源，在涵养水源、降解污染、防洪蓄水等方面起到重要作用[3]，但在汾河早期的治理中，原本属于河道的大量滩涂、湿地被开发利用，湿地资源不断减少[4]。为此，山西省委省政府实施了汾河流域生态修复中游核心区干流蓄水工程，该工程主要通过河道蓄水向堤外湿地补水，蓄水期间使流域内湿地水位抬升，蓄水水位保持在滩面1～2m。

流域内湿地水位的抬升使得在湿地补给地下水的过程中面临各种技术与科学问题，例如湿地水位的升降容易引起土地盐碱化、沼泽化等问题。目前针对蓄水工程对下游地下水的影响多从渗流量[5]、环境流量[6]角度展开，也有研究者运用不同的方法，对汾河流域的地下水位进行了动态预测[7-10]，但对湿地水位抬升影响下的地下水流场特征研究较为罕见。

数值模拟也已成为研究河水与地下水交换的重要手段[11-14]。其中，对河流与含水层系统模型的正确概化以及对模型参数的准确估算是影响研究结果的关键因素[15-16]。本文运用数值模拟的方法，研究了汾河兰村—柴村段在河流水位瞬时抬升0.5m、1.0m、2.0m后的潜水流场特征，可为山西省汾河流域湿地水资源管理和修复提供相关参考。

1 研究区水文地质条件

研究区位于太原市西张盆地，具体位置见图1。西、北两面被群山环抱，海拔多在780m以上，面积约为66km2。兰村—柴村段的河流属平原性河流，长度约10km，宽约50～300m，随着河流流程的增加，含水层富水性增大。地下水主要接受汾河干流段地表水的渗漏补给以及大气降水入渗补给，地下水位低于河水位，河水通过入渗向河两岸含水层侧向渗透补给地下水。潜水的地下水径流主要受地形、含水层岩性、补给来源等控制。总体上由北部向南部，由东西两侧向中部流动运移。

图1 研究区位置

研究区地下水为松散岩类孔隙水，含水层多为全新统砂卵砾石层，在0～160m埋深范围内分为潜水含水层和弱透水层，包气带厚度约10～22m，岩性多为粉土、亚砂土。潜水含水层顶板高程为620～824m，厚度约30～150m，岩性多为砂卵石。弱透水层厚度约10～80m，底板埋深40～160m，岩性多为黏土。

潜水的排泄方式主要有蒸发、越流排泄及含水层的侧向径流3种形式。

2 潜流带水动力特征数值模拟

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

2.1.1 边界条件

根据地形地貌、地层结构及区域流场特征，潜水含水层南边界与东边界中下部概化为定流量边界，东边界上部、西边界以及北边界概化为隔水边界，弱透水层均概化为隔水边界，见图2。

图2 潜水含水层边界条件

2.1.2 含水层空间概化

本文只关注与河流水力联系密切的潜水含水层，根据已有的地质图、剖面图、地质报告等资料将含水层概化为非均质各向同性潜水含水层，地下水流以水平运动为主。研究区含水层三维结构与剖面概化见图3。

图3 研究区含水层三维结构与剖面概化

2.2 地下水流数值模型的建立

2.2.1 空间剖分

对模拟研究区进行网格剖分，共剖分为20000个网格(100行×100列×2层)，其中活动单元格为8368个，垂向上划分为两层，见图4。

图4 网格剖分示意图

2.2.2 模型识别应力期选取

根据现有资料，模型以2012年1—12月为校正期，2013年1—12月为验证期，以月为应力期，每个应力期为1个时间步长。

2.2.3 源汇项处理

本文模型中潜水含水层地下水补给项主要包括降水入渗补给、地表水渗漏补给，排泄项主要包括蒸发、径流排泄。

降水入渗补给量根据降水量乘以降水入渗补给系数得到；地表水渗漏补给根据水文站径流量测量资料，按一定渗漏比计算得到；潜水蒸发根据已有资料，潜水含水层现状地下水位埋深在临界埋深(5m)以下的区域不考虑潜水蒸发，在临界埋深以上的冲积、洪积交界洼地区潜水蒸发系数值采用0.040，冲积平原区采用0.013；边界径流排泄量采用达西定律计算得到。

2.3 模型的校正与验证

根据研究区校正期和验证期潜水流场模拟值与实测值的比较，两者的拟合程度较高，说明所建立的水文地质概念模型和数值模型基本合理，该模型能较好地模拟研究区地下水动态变化过程。

2.4 情景假设

本文以2014年1月地下水位为初始条件，预测了河流水位瞬时抬升0.5m、1.0m、2.0m时引起潜水含水层地下水位的变化情况，并且在每个河流抬升幅度下选取5日、10日、30日、60日、90日来举例说明。

预测期模型主要考虑的补给项有降水入渗补给、地表水渗漏补给，排泄项主要有蒸发、径流排泄。其中降雨量采用2004—2014年降雨量的平均值作为降水入渗补给输入模型，其他参数与2013年保持不变。

3 结果与讨论

研究区初始潜水流向为北部由北向南；中部受河流影响，流向变为由西北向东南；南部为由东西两侧向中间流动。北部等水位线密集，潜水流速快，水力坡度约为0.0065，流速约为0.220m/d。中部与南部等水位线稀疏，流速与北部相比较缓慢，水力坡度约为0.0009，流速约为0.041m/d。研究区内河流水位高于潜水水位，河流与潜水的补给关系为河流补给潜水。河流水位瞬时抬升0.5m、1.0m、2.0m后，在5日、10日时对潜水流场整体影响不大，但河流处的等水位线出现微小变化，河水补给两侧潜水，河流对潜水的补给作用增强，见图5、图6。

图5 5日潜水模拟流场

图6 10日潜水模拟流场

随着时间的推移，河流对潜水的侧渗补给不断加强，河流水位瞬时抬升0.5m、1.0m、2.0m时，河流附近的潜水流向逐渐变为从河流向东西两岸流动，特别是当河流水位抬升1.0m、2.0m时变化最为明显，见图7、图8。

图7 30日潜水模拟流场

图8 60日潜水模拟流场

河流水位瞬时抬升0.5m、1.0m、2.0m之后，经过90日，研究区潜水流向在由北向南的趋势上，由于河流的侧渗作用，河流附近的潜水流向明显变为由河流向东西两岸补给，见图9。北部水力坡度与潜水流速变幅较小，中部、南部水力坡度和流速与初始流场相比均有降低，水力坡度为0.0004～0.0006，流速为0.017～0.027m/d。

图9 90日潜水模拟流场

图10 河流水位抬升0.5m影响距离

图11 河流水位抬升1.0m影响距离

将相同时刻下河流水位抬升影响后的潜水流场与初始流场相比较，绘制出河流水位抬升不同幅度时对潜水的影响距离，见图10～图12。由图10～图12可以看出，河流水位抬升相同幅度下，随着时间的推移，影响距离逐渐增加，并且由于河流右岸空间有限，在河流水位抬升0.5m 60日、河流水位抬升1.0m以及2.0m 30日时影响距离已经超出研究区边界。河流水位瞬时抬升后，左岸90日时部分区域的影响距离也已超出研究区边界。

图12 河流水位抬升2.0m影响距离

根据图10～图12，得出河流抬升不同幅度时，经过不同时间，其相应的影响距离，见表1。多数情况下，河流水位抬升幅度越高，相同时刻下的潜水影响距离越大。

表1 不同时间不同河水位抬升幅度潜水水位影响距离 单位：m

4 结 论

a.研究区初始潜水流向为北部由北向南；中部受河流影响，流向变为由西北向东南；南部为由东西两侧向中间流动。北部等水位线密集，潜水流速快，水力坡度约为0.0065，流速约为0.220m/d。中部与南部等水位线稀疏，流速与北部相比较缓慢，水力坡度约为0.0009，流速约为0.041m/d。河流水位高于潜水水位，河流与潜水的补给关系为河流补给潜水。

b.随着时间的推移以及河流水位抬升幅度的增加，河流对潜水的侧渗补给不断加强，使河流附近的潜水流向最先开始改变，经过90日之后，潜水流向在由北向南的趋势上，河流附近的潜水流向已明显变为由河流向东西两岸补给。

c.多数情况，河流水位抬升相同幅度时，随着时间的推移，影响距离逐渐增加；河流水位抬升幅度越高，相同时刻下的潜水影响距离越大。河流水位瞬时抬升0.5m、1.0m、2.0m时，河流左岸60日影响距离分别为1742～2345m、1541～2613m、1809～3015m；90日时部分区域已超出研究区范围。