基于GSFLOW的镜湖湿地地表水与地下水耦合数值模拟

2022-05-23郜会彩肖玉福胡云进陈柳安周如杰

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

郜会彩，肖玉福，胡云进，陈柳安，周如杰

（1.绍兴文理学院土木工程学院,浙江绍兴 312000；2.浙江省岩石力学与地质灾害重点实验室,浙江绍兴 312000；3.浙江省山体地质灾害防治协同创新中心,浙江绍兴 312000）

湿地生态系统通过其特殊的水文物理特性，以水循环为纽带，削减地表径流、降低流速和补给地下水发挥其水文调蓄功能[1]。湿地地表水与地下水之间的交换转化，一直是水文地质研究中的热点问题[2]。以往国内外对河流的地表水-地下水交换转化研究较多，而对于水体相对静止或流动较缓滞的湖泊和湿地研究相对偏少[3]。目前研究湿地地表水与地下水耦合作用的方法主要有实地调查法（直接测量法、示踪剂法、水温-热量传递法、地球物理探测法）和数值模拟法[4]。其中湿地地表水与地下水耦合作用的数值模拟有利于理解湿地水文过程机理和预测水文情势变化，目前湿地地表水与地下水耦合作用模型的研究大多借助已有的地表水与地下水耦合模型[5]。根据模型耦合方式不同，可分为松散耦合（“四水”转化模型）、半松散耦合（SWATMD、GSFLOW 等）和紧密耦合（MIKE-SHE、IGSM、LL-II 等）[6]。一般湿地内地形起伏较小，由于湿地植物的作用，造成湿地中水流速度缓慢，水力持留时间长[7]。湿地表层流复杂，周边高地或其它水体以坡面漫流、片流、明渠流、漫滩流等多种形式进入湿地，很难建立精确的数学模型去模拟[8]。冯夏清[9]在SWAT 模型基础上重新构建了湿地模块来模拟扎龙湿地水文过程，实现了流域模型与湿地水文模型的松散耦合。Rahman 等[10]采用SWATrw 模型探讨了Barak-Kushiyara 河洪泛湿地地表水与地下水之间的水量交换，结果表明水文特征是影响湿地地表水补给的主要因素。现有模型中GSFLOW 可以将湿地作为独特的水文响应单元，重点考虑湿地水循环过程，模拟精度高、适用性强[11]。从已有研究来看，GSFLOW 已被应用于流域地表水与地下水耦合作用研究和流域水资源配置研究[12-14]；为帮助用户一站式构建模型，田勇等[15-17]基于GSFLOW 开发了具有国际先进水平的生态水文集成建模软件系统Visual HEIFLOW（VHF），并以黑河流域和韩国Miho 流域为例，验证了在大流域尺度下的适用性；但未见其在湿地地表水与地下水耦合模拟中的应用。

湿地中有大量的表生植物和沉水植物，阻碍湿地地表水流动，使得湿地的物理结构和化学环境更加复杂、异质和多变[4]；并且湿地土壤是固相和液相组成的疏松多孔体系，多年植物生长沉积形成了独特的海绵状土壤层。因此，本文针对目前研究的不足，考虑湿地海绵层作用，对湿地明渠流和片流同时进行模拟，利用GSFLOW 模型构建湿地地表水与地下水耦合模型，并将其应用于镜湖湿地实例中，检验模型可靠性，以期为准确评估湿地地表水与地下水耦合作用提供理论依据。

1 GSFLOW 模型原理

GSFLOW 模型耦合了PRMS 和MODFLOW-2005两个模型，分别进行湿地地表水文模拟和湿地地下水流模拟（图1）。地表部分将湿地划分成许多个水文响应单元（HRU），地下部分则离散成许多个有限差分单元，在地表水与地下水相连的河段，建立重力水库（GRV）将地表部分的HRU 与其对应的地下部分有限差分单元连接起来[14]，利用达西定律，根据水头差计算地表水与含水层的水交换量。

图1 GSFLOW 结构模型示意图[13]Fig.1 Structural framework of the GSFLOW model [13]

1.1 地表产流模拟

对于内陆湿地，降水首先被冠层截留，到达地面后一部分下渗到土壤成为地下水，其余的水沿着斜坡形成漫流，即其产流模式一般为超渗产流。GSFLOW模型将湿地分为透水区和不透水区分别进行计算[13]，每个HRU 的透水面积和不透水面积比例需事先给定。在HRU 透水区域采用超渗产流机制计算地面径流及水分下渗过程，其超渗产流利用贡献区域的概念计算产流面积[14,18]。在不透水区域主要计算地表径流量和蒸发量。

1.2 地表径流模拟

由于绝大部分湿地地势平坦，水流在平面大范围的自由表面流动，水流在水平方向的运动尺度远大于垂向，所以假定水流在垂向均匀分布，属于浅水流动。GSFLOW 中对地表水流运动的描述采用扩散波近似的圣维南方程，其中用来计算水流在矩形或棱柱形河槽中的运动方程如下：

qstr——单位长度河道总的侧向入流或出流量/（m2·d-1）；

xchannel——沿河道方向的河道长度/m；

t——时间/d。

对于明渠流，GSFLOW 采用曼宁公式计算流量：

式中：Q——河流流量/（m3·d-1）；

A——过水断面面积/m²；

Cm——转换常数，国际单位制中值为1；

R——水力半径/m；

S——明渠的坡度；

n——糙率。

选用湿地植物综合糙率公式[19]确定n，如下式所示：

式中：H——植物高度/cm。

其中沉水植物是湿地生态系统的重要组成部分，其生长特征和生物量积累随湿地水位变化。通过对绍兴镜湖湿地现场调研测量，湿地沉水植物平均高度为3 cm，故n为0.066。

对湿地而言还普遍存在片流形式，对片流的模拟详见下文2.1 节。

1.3 蒸散发模拟

蒸散发是湿地水循环的重要环节，也是区别于河流、湖泊水均衡的重要特点[5]。湿地地表蒸散发量是影响湿地水热平衡的主要因素，也是能量和水分的主要消耗途径。湿地植被供水良好，可视为植被实际蒸散发等于其潜在蒸散发。本文采用1998年世界粮农组织（FAO）推荐的Penman-Monteith 公式计算潜在蒸散发[20]，该公式受季节、站点情况影响小，模拟效果好。公式中作物系数基于假定的参考作物，各项系数和变量结合利用气象数据和文献[21]计算得到。

1.4 地下水流模拟

湿地地表水与地下水交互渗流分为非饱和流与饱和流[5]。非饱和渗流常以多孔介质流体达西定律与质量守恒定律为基础形成二相渗流运动的微分方程进行模拟[22]。在GSFLOW 中非饱和渗流使用垂向一维Richards 方程运动波近似形式，公式如下：

式中：θ——含水率；

Lz——垂直运移距离/m；

K(θ)——水力渗透系数/（m·d-1）；

i——地下水蒸发速率/（m·d-1）。

饱和渗流过程将湿地地表水与地下水视为连续的整体，符合多孔介质中潜水渗流运动的Boussinesq方程，如式（5）所示：

式中：Kxx——x轴方向水力渗透系数/（m·d-1）；

Kyy——y轴方向水力渗透系数/（m·d-1）；

Kzz——z轴方向水力渗透系数/（m·d-1）；

x、y、z——空间坐标；

Ss——储水率/m-1；

W——源汇项/d-1。

1.5 地下水的补给与排泄

湿地地表常年有积水，地表水是地下水的稳定补给来源。湿地地表水与地下水的交互作用主要受湿地植被、地形坡度和湿地土壤的含水量的影响。在PRMS 模型中，土壤层分为毛细水库、重力水库和优先流三个水库[14]，这些水库反映不同的土壤含水量（图2）。地下水补给与排泄的水量通过重力水库转换。当满足时，地下水排泄到重力水库，反之亦然。计算方程如式（6）所示：

图2 PRMS 模型土壤层各水库的输入和输出关系图Fig.2 Relationship between the input and output of each reservoir in the soil layer of the PRMS model

CNDsz——重力水库区域的导水系数/（m2·d-1）；

celtop——有限差分单元格的顶板海拔/m；

Dusz——土壤带底部的地形起伏变化深度/m。

2 湿地片流和湿地海绵层模拟

2.1 湿地片流模拟

降水到达湿地表面后，部分水流充填洼地和地表下渗，大部分直接形成湿地地表径流。明渠流和流过浓密植被的片流是湿地主要的地表径流形式[23]（图3）。在湿地植被茂密的区域，地表水流的水平速度较低。当水位继续升高，湿地表面完全被水淹没，在一定的水力坡度下沿沟渠方向产生表面流，其速度一般比片流快。片流的方向和速度受地形坡度、水深、植被类型、降水、蒸散发等作用，其水文过程的模拟和计算相对复杂。本文采用圣维南二维扩散方程[24]模拟湿地片流，见式（7）（8）：

图3 湿地主要地表水流成分示意图[24]Fig.3 Schematic map of the main surface water flow components in wetland[24]

式中：Sfx、Sfy——分别为x、y方向的摩擦坡度；

qx、qy——分别为x、y方向的片流通量/（m2·s-1）；

P——降水率/（m·s-1）；

I——入渗速率/（m·s-1）；

ET——片流蒸发速率/（m·s-1）；

Vw——其它水源或汇的流速/（m·s-1）。

2.2 湿地海绵层模拟

湿地植物根系与腐朽植物的分解很困难，容易在潮湿缺氧的条件下形成泥炭[25]。在水淹的条件下，湿地地表的泥炭藓、植物凋落物不易腐烂，在湿地底部长期积累堆积形成一层柔软疏松结构，仿佛巨大厚重的海绵[26]，与泥炭层一起构成了湿地海绵层（图4、图5）。该层突出特点是有较大的饱和含水量和持水量，孔隙较大，孔隙度72%～93%，既能蓄水又能透水[27]，在湿地地下水补给和排泄过程中具有重要意义。本文通过GSFLOW 的优先流水库模拟湿地海绵层，在设置优先流阈值时，输入湿地海绵层对应的孔隙度。对于湿地海绵层渗透系数，本文通过野外现场竖管法（图6）来进行测定[28]，该方法简便易行，精度高。在选取的试验点上将直径6 cm 的PVC 竖管垂直插入湿地，尽量不破坏湿地海绵层的原有结构。湿地海绵层垂向渗透系数计算公式为：

图4 湿地地表水与地下水相互作用地点及湿地海绵层位置 [4]Fig.4 Location of wetland surface water and groundwater interaction and location of the wetland spongy layer [4]

图5 现场采集湿地海绵层样本Fig.5 In situ collection of wetland spongy layer samples

图6 竖管法测试湿地海绵层渗透系数示意图Fig.6 Measurement of the hydraulic conductivity in the wetland spongy layer using straight standpipes

式中：Kv——垂向渗透系数/（m·d-1）；

Lv——管中湿地海绵层厚度/m；

h1、h2——分别为测管中t1、t2时刻水头值/m。

测量结果如表1所示，经过计算得到湿地海绵层渗透系数的范围为15～45 m/d。

表1 湿地海绵层垂向渗透系数Table 1 Vertical hydraulic conductivities in the wetland spongy layer

3 应用实例

3.1 研究区概况

研究区（绍兴市镜湖国家城市湿地公园）位于浙江省绍兴市镜湖新区，面积为5.13 km²。研究区数字高程模型（DEM）如图7所示。区内河流属萧绍曹运河水系，河湖密布，水网交织。研究区处于江山—绍兴深断裂东北端两侧，晚燕山早期白垩纪构造沉积盆地内[29]，区内构造被大面积的第四系松散堆积层所覆盖，未见大的断裂和皱褶构造，岩石出露区主要发育节理裂隙。根据浙江省水文地质勘察资料，将研究区含水层等厚划分为三层：上层为湿地海绵层，厚度为0.5 m；中间为潜水含水层，厚度为3 m；下部为稳定分布的粉质黏土，厚度为45 m。

图7 研究区DEM 图Fig.7 DEM map of the study area

3.2 数据收集与整理

GSFLOW 模型的主要输入数据包括：DEM 数据、土地利用数据、土壤类型、气象驱动数据、水文数据及土壤属性数据，各数据来源如表2所示。

表2 GSFLOW 模型所需数据类型Table 2 Data types for the GSFLOW model

借助ArcGIS10.2，基于研究区DEM 数据，划分流域。将土地利用类型重分类，根据研究区现场情况概化为建筑用地、绿地、水域和山林地等四种；对土地利用类型进行编码赋值（图8）。由于研究区域范围较小，为便于区分将土壤类型概化为黄棕壤、高钙土、褐土、潮土和其他等五种类型（图9）。

图8 研究区土地利用类型分布图Fig.8 Zonation of the land use types in the study area

图9 研究区土壤类型分布图Fig.9 Zonation of soil types in the study area

3.3 数值模型构建

根据研究区的地形、地貌、水文地质概况进行空间离散。按照DEM 图（12.5 m×12.5 m）和河网水系进行HRU 划分，一共生成4 632 个HRU。本次模拟的地下部分主要是第四系潜水含水层，先给定一个初始渗透系数和补给率，经过率定得到最终渗透系数和补给率。湿地地下水模拟采用MODFLOW-2005，对整个研究区进行网格剖分，共分为5 599 个单元格。

研究区上边界为地表面，接受降水入渗补给、河流补给、蒸发、排泄等水量交换，在GSFLOW 中处理为可自由移动界面。下边界是凝灰质砂岩基岩，处理为隔水边界。对于模型的四周边界，高低起伏较大的山坡处取至分水岭，作为隔水边界；多数边界附近地势平坦，为减少计算工作量，根据研究区DEM 图、水域分布图及实测水文情况，做如下处理：河道横断面较宽的河流，常年水位较恒定，作为定水头边界（图10的Γ1）；河道横断面较窄的河流，常年流量较稳定，作为定流量边界（图10 的Γ2）；其它周界由于水量交换较少，按隔水边界处理（图10 的Γ3）。

图10 研究区边界条件示意图Fig.10 Schematic diagram of the boundary conditions in the study area

尽管研究区面积较小，但由于渗透系数是影响湿地地表水与地下水交互作用的重要因素，因此根据前人研究结果及《浙江省绍兴市环境地质调查评价报告》等，结合研究区水文地质条件、钻孔资料和工程地质勘察资料，对地下水模型中水平渗透系数Kh、垂直渗透系数Kv和储水率Ss进行详细分区（图11）。表3 为各分区水文地质参数取值范围。

表3 研究区各分区水文地质参数取值范围Table 3 Range of hydrogeological parameters in each subdomain of the study area

图11 研究区水文地质参数分区示意图Fig.11 Zonation of hydrogeological parameters in the study area

3.4 模型参数率定及验证

以2005—2019年逐日降水量、最低和最高气温、蒸发量等气象数据作为驱动数据，对模型参数进行率定和验证。考虑到模型初始值对模拟结果的影响以及湿地现场水文资料的限制，将模拟期（2005年1月1日—2019年12月31日）划分为两个阶段：阶段1 从2005年1月1日—2014年12月31日，为湿地地表水与地下水耦合模型参数的率定期；阶段2 从2015年1月1日—2019年12月31日，为模型参数的验证期。在模型率定过程中，采用相关系数评价模型率定结果。

在模型率定过程中，受湿地独特水力特征和气候条件的影响，土壤表层饱和含水率、快速壤中流非线性系数、线性贡献面积算法指数、HRU 不透水区最大滞留储量等参数为湿地地表水与地下水耦合模型的敏感参数，在参数变动范围内调整参数值，来率定模型使其达到良好的模拟效果。利用率定好的参数，选取研究区内地下水位对模型模拟结果进行验证（图12）。验证期湿地地表水与地下水耦合模型平均水位与模拟区水文站观测地下水位的相关系数R2=0.671。将模拟结果与实测湿地河道水位进行比对（图13），采用均方根误差（RMSE）来衡量模拟水位和实测水位两个数据序列之间的误差，计算可得RMSE=0.086，表明模拟水位与实测水位误差很小。

图12 湿地地下水水位对比图Fig.12 Comparison between the observed and calculated groundwater levels in the wetland

图13 湿地河道2019年模拟水位与实测水位对比图Fig.13 Comparison between the observed and calculated surface water levels in the wetland

从以上率定和验证结果可知，湿地植被茂盛、地势起伏较小，片流和明渠流是湿地地表径流的主要形式，考虑湿地多年植物生长沉积形成的独特海绵层的调蓄和渗透作用，将湿地片流与明渠流同时进行模拟是合理可靠的。

3.5 研究区水均衡分析

为进一步揭示湿地地表水与地下水交互作用的动态特征，利用校正后的模型模拟长时间尺度下湿地的地表水与地下水水量交互过程，并对研究区进行水均衡分析（图14、表4）。

表4 研究区2005—2019年的年平均水均衡计算表Table 4 Average annual water balance of the study area from 2005 to 2019

图14 研究区2005—2019年水均衡模拟示意图Fig.14 Simulated water budgets of the study area from 2005 to 2019

由图14 可知，研究区域水量输入有降水量、地表水流入、地下水定水头边界入流；输出项包括地表水流出、蒸散发和地下水定水头边界出流。研究区水量主要输入途径是降水和地表水流入，分别占总水量输入的41.12%和52.71%；主要输出途径是蒸散发和地表水流出，分别占总水量输出的50.98%和39.72%。研究区总蓄水量反映水量的盈余和亏缺，经计算总蓄水量为正，代表湿地水储量增加。湿地地表水储量变化不大，表明湿地海绵层蓄水能力较强；但湿地地表水与地下水交换频繁，说明湿地海绵层透水性好、下渗能力强，因此在湿地地表水与地下水耦合模拟中考虑湿地海绵层非常必要。另外，研究区湿地地表水与地下水交互作用显著且复杂，在区域水循环中发挥了一定作用。

3.6 湿地地表水与地下水转化关系

经过率定验证后的湿地地表水与地下水耦合模拟模型可较好地模拟湿地地下水位的动态变化。因此，可用来研究湿地地表水与地下水的转化关系，评价流域内2005—2019年的湿地地下水与河水水量补排关系（图15），图中正值表示湿地地表水对地下水的补给量，负值表示地下水向地表水的排泄量（即基流量）。流域内湿地地下水与河流的补排关系以地下水向河流排泄为主，但在每年6—9月份雨季期间，大量降水使得湿地地表水水位升高，造成湿地地表水对地下水的补给。