典型内陆湖区地下水数值模拟及其主控因子识别

于晓露，宋 健，林 锦，吴剑锋，吴吉春

（1.南京大学地球科学与工程学院，江苏 南京 210023；2.南京水利科学研究院，江苏 南京 210029）

在我国西北干旱地区，地下水资源是生态环境中最活跃的因素[1-2]，同时也是当地经济发展的基本依靠[3-4]。根据2020年中国水资源公报，全国水资源总量为31 601.2×108m3，内蒙古水资源总量为503.9×108m3，仅约占全国总量的1.6%，水资源短缺现象明显，严重制约当地社会经济发展。本次研究选取内陆河流域典型研究区—黄旗海盆地，位于内蒙古自治区乌兰察布市察哈尔右翼前旗境内，整体符合西北地区自然环境特征，即降水量小，蒸发量大，地表径流多季节性变化，流域对地下水资源依赖严重，地下水供水量约占总供水量的90%，导致当地地下水超采现象严重，对当地影响最为明显的是黄旗海湖水面积的变化。自1976—2001年黄旗海湖水面积下降三分之一，到了2010年，已经处于全年基本干涸状态，湖泊、湿地面积萎缩直接导致周边生态环境的退化[5]。黄旗海流域地下水动态监测信息缺乏，建立可靠的地下水模型是开展内陆河流域典型区地下水动态研究的重要基础。

地下水数值模拟技术是定量评价地下水水量的重要方法[6]，MODFLOW[7]作为最早开发的地下水数值模拟程序，历经多次更新后功能丰富，再加上开源、模块化的优势，当前被广泛运用于地下水运动及溶质运移、地表-地下水交互、海水入侵、地面沉降等问题的研究[8-11]。

敏感性分析是指模型参数对于模型结果不确定性的影响分析，具有量化输入的不确定性以及利用模型输出结果计算敏感性度量的功能[12-13]。敏感性分析通常分为局部敏感性分析（LSA）和全局敏感性分析（GSA）[14]。局部敏感性分析侧重针对某一参数改变对模型输出的影响，优点在于简单快捷，可操作性强，但难以表征多个参数对结果的整体影响程度[15]。全局敏感性分析又分为定性敏感性分析，如Morris 筛选法[16]等，以及定量敏感性分析，包括Sobol 方法、回归分析法、傅里叶振幅敏感试验（FAST）等。全局敏感性分析一并分析模型中的所有参数对模型输出的综合影响，并且可以深入分析参数间的相互影响, 但相应计算压力较大。随着敏感性分析方法在水文地质领域的不断发展，当前全局敏感性分析是识别地下水流模型影响因素与模型输出之间复杂关系的首选方法[17]。

对黄旗海盆地及其周边地下水的研究开始于1973年。自此以来，多个研究单位对周边水源地的地质条件、水文气象特征、土地利用状况、水资源承载力等方面进行了详细调查，发现研究区内持续存在缺水问题并定性给出导致黄旗海湖泊萎缩的可能因素，但未计算黄旗海湖泊与地下水在不同时期的相互交换量，进而难以准确分析导致湖泊萎缩的各类影响因子。所以，通过建立黄旗海盆地地下水流数值模拟模型，与实际监测数据相互佐证，才能定量分析地表水与地下水的相互转化关系，利用全局敏感性分析方法表征多个主要输入参数对模型输出结果的综合影响，并进一步分析获取黄旗海湖泊萎缩的主控因子，以期为黄旗海湖周边经济社会的可持续发展以及生态环境的逐步恢复提供理论依据，实现人与环境和谐共处。

1 研究区概况

黄旗海盆地位于内蒙古自治区乌兰察布市察哈尔右翼前旗，属典型内陆河流域，其主要地形为丘陵、高平原和盐湖平原，总流域面积为4 511 km2，见图1（a）。研究区位于阴山山地水文地质区，本次研究的主要供水含水层为第四系全新统（Q4）和更新统（Q3）孔隙潜水含水层组，见图1（b）。流域降水总量小且年内分配极不均匀，多年平均降水量为327.6 mm，蒸发量大，多年平均蒸散量（除黄旗海湖面）为237.93 mm。

图1 黄旗海流域地形及典型水文地质剖面图Fig.1 Map of the Huangqihai watershed (a) and a profile (b)

黄旗海流域是包括霸王河、泉玉林河等内陆河的最终汇入地，地表径流多季节性变化，且由于上游水库的建设，仅雨季或开闸放水农灌期内，河道有水渗漏补给地下水。黄旗海原属于桑干水系外流湖，由于陷落较深，隆盛庄河的出口被切断，演变为现在的内陆闭塞湖。研究区地下水的补给来源主要为降水入渗补给和周围玄武岩裂隙侧向补给，此外还接受农业回灌补给等其它人工补给；受地形地势影响，地下水径流主要为自西北向东南径流；主要排泄方式以人工开采和蒸散发为主。近些年，随着地下水人工开采量的增加，地下水水位埋深下降，黄旗海盆地地下水动态变化比较剧烈。

2 研究区模型构建

2.1 水文地质概念模型

结合研究区的地层岩性和含水层的供水条件，确定本次数值模拟的主要对象为第四系松散岩类潜水含水层组，组间存在少许弱透水层，但研究区内大部分弱透水层较薄甚至尖灭，各含水层具有较好的水力联系。从不同精度的地质调查与评价工作积累的钻孔资料中选择了45 个标志钻孔作为生成研究区水文地质模型顶底板的基础数据。黄旗海流域为四周高，中间低的碗状盆地，见图1（a）。依据研究区地质、水文地质条件以及相关水文资料，可将底部边界设为隔水边界；上边界设为开放边界，考虑降水入渗与蒸发等外部源汇项补给与排泄；基于研究区地下水水位长观数据，南北部边界通过达西定律计算，设置为给定流量边界；模型西侧有少量浅层地下水观测井，根据其观测深度可处理为给定水头边界；地下水径流较弱地区设置为隔水边界（图2）。研究区含水层结构和边界条件主要反映了黄旗海盆地的地下水流趋势，区内地下径流自盆地南北向中央汇聚，形成了对黄旗海的基流补给作用。

图2 水文地质概念模型Fig.2 Schematic model conceptualization

研究区含水层岩性主要为砂砾石、含砾粗砂、含砾中粗砂及粉细砂以及小部分玄武岩孔隙形成的裂隙含水层。初始水文地质参数取值参考了前人开展的39 个勘察孔抽水试验结果，采用克里格插值后生成的渗透系数作为模型的初值，并在模型识别过程中不断修正。

研究区含水层地下水补给项包括大气降水入渗、边界侧向流入、农业回灌补给、湖泊水入渗。研究区的主要排泄方式包括人工开采、蒸散发、地下水排泄至湖泊。大气降水数据来自TRMM 遥感数据，并结合水文地质调查报告以及土地利用数据，初步刻画出该区的降水入渗补给速率分区（图2）。蒸发数据基于MODIS 遥感数据、气象数据反演得到。基于研究区地下水蒸发量及 2015年土地利用分区图划分2010—2016年研究区潜水最大蒸发速率分区。模型中实际潜水蒸发量受含水层极限蒸发深度以及实时地下水水位控制。另外，研究区农业灌溉用水占比大，设置了田间灌溉补给系数计算农业回渗补给量。研究区内总共有开采井1 312 眼（图2），开采强度为 0.2×104～58.3×104m3/a，开采时间主要集中在4—8月。总开采量2010—2016年逐年递增，从3.292×107m3/a 增加到4.298×107m3/a。研究区内4 个水源地范围内抽水井较为密集且抽水量较大，对地下水水位变动影响较大。

2.2 数学模型

地下水流数值模型可概化成非均质各向同性的潜水二维非稳定流模型，相应的数学模型可表述为：

式中：K—含水层渗透系数/(m·d-1)；

b—潜水含水层底板标高/m；

W—大气降水、蒸散发、湖泊、人工开采和农业灌溉等源汇项/(m·d-1)；

μ—给水度；

t—时间/d；

h—含水层水位标高/m；

h0—初始水位/m；

h1—给定水头边界水位/m；

Ω—模拟范围；

Γ1—给定水头边界；

Γ2—侧向给定流量边界；

q—流量边界侧向单宽流量/(m3·d-1·m-1)，流入为正，流出为负，隔水边界值为0。

2.3 模型识别及验证

模拟区域总面积约为798 km2，网格剖分为100×100，共计 10 000 个有效单元。模拟期为2010年4月—2016年10月，根据当地春耕春灌特征以及降雨情况，总共划分为20 个应力期。其中识别期为2010年4月—2015年10月；验证期为2015年10月—2016年10月。

2.3.1 参数识别

模型计算过程中结合模拟情况反复调整渗透系数、降雨入渗补给系数等模型参数，以求模型模拟流场与真实流场拟合情况整体良好。最终模型渗透系数为7.61～70 m/d，具体分布情况见图3。给水度范围为0.16～0.20。具体各源汇项参数取值范围见表1。

表1 模型源汇项参数取值范围Table 1 Values assigned to the source and sink terms

图3 黄旗海盆地渗透系数分区图Fig.3 Map of hydraulic conductivity in the Huangqihai Basin

2.3.2 模型识别与验证

数值模型基于美国地质调查局的地下水水流模拟程序MODFLOW 通过试估-校正法进行参数识别过程。根据内蒙古地质环境监测总站提供的黄旗海盆地地下水水位观测资料，其中共有18 眼有效观测井用于模型的识别与检验，井位见图4。选取决定系数（R2）和均方根误差（RMSE）2 个性能指标作为模型拟合情况的统计指标，用于评估模型校准和模型验证。如图5 所示，模型的散点图和拟合曲线表明，在模型识别和验证阶段，观测值与模拟地下水水位之间的决定系数（R2）和均方根误差（RMSE）分别为0.994，0.870，表明该地下水模型拟合度较高。随后，选取6 组典型观测井绘制在整个模拟期内观测与模拟水位对比曲线图，观测孔91#和164#存在部分数据缺失。从图6中可知，整体来说，计算水位与实际观测值差距较小，大部分观测值与模拟值误差约为0.5 m，并且能够较好地刻画观测孔地下水水位动态变化趋势，表明该模型能够真实刻画区域含水层地下水水位变动情况。在春耕春灌期间，部分观测井水位与模拟水位部分偏差较为明显。主要是因为处于水源地，地下水开采强度变化剧烈；同时，根据总体水位变化趋势及规律，测量水位的周期可能存在偏差，无法获得本区域地下水水位的真实动态变化。

图4 观测孔位置分布图Fig.4 Locations of the observation wells

图5 模拟水位与观测水位的相关性Fig.5 Correlation between the computed and observed groundwater levels

图6 典型井观测水位与模拟水位对比曲线图Fig.6 Comparison of the observed and computed groundwater levels at the typical observation wells

3 水均衡分析

根据模拟计算的2010年4月—2016年10月地下水均衡组分分析，第四系潜水含水层接受补给量为2.550 7×108m3/a，主要补给项为降雨及灌溉补给入渗，达到1.509 4×108m3/a，占总补给量的60.4%；其次为地下水侧向补给1.040 6×108m3/a，占39.5%。含水层的主要排泄方式为人工开采地下水，开采量为1.662 7×108m3/a，占总排泄量的50.9%；其次为蒸散发，为1.204 6×108m3/a，占36.8%，具体水均衡组分评估见表2。研究区2010—2016年地下水开采量逐年增加，尤其是在每年的4—6月份的春耕春灌期内，导致含水层地下水超采现象严重。根据模拟结果绘制模拟期始末的地下水流场对比图（图7），可以直观看到地下水水位的下降态势，而湖周边地下水水位的下降也将直接导致湖水位降低。

图7 黄旗海盆地模拟期始末流场对比图Fig.7 Comparison of groundwater levels at the beginning and end of the simulation period in the Huangqihai Basin

表2 地下水均衡表Table 2 Water budget of the aquifer system

根据总流入量与总流出量对比可知，地下水的储蓄量总体呈减少趋势。2010—2016年地下水储量累计变化见图8（正值表示储量增加，负值表示储量减小）。整个研究期间地下水储量总体变化约为-1.5×108m3，年平均储量变化约为-2.188 3×107m3，尤其明显的是在降雨补给最为匮乏的2011年，年内地下水储量下降4.566 3×107m3，占整个研究期间地下水储存量总减少量的30%。湖泊渗漏量于2010—2014年间逐年递减，在2015年及2016年呈现微弱反弹趋势；地下水补给湖泊量整体呈上涨趋势；研究期间年平均湖泊渗漏量约为3.3×104m3，年均地下水补给湖泊量约为1.179×107m3。总体地表水与地下水转化关系为地下水补给湖泊水。

图8 黄旗海盆地地下水储量累计变化Fig.8 Groundwater storage changes in the Huangqihai Basin

基于1975—2015年湖水面积数据显示[5,18]（图9），黄旗海总体上呈现收缩下降态势，与近几十年全球性湖泊退化趋势相匹配，模拟期内黄旗海湖水位变化趋势与湖水面积变化趋势基本一致，拟合较好。研究采用皮尔逊相关系数分析模拟期内黄旗海湖水面面积萎缩驱动因子。在自然科学领域中，皮尔逊相关系数广泛用于度量2 个变量间的相关程度，计算值范围为-1～1。根据人工地下水开采量与黄旗海湖水位间的皮尔逊相关系数结果为-0.65。可知，两者在0.01（双尾）级别负相关性较为显著，所以模拟期间地下水水位的下降直接导致黄旗海湖的补给量不足，湖泊水位常年保持较低水平，难以保持地表水供水能力以及生态功能。

图9 黄旗海湖面面积年际变化Fig.9 Interannual changes of the lake area of the Huangqihai Lake

黄旗海盆地的地下水过度开采是该区水均衡呈现出负均衡以及黄旗海湖面面积持续萎缩以致消失的主要原因。因此，为了保护当地水资源并逐步恢复黄旗海湖泊及湿地面积，需要优化研究区地下水开采以及提高农业灌溉效率。

4 Sobol 法全局敏感性分析

4.1 Sobol 方法

Sobol 方法是基于方差分解的全局敏感性分析方法[19-21]，最早由Sobol 于1993年提出[22]，是一种应用较为广泛的定量全局敏感性分析方法。在本研究中，使用一阶和总阶Sobol 指数作为定义模型参数敏感性的指标。一阶敏感度（Si）描述单个参数对模型结果的影响；总阶敏感度（Sti）描述某一个参数与其他参数结合对模型结果的共同影响。

式中：S～i—所有不包含参数Xi的项，即与Xi互补的总方差。

基于本模型的特点，表3 总结了选取的8 个模型参数或输入量，包括水文地质参数（如渗透系数和给水度）以及关键的驱动数据（降水入渗补给、蒸发、湖床渗漏速率和地下水开采等源汇项）。考虑到大部分变量具有时空异质性，所以在敏感性分析过程中并不针对原始变量，而是给原始变量设定一个乘数，称为“调整因子”[23]，采用拉丁超立方抽样方法（LHS）在[0.8,1.2]中随机抽取参数样本，最后通过调整因子与地下水储量、湖水位、湖泊渗漏量、地下水补给湖水量等4 个目标函数的关系进行敏感性分析。

表3 敏感性分析的模型参数与输入量Table 3 Model parameters and their ranges of the input values for the sensitivity analysis

4.2 敏感度计算结果与讨论

采用地下水储存变化累计量和其余3 种输出变量在2010—2016年间平均值作为响应变量，采用基于全局敏感性分析的Sobol 法计算敏感度，分析在4 种目标函数下8 个输入参数的敏感性。图9 显示了参数样本为8 000 组情况下4 种目标函数对应的一阶和总阶敏感度及置信区间，其中部分敏感度可能出现负值，这是由于简化公式引起的数值误差。这种误差不影响分析结果，因为敏感度出现负值一般表示敏感度非常低[24]。

敏感度结果与置信区间见图10。除个别参数外，大部分参数的置信区间极值差小于0.1 且总体上不改变主控因子的敏感性排序，说明Sobol 敏感度计算结果可靠。一阶敏感度反映单个参数对目标函数的影响能力；总敏感度反映一个参数与其他参数共同对目标函数的影响能力；总敏感度和一阶敏感度的差值为参数的交互敏感度，反映了参数间相互耦合作用的间接影响，实际表征模型“异参同效”的可能性[20]。

图10 不同目标函数下参数敏感度结果及置信区间Fig.10 Sensitivity analysis results by the Sobol method with the confidence bound for different objective functions

图10（a）显示黄旗海地下水储量目标函数中人工开采与降水灌溉补给的一阶敏感度分别为0.379 6，0.262 6；总敏感度分别为0.381 3，0.439 9；交互敏感度分别为0.001 7，0.177 3。因此可知，地下水储量变化主要受人工开采和降雨灌溉补给2 个参数的控制，且降雨灌溉补给的交互敏感性较大，可能存在模型参数的“异参同效”性。

从图10（b）中可知，湖水位函数中渗透系数与侧向补给的一阶敏感度分别为0.443 7，0.168 0；总敏感度分别为0.709 6，0.332 7；交互敏感度分别为0.265 9，0.164 7。渗透系数与湖床渗漏速率2 个参数的敏感度显著高于其他参数，表明湖水位受两者影响并受到参数交互影响明显。

图10（c）显示黄旗海湖泊渗漏量目标函数中渗透系数与侧向补给的一阶敏感度分别为0.532 7，0.230 7；总敏感度分别为0.625 5，0.213 2。图10（c）与图10（b）对比可知，湖水位函数和湖泊渗漏量函数均受渗透系数与湖床渗漏速率2 个参数的显著影响，除主控参数外的其余参数敏感性排序也具有相似一致性。

图10（d）显示黄旗海地下水补给湖水量函数中人工开采与降水灌溉补给的一阶敏感度分别为0.526 0，0.180 7；总敏感度分别为0.572 7，0.349 5；交互敏感度分别为0.04，0.168 8。可知地下水补给湖水量主要受到人工开采和降雨灌溉补给2 个参数影响，进而与地下水储量目标函数相似。

在气候条件稳定的情况下，人工不合理开采地下水是影响水资源储量以及地下水补给湖水量的主控因子，这与前文中水均衡分析结果一致。现有的地下水开采强度使得研究区内的地下水资源不能实现可持续利用并加剧黄旗海的干涸速度，对湖泊及其生态环境造成较大威胁。因此从长远角度看，黄旗海水资源的利用要全流域统筹考虑，建立地下水资源优化管理模型，通过上游大坝按期放水及严格规范湖泊周边的人工取水量等切实措施，逐步恢复黄旗海湖水位，保护黄旗海湖泊湿地等珍稀生态资源，实现人与水的和谐相处。

5 结论

（1）构建的黄旗海盆地地下水流数值模型显示，长期观测井的地下水水位观测值与模拟值的决定系数（R2）为0.994，均方根误差（RMSE）为0.87m，且模拟水位能较好地反映观测孔水位动态变化趋势，说明该地下水水流数值模型能基本再现研究区地下水流运动变化状况，具有较好的可信度。

（2）模拟分析表明，第四系潜水含水层主要补给项为降雨及灌溉补给入渗，约占60.4%；其次为地下水侧向补给，约为39.5%。含水层的主要排泄方式为人工开采，占总排泄量的50.9%；其次为蒸散发排泄，约为36.8%。2010—2016年，湖泊渗漏量约为3.3×104m3/a，地下水补给湖泊量约为1.179×107m3/a，地表地下水转化关系整体呈现地下水大量补给湖泊水，且因地下水过度开采，研究区地下水储量累计亏空达1.5×108m3，这是黄旗海盆地地下水系统呈现为负均衡的核心原因。

（3）基于Sobol 分析法的参数敏感性分析表明，湖水位函数和湖泊渗漏量函数均受渗透系数与湖床渗漏速率2 个参数的显著影响。地下水储量函数和地下水补给湖水量函数主要受到人工开采和降雨灌溉补给2 个参数影响。说明超采地下水是造成黄旗海盆地地下水资源枯竭以及黄旗海湖泊面积萎缩的主控因子。

（4）研究区内观测井分布不均匀，大部分集中于黄旗海北岸水源地附近，而湖泊周边及南部区域缺乏长期地下水水位动态监测井，因此全区地下水模型精度仍有待提高。敏感性分析中忽略了各敏感参数的时序影响，未来需对此进行补充研究和完善。现有地下水开采模式会加剧黄旗海湖泊的干涸速度，不利于该区地下水资源的可持续利用。有必要严格控制地下水开采，提高农业灌溉效率，优化地下水开采方案，在保证当地社会经济发展的同时，逐步恢复黄旗海湖泊及湿地面积，促进人与自然和谐发展。