MODFLOW在库区浸没预测上的应用

2014-03-07周子东河南省水利勘测有限公司

河南水利与南水北调 2014年8期

□陈鹏 □孙刚 □周子东（河南省水利勘测有限公司）

0 引言

预测水库蓄水后库周水位分布情况是研究浸没发生范围、程度的关键一步，在查清地块水文工程地质条件的前提下运用MODFLOW可以帮助我们较准确的预测蓄水后水位的变化，从而对比临界埋深得出可靠的结论。

1 MODFLOW简介

MODFLOW是地下水模拟系统 (GroundwaterModeling Systems)的一个计算模块，是美国地质调查局在20世纪80年代开发出来的专门用于孔隙介质中三维有限差分地下水流数值模拟的软件，问世以来，已经在全世界范围内科研、生产、环境保护和水资源利用等领域得到了广泛应用，成为最为普及的地下水运动数值模拟的计算软件。应用MODFLOW进行地下水流模拟的步骤如下：

渗流区的离散化→建立地下水流动问题的差分方程组→求解差分方程组。

2 研究地块的工程水文地质条件

拟建水库位于淮河干流上游。因水库左岸地形平缓开阔，第四纪冲洪积物分布广泛，地下水位埋深较浅，初判将遭受较大范围浸没灾害的影响。选择了龙井～下土城地块为作为浸没勘察研究的典型地块。本区气候属华中副热带气候区。多年平均降雨量为1059mm。多年平均水面蒸发量790mm，陆面蒸发量650mm。

该地块地形总体由东北向西南河谷缓倾，地块以北为由中更新世冲洪积物组成的岗地，向南依次降为二级阶地、一级阶地及河床漫滩，地表坡降多在2.00×10-3～2.00×10-2之间。

据现场调查及钻探资料可知，本区分布的地层中有第三系和第四系，由老到新简述如下：一是第三系毛家坡组（Em）以紫红色粘土岩为主，夹有泥质填充的砂砾岩；二是中更新统（alplQ2）为褐红色低液限粘土；三是上更新统（alplQ3）上部为低液限粘土，下部为砂砾石；四是全新统（alQ4）上部为低液限粘土，下部为细砂土，构成一级阶地的主体。水库水边线位于一级阶地。

表1中的水文地质参数是通过现场和室内试验，结合类似工程和经验确定。底部第三系粘土岩和砂砾岩，透水性小，可作为本区相对隔水层。场区等水位线近似平行于河床，其补给来源以后缘中更新统组成的岗地以潜流形式补给，其次为降水补给。地下水的排泄方式有两种：一是以潜流和泉的形式向游河排泄；二是通过民井排泄。

表1 龙井～下土城地块水文地质参数建议值表

龙井～下土城地块可能浸没范围内主要为一级阶地Q4低液限粘土，通过野外试验确定土壤毛管水上升带的高度取62.50 cm。

3 地块的地下水流动数学模型

研究区的地下水主要储存于砂砾层中，存在统一的地下水面，地下水流动系统遵从水均衡原理和达西定律。研究区域地下水流从空间上看以水平运动为主，含水介质的岩性特征在垂向上变化较显著，为保证计算准确性将水库蓄水前后的地下水运动均概化为非均质各向同性三维稳定流。

根据上述的概念模型，结合研究区的边界条件，可得出研究区地下水流动数学模型：

式中：h0—初始水头；H—含水层厚度；n—含水层边界上内法向单位矢量；Ω—模拟区范围；B1—模拟区上下游定流量边界；B2—模拟区两侧零流量边界；q—含水层边界单宽流量，流入为正，流出为负。

4 地块的空间离散化与边界条件的确定

龙井～下土城地块水文地质结构主要为非均质多层水文地质结构。模拟区域长1980m，宽1200m的地块，模拟区面积为1.90 km2，见图1。模型两侧边界分别取在2个勘探线外侧，并与地下水等水位线垂直，为零流量边界；下游边界取已知水头边界，其值为各处的河水位；上游边界取在Q2和Q3地层界线附近，为定流量边界，其值可根据达西定律确定：

依据实测的地下水等水位线，分别计算出弱透水层和强透水层流量，最后叠加得到模拟区上游边界流量约2000m3/d。

利用GMS(Groundwater Modeling System)中Modflow模块建立栅格模型。按模拟区地层岩性将模拟区分成4层，每层33×20个单元格，其中每层有效单元格为577个，总共有效单元格为2308个,见图2。分别给出4层的顶底板高层，然后赋给各单元格水文地质参数，并进行地下水流动模拟。实测水位时间为3月份，因降雨较少，所以模型识别时不考虑降雨。

图1 龙井～下土城地块模拟区边界条件示意图

5 典型地块数值模型的识别与校验

由对水文地质结构的认识到水文地质概念模型再到地下水数值模型的整个过程经历了一系列概化处理，这些因素导致模型不可能细致地、完全真实地刻画出各处的地下水运动，而只能对模拟区内地下水运动进行宏观的、趋势性的拟合。为了判别所概化的模型与实际地质原型的符合程度，必须对模型进行识别与检验。

采用试算法进行校验，不断调整含水层参数，使模拟得到的水位与实测水位吻合。模拟的水位与实测水位有一定的差别，这是因为：①概化的水文地质概念模型是理想化的表述，例如岩层中存在的透镜体难以一一刻画；②边界条件都是人为边界，与实际的边界条件也有一些差别；实测的等水位线是根据离散的水位值经插值而得，其精度有一定局限。

尽管如此，模拟所得的水位和观测水位形状相似，趋势相近，两者误差最大不超过0.50m，地下水流动方向和观测的流动方向基本一致。模型可以用于预测典型地块水库蓄水至90m时地块的等水位线变化情况。计算参数采用经模型识别调整后的计算值，如表2所示。

表2 龙井～下土城地块浸没预测计算采用水文地质参数表

6 水库蓄水后地块边界条件的确定

水库蓄水后，在龙井～上土城形成三面临水的半岛，其地下水流动方向将发生变化。模拟区下游边界为库水位定水头边界，上游边界偏安全考虑，仍然为定流量边界，但其水位将随库水位上升有所壅高。取平行于流线的一个剖面，下游边界水位初始高程h为85m，蓄水后至90m；上游边界初始水位h1为90.50m，距离库水边界约900m，可以近似计算得到库水位上升后上游边界的水位。当库水位上升至90m时，上游边界水位上升至90.89m，偏安全考虑取91m。蓄水后龙井～下土城地块模拟区面积为1.60 km2。

降雨入渗作为模拟区上部边界条件，在模拟初始水位时，由于实测水位的时间为3月下旬至4月上旬，降雨量小，所以在模型识别与校验时未考虑降雨影响。在预测水库蓄水后，典型地块地下水位壅高时，应当考虑降雨入渗的影响。

模型计算时，偏安全考虑取6-8月降雨集中时间，平均日降雨量为2.50mm/d，两典型地块地表都以低液限粘土为主，其入渗系数取经验值0.01，得到降雨入渗量为2.50×10-5m/d。

7 水库蓄水后典型地块地下水位预测

如前所调整的模型边界及计算参数，代入模型计算，可以得到地块模拟区在水库蓄水以后的等水位线图，进而作出模拟区的地下水埋深等值线图。在龙井～下土城地块模拟区，在蓄水后模拟区的渗流场特征比蓄水前发生了明显变化。其一是随着库水位上移至高程90m，从宏观上看，全区地下水位上升后的运动总趋势大致仍由东北向西南流动。但在上土城和下土城将变成三面临水的半岛状库岸，地下水呈扇状流动，分别向淮河和半岛两侧的冲沟或低地（都是库区）排泄。其二是随着库水位上升，地下水位发生不同程度的壅高，升幅多在0.40～4.50m。其中以沿库周边线附近上升较明显，向后缘升幅逐渐减小，至后缘上游边界附近升幅仅40～50 cm。其三是地下水埋深发生不同程度的变浅，最明显的是沿库边线附近地下水埋深<1m的范围有所扩大。

8 临界埋深的确定

依据《水利水电工程地质勘察规范》（GB50487-2008），评价浸没的地下水临界埋深Hcr=土壤毛管水上升带的高度（m）+安全超高值（m）。

库区浸没影响主要有农作物、居民，农作物主要为小麦、水稻、油菜等，最大根系发育深度约0.60m，一般民房基础砌置深度为0.50m左右。故安全超高值取0.60m。龙井～下土城地块Q4低液限粘土毛细上升高度取0.63m，可以得到浸没的地下水临界埋深为1.23m，偏安全取1.25m。

9 水库蓄水后典型地块浸没预测

通过预测的地下水埋深Hm与浸没的地下水临界埋深Hcr进行比较，可判断模拟区的浸没程度。当地下水埋深Hm>Hcr时，不会产生浸没。当Hm≦Hcr时，将产生浸没危害，即建筑物地基承载力下降，威胁建筑物安全稳定；因过度潮湿使当地群众生活环境恶化；因盐渍化危及农作物的正常生长等。其中当Hm<0m时，将产生沼泽化，危害更大。

龙井～下土城地块蓄水后模拟区地下水埋深在0～3m，地下水临界埋深Hcr为1.25m，以此可以得到该地块的浸没范围。当库水位达到90m时，该地块模拟区浸没宽度（90米库水位至浸没边界）在60～780m之间不等，90m库边线长3300m，浸没面积约为0.60 km2。受灾对象以农田为主，其次沿库岸边线的村民和农田也受到一定程度的危害。