引言

民勤绿洲位于甘肃西北部武威市石羊河流域下游，干旱少雨，绿洲东、北、西3面都被沙漠包围。绿洲唯一的地表水来源为南上游红崖山水库。红崖山水库来水主要由上游6大河的水量转化而来。由于流域内经济发展，上游用水量增加，地表水被拦蓄或开发利用，致使民勤灌区内地表水来水量锐减⑴，直接导致地下水开采量增加，地下水抽水井数量从1980年的3 000眼增加到了2000年的9 010眼^[2]，超采严重。长期超采使得地下水位持续下降，生态环境恶化，并陷人恶性循环。民勤的生态危机引起了多方的关注和研究^[3-4]。许多学者就如何提升地下水位，恢复生态环境等问题进行了一系列的探讨研究。其中，马兴旺、霍再林等均曾借助Feflow软件使用早期(2000年前)的观测数据进行了民勤地下水位预测^{[2, 5-6]}，并分析了一些相关措施对水位的影响。近几年间，当地实施了一系列节水改造工程和流域调水工程，并采取了“关井压田”措施，以防止生态环境进一步恶化。种种措施后，当地的用水环境已经改变，因此，有必要对用水现状条件下地下水位动态变化进行重新预测。本文基于观测数据(2011-2012年），构建民勤绿洲地区的潜水-承压水三维数值模型，应用Feflow软件，对该地区地下水位进行预测分析，以期为未来水资源管理和开发利用提供参考和依据。

1 研究区概况

民勤绿洲处于武威盆地，按照地形、水利条件及历史习惯，当地又将该区分为坝区、泉山区和湖区3个灌区(图 1)。南至红崖山水库，北到青土湖，长约100 km。民勤地处内陆腹地，属典型的大陆温带干旱气候，多年平均降雨量110 mm，降水多集中于7~9月。因灌区内降雨稀少，降雨量多为5 _以下，不产生地面径流，降水除对绿洲边缘地带地下水有补给作用外，对农业生产作用不大^[7]。而用于农业灌溉的地表水唯一来源是红崖山水库调蓄供水。河水经跃进总干渠输人各干支渠，最后进人田间，用于农田灌溉。2001年后，景电二期工程即开始向民勤调水^[8]。2010年，西营河向民勤蔡旗专用输人渠道工程通水并开始稳定输水，改善了红崖山水库断面来流逐年锐减的局面。并且，红崖山水库为青土湖生态恢复供给一定的生态水量。另外，2008年到2010年间，当地进行了一系列节水改造工程，对渠系进行了改建衬砌，提高了地表水利用率^[9]。与此同时，针对地下水超采严重的情况，当地在2006年到2010年间实施了“关井压田”的措施，通过关闭部分机井来达到减少地下水开采量的目的。至2010年，民勤机井数量从9 010眼削减到了7 340眼，民勤地下水年均开采量由此前的5.20×10⁸ m³降至1.05×10⁸ m³。

研究区包括民勤绿洲核心区域。南部边界为红崖山水库，东、西、北为绿洲边缘。具体如下：西部以黑山头向北余脉及民勤绿洲与昌宁盆地地下水分水岭为界，北部沿北部诸山的山前地形走向确定，东北沿民勤绿洲与流沙接壤区为界，东南沿腾格里沙漠与绿洲接壤为界，区域总面积为2 868 km²。

2 潜水-承压水三维数值模型 2.1 水文地质概念模型

根据民勤绿洲水文地质构造，民勤盆地第四系含水层在地表以下120~150 m处由南向北存在多个总厚度约50 m的分布比较完整的弱透水黏土层^[10]。因而可将研究区含水层在垂向上概化为3层：上层为潜水含水层，底板高程从南边界1 240 m渐变为1 210 m; 中间层为弱透水层，厚度为50 m; 下层为承压含水层，承压含水层底板高程设定为1 060 m。承压水含水层以下视为基岩^[11]。根据当地地质资料，处在研究区南部的苏武山和北部的狼刨泉山是基岩抬升出露地表的产物，这两个山体下的岩层空间应视为不透水体。

2.2 地下水流数学模型

根据上述分析，将研究区地下水流概化为非均

$ \left\{ \begin{array}{l} S\frac{{\partial h}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {{k_x}\frac{{\partial h}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{k_y}\frac{{\partial h}}{{\partial y}}} \right) + \frac{\partial }{{\partial z}}\left( {{k_z}\frac{{\partial h}}{{\partial z}}} \right) + \varepsilon \\ \mu \frac{{\partial h}}{{\partial t}} = {k_x}{\left( {\frac{{\partial h}}{{\partial x}}} \right)^2} + {k_y}{\left( {\frac{{\partial h}}{{\partial y}}} \right)^2}-{k_z}\frac{{\partial h}}{{\partial z}} + P\\ h(x, y, z, t)\left| {_{_{t = 0}}} \right. = {h_0}(x, y, z)\\ {\left. {h(x, y, z, t)} \right|_{{h_{{\rm{osp}}}}}} = {h_{{\rm{sp}}}}\\ {\left. {h(x, y, z, t)} \right|_{\Gamma 1}} = {h_1}(x, y, z, t)\\ {k_n}\frac{{\partial h}}{{\partial \vec n}}\left| {_{\Gamma 2}} \right.{\rm{ = }}q\left( {x, y.z.t} \right) \end{array} \right. $

式中：S-承压含水层的储水系数，无因次；t-时间，d; h-地下水位，m; k_x，k_y，k_z-x，y，z方向上的渗透系数，m/d; ε-源汇项，流人为正，流出为负，m³; Ω-地下水研究区；µ-水层给水度，无因次；h₀-地下水初始水位，m; h_osp-泉口处位置，m；h_osp-泉口标高，m; h₁第一类边界水位，m; Γ₁一第一类边界；Γ₂-第二类边界；P-潜水含水层的蒸发、降雨、下渗等，d^-1; k_n-第二类边界处的侧向渗透系数，m/d; $\mathop n\limits^ \to $二类边界处的侧向法向量; q-渗流通量，m³/d。

2.3 三维数值模型的建立

研究区边界在GIS中绘制生成，导人Feflow中进行建模。模型采用三角形单元剖分，并对地下水位变化剧烈区域（主要是井渠灌溉结合的绿洲农田）进行有限元网格的加密。研究区剖分为484 206个三角形单元，共计324 284个结点。

民勤绿洲南部边界视为二类流量边界，在前人资料基础上^[10]，综合考虑调水工程以及水库扩容的条件下，水库蓄水量扩大约一倍，引起水位上升约一倍，利用达西公式确定其南边界渗流通量为0.24×10⁸ m³/a。西部边界定为零通量边界，北部为地下水流出边界，东部为地下水流人边界。其通量数据参见相关文献。

地下水在垂直方向上的源汇项主要考虑自然条件(降水、蒸发)和人类活动(抽取地下水、灌溉水回灌)等因素^[12]。根据文献[2]和文献[6]研究结果，草地、裸地等利用类型土地蒸散、降水和凝结水人渗补给地下水仅限于地下水埋深小于5m地区。目前，民勤地下水埋深均大于5 m，所以源汇项中不考虑降水、凝结人渗及草地、裸地上的蒸发，主要以地下水开采、渠系和田间灌溉水人渗为主。其中，地下水开采量以时间序列方式在井边界处进行赋值; 质各向同性中的三维潜水承压水非稳定流，采用的数学模型（含边界条件）为

$ \begin{array}{l} x, y, z \in \Omega, t \ge 0\\ x, y, z \in \Omega, t \ge 0\\ x, y, z \in \Omega, t \ge 0\\ 一类边界, t > 0\\ x, y, z \in \Gamma 1, t > 0\\ x, y, z \in \Gamma 2, t > 0 \end{array} $

(1)

渠系下渗量的计算不仅考虑了地表水灌溉水渠系人渗量，还考虑了青土湖生态泄水渠道的额外下渗，并且将下渗计算范围覆盖到支渠；田间下渗量则进一步考虑当地种植结构及作物实际灌水时序。

水文地质参数是地下水运动数值模拟不可或缺的输人项。根据有关水文地质资料，将研究区水文地质参数分为21个区域(图 2)。

为便于调节，将含水层、弱透水层和承压水含水层的参数分区设为一致，并依据文献资料对各个小区域赋予初值。每一层的水文地质参数最后都通过结果调参确定。

在Feflow中，模型的3层由切片隔开，每2个切片构成1个层。有限元网格、边界条件、初始条件均只对应于切片进行设定，而给水度、渗透系数等水文地质参数则对应于层进行设定。

2.4 模型识别与校验

采用2011年数据对模型进行识别，通过反复调试，最终确定模型选用的各层水文地质参数如表 1所示。使用2012年数据进行模型验证。经计算，大部分观测点模拟水位与实测水位误差小于0.7 m，模型满足精度要求。

3 模拟结果分析

以2010年12月实测水位为初始水位，采用上述地下水数值模型，预测民勤绿洲5 a、10 a后地下水水位变化。时间步长设为1 d。

3.1 不同地理位置地下水位变化

图 3为分别从坝区、泉山区和湖区各取一观测点绘制的地下水位随时间变化曲线图。通过分析研究区内不同位置地下水位变化过程，可以得到:相比2010年以前，民勤绿洲整体水位下降趋势变缓。其中，民勤绿洲腹部地带地下水位下降速度较灌区边缘和荒漠区边缘缓慢，这是由于近年来地表水来水量增加，增加水量主要集中补给在绿洲中心渠系密集地区，并且靠近水库的区域地下水位反而有上升趋势。泉山区由于缺少水库渗漏补给而地下水位下降速率高于坝区。下游湖区十一干渠附近区域由于青土湖生态泄水补给而地下水位下降幅度较小。

3.2 绿洲地下水位时空变化预测

图 4是民勤地区现状用水条件下2010年、2015年、2020年的地下水位分布图。可以看出，该地区整体水位趋于稳定，从坝区至湖区，地下水水位呈现递减的趋势。湖区地下水位偏低并非完全是因为地下水埋深较深，与湖区地势偏低也有关系。

为分析地下水漏斗发展趋势，利用GIS做出地下水埋深插值分布图（图 5)。由图 5可知，民勤绿洲地下水埋深值呈现全局增加，局部减小的趋势。上游坝区地下水埋深值大于下游湖区埋深值。坝区上游靠近水库区域由于截取了较多的地下水呈现出区域上升的趋势，而坝区中心位置的苏武，薛百地区的地下水漏斗随着时间推移呈现扩大趋势，地下水埋深值也随之增加。而灌区下游，在红沙梁与西渠镇，分别有一大一小两个地下水漏斗。随着青土湖下泄水的补给，2010-2015年，西渠的地下水漏斗已经恢复，红沙梁地区的地下水漏斗呈现缩小的状态。但是到2020年，随着地下水埋深进一步加大，红沙梁地区的地下水漏斗又呈现扩大趋势。总之，当地地下水仍然属于超采状态，地下水埋深持续增加，青土湖生态泄水虽能短时间恢复部分区域地下水水位，但若要持久地保持地下水水位仍需采取更多措施。

3.3 绿洲植物生存适宜区域面积变化

民勤地区的生态需水植被包括人工植被和天然植被两类。人工植被主要是以杨树、沙枣为主的农田防护林和以梭梭为主的人工固沙林。天然植被，只考虑绿洲与荒漠带1 km范围以内天然怪柳沙包和白刺沙包（占总天然植被的90%)。这些植物生长存活的地下水警戒埋深值为15 m^[13]。通过对地下水埋深插值图计算分析可得，2010年研究区植物适宜生存面积为1.91×10⁴hm², 占总面积的比例为6.8%; 2015年降至1.68×10⁴ hm², 占总面积的6.0%; 2020年降至1.56×10⁴ hm²，占总面积的5.6%;适宜植物生存的面积在不断减小。而其中，地下水埋深小于10 m的区域面积从2010到2020年减少了1 876 hm², 占适宜面积减少量的50%以上。

4 结语

利用建立的潜水-承压三维数值模型对民勤地区现状条件下地下水位进行了预测分析，结果表明，经过一系列综合措施，地下水位下降趋势虽得到有效遏制，但在未来5~10 a内，民勤地下水位依然整体以低幅度持续下降。因此，该地区仍需进一步采取相关措施解决研究区地下水位持续下降的问题。