我国西北干旱区水资源匮乏，随着社会经济的发展，区域用水需求增大，生产生活用水不断挤占生态用水，地下水水位下降导致植被退化、湿地面积衰减、生态系统退化等一系列生态问题^[1]。为缓解区域各部门用水矛盾，遏制生态环境恶化的趋势，近年来生态用水相关研究已成为当前西北干旱区水资源开发利用的重点。

青土湖是我国西北河西走廊三大流域之一——石羊河流域的尾闾湖^[2]，作为腾格里沙漠和巴丹吉林沙漠之间的生态屏障，青土湖湿地在防止沙漠合拢、遏制流域生态恶化趋势上具有重要作用。由于气候变化及人类活动的影响，青土湖于1959年完全干涸，湿地消失^[3]；自2010年开始，红崖山水库有计划地沿渠道向下游进行生态输水，青土湖开始形成季节性水面，周边生态环境明显好转^[4]。因此，研究青土湖生态输水量与湖水面面积的关系，对于确定合理的生态输水量尤为重要。石羊河的地表来水及调水工程的水汇入红崖山水库后，通过渠道供给下游民勤盆地用水（主要为农业灌溉用水）和青土湖生态输水（图1）。当前针对生态输水对青土湖生态环境的影响已有一些研究成果^[3,5-7]，然而这些研究多是结合遥感解译和定位观测的方法，对获得的数据进行统计分析得到输水量-湖水面积的规律，如果用这个统计规律外推，幅度偏大，其结果会存在着很大的不确定性。

图  1  红崖山水库、民勤盆地与青土湖

Figure  1.  Location of the Hongyashan Reservoir, Minqin Basin and Qingtu Lake

下载: 全尺寸图片 幻灯片

地表水-地下水耦合数值模拟方法是定量分析输水与湖区面积关系最有效的方法。当前地表水-地下水的耦合模拟模型大致分为独立型、联合型和集成型3类。独立型是在相对成熟的地下水模型的基础上拓展部分地表水模拟的功能，如模块化三维有限差分地下水流动模型（Modular Three-dimensional Finite-difference Ground-water Flow Model，MODFLOW）的湖泊（LAK）、河流（RIV）、蒸散发（EVP、EVT）和入渗补给（RCH）模块等，这些模块已经在实际地下水模拟研究中被普遍应用。联合型是将成熟的地表水和地下水模型通过一定手段进行连接，共同构建地表水-地下水的模拟系统，模型中各子系统按一定顺序进行独立计算，而子系统之间只进行单向或双向传输^[8-9]。集成型模型是一种完全耦合模型，该类模型将地表水和地下水作为一个系统，通过同时求解各个水文过程的控制方程描述地表水和地下水水分的交换和动态变化过程，每一步都有其具体的物理意义，机理性强^[10-11]。独立型模型虽加入了部分地表水模拟功能，但对大区域地表水流动过程的刻画仍不足；联合型模型较独立型模型机理性更强的同时，又不需要集成型模型中过多的参数，是当前地表水-地下水耦合模拟中最为实用的模型。

由美国地质调查局2008年公布的用于模拟地表水-地下水相互作用的模型（Coupled Groundwater and Surface-Water Flow Model，GSFLOW）是近年来应用较多的联合型耦合模型，该模型将降雨径流模型系统（Precipitation Runoff Modeling System，PRMS）和MODFLOW进行耦合，能够同时模拟气候、地表径流、地下潜流以及融雪、湖泊、溪流和湿地等与地下水之间的相互作用，已在国内外部分地区的地表水-地下水相互作用研究中得以成功应用^[12-14]。当前GSFLOW的模拟技术多应用于流域尺度的地表水和地下水关系的模拟中，多研究河流与地下水的相互作用。GSFLOW很好地将MODFLOW的湖泊（LAK）模块和地表水模块进行耦合，适用于本次针对青土湖的模拟研究。本次研究利用地表水-地下水耦合数值模拟的方法，定量化研究了流域中上游生态输水对下游尾闾湖水域面积的影响。

本文基于高分辨率的DEM数据、利用ArcGIS分析功能确定了水面面积与水位、蓄水量的转换关系；进而运用GSFLOW软件构建地表水-地下水耦合数值模型，通过模拟得到不同生态输水量情景下青土湖水面面积，论证青土湖的适宜生态输水量。

1.   数据来源及分析

1.1   数据来源

本次研究针对青土湖生态输水量-湖水面积的关系，以数值模型为基础，其中涉及的数据主要包括钻孔资料、补给项（降雨、生态输水数据）、排泄项（蒸发数据）以及模型校准所需的验证数据（水位数据）等。具体数据如下：

（1）研究区不同分辨率的高程数据，包括90 m×90 m的DEM数据及1 m×1 m的无人机遥感地形数据^[15]，低分辨率数据作为数值模型的顶板标高，高分辨率数据用于后续确定不同输水量下水面面积的实际分布。

（2）研究区2010—2019年多期Lansat5及lansat8及国产高分2号（GF-2）遥感影像数据，用于识别植被及水体面积。参考NDVI及MNDWI进行分类^[16-17]，利用ArcGIS对其分布进行矢量化处理，获得研究区年内较为连续的湖面及芦苇分布及面积信息^[18]。

（3）本次研究主要采用民勤站2010—2019年日尺度的蒸发和降水数据，气象数据来自中国气象数据网的中国地面国际交换站气候资料日值数据集（data.cma.cn），经处理后主要作为数值模型中的源汇项。

（4）生态输水数据来自民勤县水务局2010—2019年的统计资料；地下水水位数据来自位于青土湖的地下水水位观测井的记录值。

（5）含水介质参数和水文地质参数初值根据研究区地层岩性及前人研究成果确定^[19]。

综合以上资料，青土湖水面面积和芦苇面积的年内变化特征见图2。从数据较完整的2013—2018年看，水面面积增加趋势与输水时段基本重合，随着输水的进行，水面面积持续增大，直至输水结束，在水面蒸发的作用下水面面积随即减小。芦苇面积逐年上升，且随着多年生态输水的进行，芦苇面积趋于稳定。综上，青土湖区水面面积与芦苇面积的年内变化特征为：水面面积变化在生态输水期主要受生态输水影响，非生态输水期主要受蒸发影响；生态输水后的水面面积增大对芦苇面积增大有促进作用。

图  2  水体、芦苇面积及入湖水量变化

注：入湖水量及输水时段中纵坐标表示当年入湖水量。横线长短表示输水时长。

Figure  2.  Variations in the surface water area, reed area and water inflow into the lake

下载: 全尺寸图片 幻灯片

1.2   水面面积与水位、蓄水量的转换方法

建立水面面积与水位、蓄水量的转换关系是后续湿地地表-地下水耦合模拟验证的重要参考依据，转换方法具体如下：

（1）利用高分辨率的DEM数据识别湖泊区域；

（2）利用ArcGIS的Surface Analysis工具以0.1 m的水位差为单位，计算不同湖水位时对应的淹没面积及淹没面以下的体积；

（3）对得到的水位与其对应的面积及体积数据进行趋势分析，即可得到水面面积与水位、蓄水量的转换关系公式。

遥感解译结果为湖面面积的变化，需要建立水面面积与水位、蓄水量的转换关系，并将湖面面积转换为湖水水位，作为验证湖泊模拟结果合理性的验证数据。水面面积与水位、蓄水量的转换关系见图3。

图  3  青土湖水面面积与水位关系

Figure  3.  Relationship between the surface water area and water level of the Qingtu Lake

下载: 全尺寸图片 幻灯片

如图3所示，水位与水面面积的关系式可表示为：

式中：S——水面面积/km²；

${h}_{{\rm{l}}}$——湖水水位/m。

式（1）中R²分别达0.9967和0.9999。

2.   数值模拟方法

本次研究采用GSFLOW进行青土湖地下水-地表水的耦合模拟。GSFLOW将MODFLOW中使用的有限差分单元和PRMS中使用的水文响应单元（Hydrological Response Unit，HRU）进行空间链接是PRMS和MODFLOW模型耦合中的一个关键步骤。这个过程是通过生成重力储层（GRV）实现的，GRV的作用是进行HRU和有限差分单元间的水量传输。由于HRU和有限差分单元具有不同的空间范围，所以每个GRV的空间范围由HRU和有限差分单元的交叉部分确定。每一个GRV被分配一个唯一的标识号，每个URU和有限差分单元内可以有多个GRV。如图4中A部分所示，第一层为MODFLOW中的有限差分单元，第二层为PRMS的HRU，第三层显示了GRV分布。具体关系如图4中B部分所示，包括4个 水文响应单元、6个有限差分单元和9个GRV。在实际应用中，需要为每个GRV指定拓扑参数，将每个GRV与相应的HRU和有限差分单元联系。每个GRV的重力排水通过GSFLOW模块gsflow_prms2mf被添加到对应的有限差分单元中。类似地，每个有限差分单元的地下水排出量也可通过GSFLOW模块gsflow_mf2prms被添加到GRV中。

图  4  PRMS与MODFLOW耦合模式^[20]

Figure  4.  Coupling model between PRMS and MODFLOW^[20]

下载: 全尺寸图片 幻灯片

利用GSFLOW中湖泊（LAK）模拟功能对湖泊随生态输水的变化以及与地下水的交互进行模拟，GSFLOW中的湖泊在PRMS中表示为湖泊HRU，在MODFLOW-2005中表示为一组有限差分单元。在MODFLOW中直接在湖泊模块中输入作用于湖泊的降水量、蒸发量和地表径流数据 ^[21]。然而在GSFLOW中，这些过程以及通过土壤带的壤中流是在PRMS中计算的，因此LAK模块的这些输入变量应设置为0。流入湖泊HRU的地表径流及壤中流计算公式为：

式中：${V}_{{\rm{lakeHRU}}}^{m,n}$——第m时间步长第n次迭代中从贡献 HRU到湖泊HRU的地表径流和壤 中流体积/m³；

${F}_{J,{\rm{lakeHRU}}}$——HRUJ中为湖泊HRU提供地表径流 和壤中流的面积占总面积的十进制 分数，在GSFLOW的参数hru_pct_up 中进行定义；

$J$——HRU的计数；

$JJ$——某一段河段贡献地表径流和壤中流的 HRU的总数。

进入MODFLOW中湖泊的量表示为：

式中：${Q}_{{\rm{lakeHRUtolakemf}}}^{m,n}$——第m时间步长第n次迭代中从 湖泊HRU到MODFLOW-2005 定义的湖的体积流量/（m³·d⁻¹）；

${C}\text{'} _{{\rm{prms2mf}}}$——prms的单位到modflow中单位的换 算系数；

${P}_{{\rm{lakeHRU}}}^{m}$——第m时间步长的湖泊HRU上的降雨 量/m；

${EVAP}_{{\rm{lakeHRU}}}^{m}$——第m时间步长的湖泊HRU上的 蒸发量/m；

${A}_{{\rm{lakeHRU}}}$——湖泊HRU的面积/m²。

假设湖床存在于湖泊单元和含水层有限差分单元之间，并且具有与底层有限差分单元不同的特性，则湖泊和地下水之间交换量的计算公式为：

式中：${Q}_{{\rm{lakeleak}}}^{m,n}$——第m时间步长第n次迭代中穿过湖 床至含水层有限差分单元中心的体 积流量/（m³·d⁻¹）；

${K}_{{\rm{lkbd}}}$——湖床的水力传导系数/（m·d⁻¹）；

${thick}_{{\rm{lkbd}}}$——湖床厚度/m；

${thick}_{{\rm{aq}}}$——含水层厚度/m；

${A}_{{\rm{aq}}}^{m,n}$——湖床覆盖有限差分单元的面积/m²；

${K}_{{\rm{aq}}}$——靠近湖单元的含水层有限差分单元的水 平或垂向渗透系数/（m·d⁻¹）；

${h}_{{\rm{lake}}}^{m,n}$——第m时间步长第n次迭代中的湖水位/m；

${h}_{{\rm{aqfdc}}}^{m,n}$——第m时间步长第n次迭代中有限差分 单元靠近单元节点处湖泊的地下水水 头/m。

3.   结果与分析

3.1   模型建立与识别

依据研究区的地质和水文地质条件等，确定本次模型模拟范围包括青土湖及其周边部分沙漠地区，南部以隐伏断层为界，设定为隔水的零流量边界；东北侧以沙漠边缘为界，设定为通用水头边界；考虑到裸土区地下水的极限蒸发埋深，以地下水埋深2.5 m为界，西北部以等水位线1309 m和1311.5 m为界，西南部以1307.5 m为界，设定为定水头边界；总面积约159 km²。模拟区的补给项主要为降水补给和人工生态输水补给，其中上游红崖山水库对青土湖的人工输水是整个区域的主要补给来源。模型的主要排泄项为湖水及地下水的蒸散发。模拟期为2015年1月1日—2019年12月31日，1 d为一个应力期进行模拟。

地下水模型区剖分为500 m×500 m的网格，如图5（a）。研究区地表标高采用90 m×90 m的DEM数据及1 m×1 m的无人机遥感地形数据拼接作为研究区的地表高程使用。

图  5  模型范围与空间离散

Figure  5.  Simulation area and spatial discretization

下载: 全尺寸图片 幻灯片

含水层结构如图6所示，湖区西部、东南部及南部均有黏粒含量高的沉积物分布，透水性差，东北部主要以沙土为主。为刻画含水层水文地质条件，反映湖水与地下水间的补排特征，将整个含水系统划分为一层，含水层厚度定为15 m。

图  6  湖区钻孔分布及地层剖面^[19]

Figure  6.  Borehole distribution and stratigraphic profile of the Qingtu Lake^[19]

下载: 全尺寸图片 幻灯片

地表水模型部分，对研究区的地表空间、地下空间及模型接口进行划分。研究区内地形相对平坦，降水稀少，地表产汇流较少，因此将研究区整体划分为陆地HRU和湖泊HRU，图5（b）；利用GSFLOW模型的GRV接口将地表HUR和地下水模型的有限差分单元网格连接起来，共划分769个GRV，图5（c）。

利用民勤县气象站监测的2015—2019年降雨量数据、气温数据作PRMS的数据文件，对研究区的降雨入渗和蒸散发进行计算；利用输水渠道模块SFR2模拟向青土湖的渠道输水，利用湖泊模块LAK模拟湖泊的变化。

青土湖水位变化拟合效果见图7。模型模拟水位变化与遥感解译水位变化趋势基本一致，模型模拟结果基本上能够反映湖水实际的水位变化情况，纳什指数NSE为0.766、决定系数R²为0.772。然而，由于青土湖区实际的地形相对平坦，但也存在微小的起伏，网格剖分的空间分辨率不足以反映某些微小的地形起伏变化，因此模型模拟的湖泊水位与实际水位存在一定的偏差，特别是当水位较低时，偏差相对较大。

图  7  青土湖模拟水位与实际水位变化

Figure  7.  Changes in the simulated water level and actual water level of the Qingtu Lake

下载: 全尺寸图片 幻灯片

模拟区观测孔位置如图8所示。为了验证模型模拟效果，2018年和2019年观测的实际水位与模拟水位进行对比（图9），结果表明，观测孔的模拟水位变化与实测值拟合较好。

图  8  观测孔分布图

Figure  8.  Distribution of the observation wells

下载: 全尺寸图片 幻灯片

图  9  观测孔拟合图

Figure  9.  Fitting of the predicted water level and measured water level

下载: 全尺寸图片 幻灯片

3.2   生态输水方案预测与结果分析

生态输水量是青土湖湖面变化最主要的影响因素，为研究生态输水量对湖面面积变化及地下水水位变化的影响，分别设置不同生态输水量：2 000×10⁴~4500×10⁴ m³/a每隔100×10⁴ m³/a设置1个输水方案，4500×10⁴~6000×10⁴ m³/a每隔500×10⁴ m³/a设置1个输水方案。输水时段均在8—10月，共设置29个输水方案。运用已建立的地表水-地下水流耦合模型，预测未来20年模拟区内地表水和地下水的变化情况。模型的结构、水文地质参数及边界条件均保持不变，模型模拟期设置为2020年1月1日—2040年12月31日。模型的降雨蒸发项需输入逐日数据。本次研究搜集到的逐日降雨蒸发数据为2010—2019年。从蒸发数据看，2000年后研究区年蒸发量在1876~2924 mm范围内波动（图10）。故在预测模型中降雨蒸发项取对应2010—2019年每日的多年平均值，年总蒸发量约2300 mm。

图  10  民勤盆地1960—2019年降水、蒸发量变化

Figure  10.  Variation of rainfall and evaporation from 1960 to 2019

下载: 全尺寸图片 幻灯片

不同输水量情景下预测湖水水位变化情况见图11。随着生态输水量的增大，水位变化呈现由大到小的过程。预测初期5~6年，湖水水位快速增高，变化幅度与生态输水量的大小有关；6年后湖水水位缓慢变化。当生态输水量为2 000×10⁴ m³/a时，湖水水位逐渐下降，最高水位保持在1210.4 m（平均水位1210.0 m），对应水面面积6.68 km²；当年输水量保持在3 100×10⁴ m³/a时，与现状模型2019年相比，水位变化不大，最高水位保持在约1212.2 m（平均水位1211.7 m），说明现状3100×10⁴ m³/a的生态输水量可以保证青土湖维持当前水面面积的生态需水。随着输水量的增大，到2040年，生态输水量达6000×10⁴ m³/a时，湖泊最高水位可达1217 m（平均水位1215.9 m），年内变化幅度达2.45 m。

图  11  不同输水方案下湖水水位变化图

Figure  11.  Water level variations of the lake under different water conveyance schemes

下载: 全尺寸图片 幻灯片

由于各观测孔水位变化趋势相近，故以V01孔为例（图12）显示不同方案下地下水水位的变化情况，该位置地表高程约为1 211 m，见图12。与湖泊水位变化相似，随着生态输水量的增大，预测初期5~6年，地下水水位快速增高， 6年后地下水水位缓慢变化。当输水量保持在3100×10⁴ m³/a时，地下水水位的年际变化不大，与现状年相比基本保持稳定，最高水位维持在约1210.8 m，最低水位维持在约1210.6 m；随着生态输水量的增大，湖泊对地下水的补给量逐渐增大，地下水的稳定水位逐渐提升，当水量达4000×10⁴ m³/a时，地下水水位达到地表高程；当输水量达5000×10⁴ m³/a时，水位常年高于地表以上；当输水量达6000×10⁴ m³/a时，最高水位可达1211.3 m，高出地表面约0.3 m。

图  12  不同方案下V01观测孔地下水水位变化图

Figure  12.  Water level variations of observation well v01 under different water conveyance schemes

下载: 全尺寸图片 幻灯片

生态输水量分别为3100×10⁴ ，4500×10⁴，6000×10⁴ m³/a时，预测2039年湖泊水均衡情况见表1。

表  1  不同生态输水方案下2039年预测湖泊水均衡情况

Table  1.  Predicted lake water balance in 2039 under different ecological water conveyance schemes /10⁴ m³

生态 输水量	补给项			排泄项		均衡量
	降水量	生态入湖输水量		蒸散发量	湖泊补给地下水量
3100	147	2170		2550	14	−247
4500	162	3150		2810	463	40
6000	165	4200		2862	1518	−15

下载: 导出CSV 

| 显示表格

如表1所示，当生态输水量小于4500×10⁴ m³/a时，随着生态输水量的增加，水面面积不断增大，水面蒸发量增大幅度明显，是湖泊的主要排泄项。当生态输水量大于4500×10⁴ m³/a时，随着湖泊水位的提高，地下水和地表水位差增大，湖泊向地下水的排泄量增大，湖面面积增大幅度很小，水面蒸发量增量有限。

3.3   青土湖生态输水适宜量确定

随着生态输水量的增大，湖面面积也逐渐增大，但增大的幅度有所变化。将水位转化成湖水的水面面积，则生态输水量与水面面积及水面面积变化率的关系见图13。随着生态输水量的增大，水面面积增大情况大致可以分为3个阶段：当输水量为2000×10⁴~3700×10⁴ m³/a时，水面面积随输水量增大而增大，基本为线性关系，面积变化率相对稳定，保持在8122~10796 m²/（10⁴ m³·a）；输水量为3700×10⁴~4500×10⁴ m³/a时，水面面积增大幅度减缓，面积变化率逐渐大幅度减小至2000 m²/（10⁴ m³·a）；当输水量大于4500×10⁴ m³/a时，随生态输水的增大水面面积增大幅度很小，特别是当生态输水量大于5500×10⁴ m/a时，面积变化率趋近于0，水面面积几乎稳定在26 km²。从维持当前生态水面面积不至减小的前提看，生态输水量不宜低于3100×10⁴ m³/a。另外从生态输水的效益考虑，生态输水不宜高于4500×10⁴ m³/a。因此，青土湖生态输水适宜量为3100×10⁴~4500×10⁴ m³/a，维持湖面面积大致为16.27~26.60 km²。

图  13  不同输水量条件下最大水面面积变化图

注：面积变化率指每1×10⁴ m³/a生态输水引起的水面面积（m²）变化。

Figure  13.  Maximum surface water area under different water conveyance schemes

下载: 全尺寸图片 幻灯片

针对生态输水与青土湖水面面积的变化研究已有一些研究成果，如利用统计分析和水均衡分析方法对2010年生态输水以来的生态输水数据、气象数据和水面面积变化数据进行分析，结果显示生态输水量在3145×10⁴ m³/a时可以基本维持当前的最大水面面积13.43 km^2[22-23]。而本次研究的结果显示生态输水量3100×10⁴ m³/a时水面面积可维持16.27 km²，结果与统计分析及水均衡分结果相比，相同生态输水状况下，维持的水面面积偏大。这可能是由于数值模拟得到的水面面积是生态输水多年累计达到的较为稳定的水面面积，略大于短期生态输水形成的水面面积。

4.   结论

（1）依据不同生态输水方案的模拟结果，当前3100×10⁴ m³/a的生态输水量可以保证青土湖维持当前水面面积的生态需水。

（2）生态输水量与水面面积变化的关系大致可分为3个阶段。结合该关系，同时考虑向青土湖生态输水的效益及保证当前水面面积的需求，红崖山水库向青土湖的生态输水适宜量为3100×10⁴~4500×10⁴ m³/a，湖面面积大致可维持在16.27~26.60 km²。

致谢：感谢清华大学胡宏昌教授及首都师范大学李浩乾硕士在地表水-地下水耦合数值模拟技术方面提供的支持，感谢马瑞教授、陈喜教授共享青土湖地形数据、地下水水位数据，感谢民勤县水务局提供的详细生态输水数据！