基于水文集成模型ParFlow的黑河流域下游地下水-地表水相互作用模拟研究

陆 峥,胡锦华,张 圆,李宗超，杨 晨,杨晓帆*

(1.北京师范大学地理科学学部地表过程与资源生态国家重点实验室，北京 100875；2.北京师范大学地理科学学部自然资源学院，北京 100875；3.中国环境监测总站，北京 100012；4.普林斯顿大学高草甸环境研究所土木与环境工程系，美国 普林斯顿 08544)

地下水和地表水是区域和全球水文循环的基本要素，它们在温度和化学组分上存在着普遍性的差异；且两者具有紧密耦合的水力联系和频繁的转化关系，共同构成了一个完整、复杂、密不可分的系统。因此，定量研究地下水与地表水的相互作用规律对流域水文和生态环境科学具有重要意义。水文集成模型耦合了地下水与地表水模型，能够较为准确地描述自然界中水循环变化的复杂驱动因素，从而描述各种驱动因素对地下水与地表水水文循环的影响机理。然而，由于自然系统存在高度的空间异质性，同时地下水-地表水相互转换的同时还伴随着能量转化和溶质运移，因而定量描述其物理过程及与之耦合的生物地球化学过程具有极大的挑战性。另外，针对不同尺度所采用的不同求解方案也会带来难以量化的不确定性。为了解决上述问题，需要建立和发展基于物理过程的、高时空分辨率的水文集成模型。但是，当前的水文集成模型依旧受到耦合方式、网格和离散化方案、驱动和校验数据质量、计算效率等因素的制约，因此评估地下水-地表水耦合模型在上述几方面的表现尤为重要。通常，在解决具体的流域实际问题前需要使用简单的二维剖面算例对模型进行测试、验证和评价。典型的水文地质剖面数值模拟研究可以定量描述局部地表径流和地下水的侧向及垂向运动，从而获得精确的二维地下水动力场，对于揭示区域地下水与地表水之间的转换关系和水文循环规律特征具有重要意义。在二维剖面数值模拟算例中，常常需将实际问题进行简化、抽象，即建立空间结构更加简单的几何概念模型(如图1所示，将真实流域水文循环概念模型简化为理想化的“V”字形双倾斜坡面流域模型)，并简化、提炼相应的物理过程[如图1(a)中的降水、蒸散发、植被冠层截留、下渗、地下水侧向流等]，最终利用解析解或实验观测对目标模型模拟的结果进行验证和评估。尤其在评估地下水-地表水耦合模型时，需基于简化的水文地质剖面结构，验证其数值解的准确性并量化其不确定性，评价其水文响应模拟结果的时空分布演变特征和统计特征。为此，本文使用基于物理过程的、地下水-地表水紧密耦合的水文集成模型ParFlow，依据不同案例设置，定量描述了地表水与地下水的相互作用，模拟计算了不同情境下地下水-地表水交互的水文响应。

图1 真实流域水文循环的简化建模方案示意图Fig.1 Diagram of simplified model of the hydrologic cycles in a reale-world watershed

1 研究方法

1.1 ParFlow模型

本研究中所采用的水文集成模型ParFlow (https://www.parflow.org/) 是一款开源的、可并行计算的、面向对象的三维地下水-地表水耦合模型。ParFlow模型基于三维地下水流动Richards方程和二维地表水流动Saint-Venant方程的二阶动力波(Kinematic wave)或扩散波(Diffusive wave)近似解，将饱水带、包气带和地表水视为一个完整的水文连续体，以实现全面刻画不同尺度(例如坡面、子流域、流域、国家/区域、大陆和全球尺度等)下的水文循环过程。ParFlow模型的主要特征见表1，更详细的求解处理方案和功能参见相关文献。

表1 ParFlow模型的主要特征Table 1 Summary and specification of ParFlow

图2为ParFlow (Parallel Flow)模型在美国俄克拉荷马州Little Washita流域的一个应用案例。在Little Washita流域的一系列研究中，首先进行了简单的基准测试算例研究，随后循序渐进地扩展到二维剖面的各类研究，最终实现了不同尺度下的流域水文循环三维模拟。秉承上述理念和逻辑递进关系，本研究选择黑河流域下游的典型二维剖面，建立了高精度的地下水-地表水耦合模型，验证了模型的模拟精度，并定量分析了水文响应特征，为后续的流域三维地下水-地表水耦合模拟及其他相关工作奠定坚实的基础。

图2 ParFlow模型在美国Little Washita流域的地下水-地表水交互模拟应用案例示意图Fig.2 Diagram of application of ParFlow to groundwater-surface water simulation case over the Little Washita,USA

ParFlow模型可跨系统适用于Linux/Unix/OSX，包含诸多不同的求解器，其中IMPES求解器适用于单相、完全饱和的情况，RICHARDS求解器适用于地下水变饱和的情况，用户可以自主选择是否耦合地表流动过程，因此ParFlow模型可用于求解饱和带-中间带-包气带-地表水流偏微分方程定解问题。本文主要介绍变饱和水流和地下水-地表水耦合的求解器，分别利用三维Richards方程和二维Saint-Venant方程来求解地下水流动过程和地表水流动过程，地下与地表的耦合主要是通过在地下水与地表水交界处的边界条件来实现，使地下水的三维Richards方程和地表水的二维Saint-Venant方程中的压力水头保持一致，从而保证了地下水-地表水交界面上压力的连续性。这种直接耦合的方法，避免了定义界面导水率(interface conductance)，并且降低了数值求解耦合系统过程中的难度。其中，地下水-地表水交界处的边界条件可根据有无积水及积水深度在第一类边界条件(Dirichlet型)和第二类边界条件(Neumann型)间进行切换。对于变饱和地下水流动方程和地表水流动方程分别采用中心差分格式的有限差分法和标准迎风格式的有限体积法进行空间离散，时间离散方法均为隐式后向欧拉法。此外，利用高效稳定的Newton-Krylov非线性求解器进行矩阵求解，所采用的网格为结构化规则网格，可选正交网格或地形网格。

1.2 控制方程

1.2.1 地下水流动方程

依据质量、动量守恒定律，地下水流动(包括包气带、中间带的壤中流和饱和带基流)水量平衡方程采用的是三维Richards方程，可表示如下：

(

1

)

式中：

h

为压力水头(L)；

t

为时间(T)；

S

为比容(L)；

S

(h

)

为土壤相对饱和度(无量纲)；φ为土壤孔隙度(无量纲)；

z

为垂直坐标(向上为正)；

q

和

q

为源汇项(L/T，具体的介绍见下文)；为达西通量(L/T)，可表示为

(

2

)

其中：为土壤水力传导率(L/T)；

K

为土壤相对渗透率(无量纲)。

土壤相对饱和度与土壤相对渗透率之间的关系可通过van Genuchten模型获得，其计算公式分别如下：

(

3

)

(

4

)

式中：

α

为孔隙进气值参数

(

L

);n

为土壤孔隙分布参数(无量纲)；

S

为土壤相对饱和含水量(无量纲)；

S

为土壤相对残余饱和度(无量纲)。

1.2.2 地表水流动方程

依据质量守恒定律，地表浅水流动水量平衡方程采用的是二维Saint-Venant方程，可表示如下：

(5)

式中:

h

为积水深度(L);为地表流速(L/T)

;q

为与地下水-地表水的交换通量(L/T)

;q

为降水通量(L/T)。

在ParFlow模型中，根据地表水动力波近似，若忽略动量平衡方程中的局部惯性项、对流惯性项和压力差项,则二维Saint-Venant方程可简化为

g(S

-

S

)=

0

(

6

)

式中：

g

为重力加速度(L/T)；

S

为坡面斜率(无量纲)，其中

i

为

x

或

y

轴方向；

S

为摩擦斜率(无量纲)，可表示为

(

7

)

此外，根据Manning-Strickler方程可以建立地表流速与积水深度

h

之间的关系式为

(

8

)

式中：

N

为曼宁粗糙度

(

T/L

)

。

1.2.3 地下水与地表水的耦合方法

地下水与地表水的耦合通过自由表面边界条件来实现。由于需要保证地下水-地表水边界处压力水头和通量连续，则地下水的压力水头

h

和地表水的积水深度

h

均等于地下水-地表水边界处的垂直平均压力

p

，即

:

p=h

=h

(

9

)

地下水的上边界流动通量与地下水-地表水边界处的交换通量相同，即

:

q

=q

(

10

)

根据地下水动量方程

(

2

)

可知，地下水-地表水边界处Neumann型边界条件

q

[

若积水深度大于0，则地下水-地表水边界处由第二类边界条件(Neumann型)可改为第一类边界条件(Dirichlet型)

]

为

(

11

)

地下水-地表水边界处的交换通量为

q

=τ(h

-

h

)

(

12

)

式中：

τ

为边界比例系数。

q

与公式(5)中地下水-地表水的交换通量相同，且可改写为

(

13

)

式中

:

‖

p,

0‖算子结果返回值为

p

和0中的最大值。将上式结果代入公式

(

11

)

中，可得到：

(

14

)

将公式

(

14

)

左侧项重新代回公式

(

1

)

中，可得到：

(

15

)

当

p

小于0时，地下水-地表水边界处为第二类边界条件(Neumann型)，地表水流动边界方程关闭，公式(15)为标准的地下水三维Richards方程；当

p

大于0，也就是存在地表积水时，地下水-地表水边界处由第二类边界条件(Neumann型)转化为第一类边界条件(Dirichlet型)，地表水二维Saint-Venant方程被激活，则利用公式(14)计算。

2 研究区概况与模拟方法

2.1 研究区概况

黑河流域位于97°6′～102°0′E，37°30′～42°42′N间，东西、南北各横跨约5°，发源于祁连山中段，北至中蒙边境，东与石羊河流域接壤，西与疏勒河流域毗邻，总面积为14.3×10km(见图3)。该地区高程范围为869～5544 m，气候主要受中高纬度西风带环流的控制和极地冷气团的影响，气候干燥，降水稀少而集中，多大风，日照充足，太阳辐射强烈，昼夜温差大，年平均气温约-1℃，年均降水量为500 mm。按照气候区划，黑河流域总体上属于半干旱高寒气候。黑河流域上游地处高原亚寒带亚干旱区，中下游为中温带干旱型气候，独特的气候水文条件使之成为开展流域水文过程研究的理想区域。北京师范大学和中国科学院西北生态环境资源研究院于2012年9月起在黑河中游和下游地区共同布设了水文气象观测网，其中额济纳旗核心绿洲二道桥东至七道桥典型河岸林区域为黑河流域下游核心观测区，区内布设了5个水文气象观测站点。

图3 黑河流域下游巴牙吉呼(BJH)地区二维剖面设置Fig.3 Configurations of the BJH case

黑河流域下游的巴牙吉呼(Bajajihu,BJH)地区在行政区划上隶属内蒙古自治区阿拉善盟额济纳旗，其气候类型属于中温带极干旱区，陆表下垫面多为裸地、草地。下游额济纳绿洲是天然的绿洲生态系统，结构简单且极度脆弱，属于极干旱气候区，多年平均降水量不足45 mm，地带性植被多为温带小灌木、半灌木荒漠植被。该地区因降水稀少，植被生长主要依靠地下水，因此获取区域水文循环特征与量化各要素之间的转化关系极为重要。另外，虽然该地区地形较为平坦，但含水层水文地质情况极其复杂，地下水流向(尤其侧向流)尚不明确，因而在此处进行二维剖面的模拟意义深远。同时在前期的研究中，干旱荒漠区地表水与地下水的转化机制及转化过程、浅层地下水的形成机理及运移方式等科学问题尚未有较完善的解决方案。

因此，本次在两个水文气象观测站点之间构建了长为850 m的BJH地区典型二维剖面，旨在验证ParFlow模型预测精度的同时，概算地表水-地下水交互过程及其中的水文响应。西北的胡杨林站海拔为876 m，东南的混合林站海拔为874 m，两站点之间的海拔呈现单调递减的趋势，地表覆盖为裸土并有块状的稀疏胡杨和柽柳。

2.2 模拟方法

2.2.1 二维剖面概述及模拟设置

BJH地区二维剖面的计算区域如图3(a)所示，长为850 m，高为4 m，水平分辨率(即

x

方向)为10 m，垂直分辨率(即

z

方向)为0.01 m。本次模拟总时长设置为24 h(2015年9月3日)，时间步长设为10 min，与降水观测的时间间隔相同；地下水水位原始数据时间间隔为30 min，采用线性内插方法获取时间间隔为10 min的地下水水位值；蒸散发值为日平均观测值，包括裸地蒸发量和植被蒸腾量。ParFlow模型可以输出多类地表水和地下水变量，诸如地表径流、地表水储量、包气带水储量、饱和带水储量、网格中心水位等。取二维剖面两侧(对应胡杨林站表层土壤2 cm和4 cm)的土壤饱和度换算成土壤水分，并与时域反射仪(Time Domain Reflectometry,TDR)观测值进行对比，同时分析了6个典型时刻下BJH地区二维剖面的土壤水分分布和4个不同水平位置处的土壤垂向饱和度分布曲线。

2.2.2 初始条件和边界条件设置

根据BJH地区二维剖面两端站点的观测水位设置初始水位，采用线性内插方法获得剖面内部每个网格的初始水头。模拟区域的底部、前侧和后侧均设置为零通量边界，左侧和右侧采用定水头的边界条件，其值取自剖面两端的胡杨林站和混合林站的观测水位。BJH剖面2015年9月3日的降水和地下水水位时间序列变化曲线，见图4。

图4 黑河流域下游BJH地区剖面2015年9月3日 的降水和地下水水位时间序列变化曲线Fig.4 Time-series of precipitation and water table depths on 3rd Sept.,2015 for the BJH transect

2.2.3 模拟参数设置

孔隙度、饱和导水率、曼宁粗糙度、van Genuchten模型参数均设置为均质，其中孔隙度、饱和导水率和van Genuchten模型参数取自中国土壤数据集，曼宁粗糙度修改自文献[32]，具体参数设置详见表2。

表2 黑河流域下游BJH地区案例中的参数设置Table 2 Parameters in the BJH cases

3 结果与讨论

ParFlow模型模拟的胡杨林站表层土壤水分(2 cm和4 cm)时间序列曲线与TDR观测结果的对比，见图5。

图5 ParFlow模型模拟的胡杨林站表层土壤水分与 TDR观测值的对比图Fig.5 Surface soil moisture time-series comparison between ParFlow simulation and TDR observation at the Populus euphratica site

由图5可见，利用ParFlow模型模拟得到的胡杨林站表层土壤水分模拟值与观测结果基本一致，2 cm和4 cm深度的模拟值与观测值的Pearson相关系数(

R

)分别为0.94和0.93，且均方根误差(RMSD)较低(小于0.008 m/m，即小于整体动态变化范围的10%)；降雨前后表层土壤水分ParFlow模型的模拟值与观测值有较高的一致性，但在降雨期间，表层土壤水分ParFlow模型的模拟时间序列曲线与观测值之间的增速和增幅有所差异：降雨开始后不久，表层土壤水分ParFlow模型的模拟值迅速上升到峰值，而观测值上升速度相对缓慢，推测有两个原因导致了这个现象，即土壤水分的“记忆效应”和胡杨林的冠层截留。总体而言，表层土壤水分的ParFlow模型模拟值与观测值显示出很好的契合度和相关性。

6个典型时刻下黑河流域下游巴牙吉呼(BJH)地区二维剖面的土壤水分分布，见图6。此6个时刻分别为初始状态(00∶00)、降水事件的起始点(06∶00)、降水强度最大时刻(11∶30)、降水事件结束(14∶30)、排水的中间时刻(17∶30)、模拟的最后一个时间步长(24∶00)。

图6 6个典型时刻下黑河流域下游巴牙吉呼(BJH)地区二维剖面的土壤水分分布图Fig.6 Sectional distribution plot of soil moisture in BJH transect at six typical moments in the downstream of the Heihe River Basin

由图6可见，00∶00 与06∶00时刻相比，剖面区域顶部土壤水分差异明显，说明降水开始后表层土壤水分模拟值迅速增加，至06∶00时下渗面与地表平行；随着降水的持续进行，地表积水因重力作用沿斜坡向下渗透，并随着降雨强度的持续增加而下渗加快，剖面顶部逐渐饱和，至11∶30时饱和深度达到1 m以上，此时模拟区域的最下层土壤和地表浅层土壤均呈现饱和状态，中间1～3 m处呈现出一个垂直分布的非饱和过渡带；临近降雨结束，降雨速率下降，导致地表浅层土壤水分的饱和度减小,随着降水事件的结束(即11∶30后)，中间饱和部分的水分继续下渗，但因地形作用地下水开始横向流动，侧向流使得地下水水位不再与地表保持平行；14∶30后地下水水位开始从高海拔处向低海拔处倾斜，由于持续的蒸发作用，表层土壤开始干燥，至17∶30以后整个区域达到相对稳定状态，土壤水分剖面图趋于稳定。这6个典型时刻的土壤水分变化表明，ParFlow模型模拟能够高时空精度地刻画出BJH剖面上土壤水分运移的动力学过程，可较为精准地表征降水、蒸散发驱动下剖面内的水文循环过程。

此外，为了进一步检验ParFlow模型对各时间步长下土壤水分迁移的模拟情况，本次还绘制了4个不同水平位置(

x

=1 m、28 m、56 m、85 m)处土壤垂向饱和度时间变化廓线(见图7)，时间间隔为30 min，与模拟所用的时间步长相同，并用不同颜色表示从00∶00(红色)到24∶00(紫色)时刻的土壤瞬时垂向饱和度分布。从不同水平位置处土壤垂向饱和度分布曲线的变化情况可以看出，土壤垂向饱和度时间变化廓线受降雨速率、蒸发、地下水水位变化和地形的综合影响，在各个因素的共同作用下不同水平位置处的土壤垂向饱和度产生了不同的响应。在降雨事件之前，4个水平位置处的土壤垂向饱和度时间变化廓线对蒸发的响应模式几乎是相同的；降雨事件期间的土壤垂向饱和度时间变化廓线受降雨强度控制，在降雨事件停止后(从17∶30到24∶00时刻，见图7中紫色的廓线)，土壤垂向饱和度时间变化廓线的分布开始稳定；降雨事件后不久，由于蒸发作用，近地表处饱和消退，土壤垂向饱和度时间变化廓线的分布开始趋于稳定。此外，在

x

=28 m和

x

=56 m处的土壤垂向饱和度时间变化廓线在各个时间步长上的表现近似，但在降水停止后的几个时间步长内差异较大，这表明除了降雨事件开始后的一小段时间外，这两个位置处侧向的排泄和补给几乎相同。然而，由于地下水水位的变化，剖面边缘

x

=1 m和

x

=85 m处的土壤垂向饱和度时间变化廓线表现出不同的分布模式，缘于地下水水位变化和地形的综合影响。从以上4组土壤垂向饱和度时间变化廓线可见，ParFlow模型模拟的土壤垂向饱和度分布时间序列具有一定的准确性和敏感性，受地形、蒸散发等驱动因素控制而产生不同的水文响应，可以较好地表征实际水文动力学过程。

图7 不同水平位置处土壤垂向饱和度分布曲线的对比Fig.7 Comparison of vertical saturation distribution curves of soil at different horizontal positions in different time steps

4 结论与展望

本文使用基于物理过程的、地下-地表紧密耦合的水文集成模型ParFlow，依据不同案例设置，定量描述了地表水与地下水的相互作用，模拟计算了不同情境下地下水-地表水交互的水文响应。从简单到复杂的基准测试案例，再到简化的实际问题，ParFlow模型均可较好地解决非线性的水文动力学问题，展现出对驱动因素较为精准的水文响应信息。具体的，在黑河流域下游的巴牙吉呼地区(BJH)剖面展开案例研究，获取了小时内时间步长的水文响应。通过对土壤水分和土壤饱和度模拟结果的验证以及时空分布分析，表明ParFlow模型在模拟地下水-地表水交互方面具备良好的时空准确度。

然而，目前的ParFlow模型仍存在许多缺陷和不确定性，例如水文响应信号依然受到模型驱动数据和空间异质性的影响，同时尺度效应也会影响模拟结果的代表性。此外，ParFlow亦可以耦合陆面过程模型(CLM)和大气模型(WRF或ARPS)，目前均已有较为成熟的版本，但其应用仍有待开发和利用。未来水文模型应与生物地球化学循环、植被动力学模型、溶质运移模型相结合，以更加全面地刻画自然界中真实的水文过程。