基于GSFLOW的青土湖生态输水量-湖水面积关系研究

2022-09-21郭云彤崔亚莉邵景力

水文地质工程地质 2022年5期

郭云彤，周妍，崔亚莉，邵景力

（中国地质大学(北京)水资源与环境学院, 北京 100083）

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

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

图1 红崖山水库、民勤盆地与青土湖Fig.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 水体、芦苇面积及入湖水量变化Fig.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 青土湖水面面积与水位关系Fig.3 Relationship between the surface water area and water level of the Qingtu Lake

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

式中：S—水面面积/km2；

hl—湖水水位/m。

式（1）中R2分别达0.996 7 和0.999 9。

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]Fig.4 Coupling model between PRMS and MODFLOW[20]

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

式中：—第m时间步长第n次迭代中从贡献HRU 到湖泊HRU 的地表径流和壤中流体积/m3；

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

J—HRU 的计数；

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

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

式中：—第m时间步长第n次迭代中从湖泊HRU 到MODFLOW-2005定义的湖的体积流量/（m3·d-1）；

C'prms2mf—prms 的单位到modflow 中单位的换算系数；

—第m时间步长的湖泊HRU 上的降雨量/m；

—第m时间步长的湖泊HRU 上的蒸发量/m；

AlakeHRU—湖泊HRU 的面积/m2。

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

式中：—第m时间步长第n次迭代中穿过湖床至含水层有限差分单元中心的体积流量/（m3·d-1）；Klkbd—湖床的水力传导系数/（m·d-1）；

thicklkbd—湖床厚度/m；

thickaq—含水层厚度/m；

—湖床覆盖有限差分单元的面积/m2；

Kaq—靠近湖单元的含水层有限差分单元的水平或垂向渗透系数/（m·d-1）；

—第m时间步长第n次迭代中的湖水位/m；

—第m时间步长第n次迭代中有限差分单元靠近单元节点处湖泊的地下水水头/m。

3 结果与分析

3.1 模型建立与识别

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

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

图5 模型范围与空间离散Fig.5 Simulation area and spatial discretization

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

图6 湖区钻孔分布及地层剖面[19]Fig.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、决定系数R2为0.772。然而，由于青土湖区实际的地形相对平坦，但也存在微小的起伏，网格剖分的空间分辨率不足以反映某些微小的地形起伏变化，因此模型模拟的湖泊水位与实际水位存在一定的偏差，特别是当水位较低时，偏差相对较大。

图7 青土湖模拟水位与实际水位变化Fig.7 Changes in the simulated water level and actual water level of the Qingtu Lake

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

图8 观测孔分布图Fig.8 Distribution of the observation wells

图9 观测孔拟合图Fig.9 Fitting of the predicted water level and measured water level

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

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

图10 民勤盆地1960—2019年降水、蒸发量变化Fig.10 Variation of rainfall and evaporation from 1960 to 2019

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

图11 不同输水方案下湖水水位变化图Fig.11 Water level variations of the lake under different water conveyance schemes

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

图12 不同方案下V01 观测孔地下水水位变化图Fig.12 Water level variations of observation well v01 under different water conveyance schemes

生态输水量分别为3 100×104，4 500×104，6 000×104m3/a 时，预测2039年湖泊水均衡情况见表1。

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

表1 不同生态输水方案下2039年预测湖泊水均衡情况Table 1 Predicted lake water balance in 2039 under different ecological water conveyance schemes /104 m3

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

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

图13 不同输水量条件下最大水面面积变化图Fig.13 Maximum surface water area under different water conveyance schemes

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