林芝地区复杂水文环境径流模拟研究

Study on the Runoff Simulation in Nyingchi District Under Complex Hydrological Environment

摘要

西藏林芝地区水文环境复杂，区域水循环受降水、冰川融水、闸坝调蓄等多因素共同影响。同时，近年来在气候变暖背景下，该流域洪涝灾害风险加剧。文章提出从简单水文环境的径流模拟出发，进一步完善了冰川融水、闸坝调蓄等复杂水文环境下的径流模拟方法。研究表明，以栅格模式的Navier-Stokes方程为基础，结合冰川融雪、闸坝调蓄构建的复杂水文模型，能够较好地模拟复杂水文环境下的林芝地区流域径流。研究成果可为林芝地区流域水文管理和复杂水文环境下径流模拟研究提供理论和实践参考。

Abstract

The hydrological environment of the Nyingchi District in Xizang is complex, with the regional water cycle influenced by a combination of factors including precipitation, glacial meltwater and reservoir regulation. At the same time, against the backdrop of global warming in recent years, the risk of flood disasters in river basins has been intensified. Starting from runoff simulation in simple hydrological environments,this study further refines the methods for simulating runoff under complex hydrological conditions involving glacial meltwater and reservoir regulation. The study indicates that a complex hydrological model, based on grid-based Navier-Stokes equations and incorporating mechanisms for glacial meltwater and reservoir regulation, is capable of effectively simulating runoff in the Nyingchi region's basin under complex hydrological conditions. Our findings provide important theoretical and practical references for hydrological management in this region and for research into runoff simulation under complex hydrological environments.

Graphical abstract

关键词

林芝地区 / 复杂水文 / 径流模拟 / 冰川融雪 / 闸坝调蓄

Key words

Nyingchi region / complex hydrology / runoff simulation / glacial meltwater / reservoir regulation

引 言

径流模型是一种用于模拟和描述地表水在自然环境中运动、分布和变化的数学模型^[1]，通过对径流汇流过程进行数学表达，能够分析和预测水文环境的径流行为^[2]，是流域水文特征分析与洪水预警预报的重要手段。随着计算机技术的不断发展，水文学家逐渐采用基于水文机理的物理模型，以更精确地模拟径流汇流过程^[3]。物理模型基于物理学的质量、动量守恒定律以及径流汇流原理构造水动力学方程组^[4,5]，使得模型能够更全面地考虑降水径流在时空上的变化，同时能够捕捉径流汇流的动力学机制和相邻栅格间的空间关系，在径流模拟中得到广泛应用。当前，国际上主要的物理模型有TOP MODEL^[6]、SWAT^[3]、Green-Ampt^[4]等，我国也相继发展出多种分布式的水文模型，如新安江^[5]、双超产流^[7]、河北雨洪^[8]和姜湾径流^[9]等。

西藏林芝地区水文环境异常复杂，区域水循环受冰川融雪以及闸坝调蓄等诸多因素影响，径流模拟难度较大。在冰川融雪方面，Singh等^[10]采用SWAT模型研究了喜马拉雅陡峭山脉亚流域尺度上的积雪和融雪；班春广等^[11]以雅鲁藏布江帕隆藏布上游流域为例，研究了高寒流域水文模拟与径流水源解析。已有研究通常采用传统的能量平衡模型进行冰川消融量的估算^[12]，但此类方法计算复杂且需要大量的气象和冰川监测数据，对于地形复杂、数据稀缺的偏远地区来说，推广应用存在一定难度。度日模型（Degree-Day Model）是一种简化的冰川消融模拟方法，其基于气温数据估算冰川消融量^[13,14]，由于数据资料的可获得性，度日模型在冰雪消融估算中的应用更为广泛。在闸坝调蓄方面，Bulti等^[15]在东非阿瓦什河流域利用SWAT模型对流域上游的水循环进行了模拟，其研究也指出，由于闸坝运行对水文参数影响较大，SWAT模型难以精准模拟流域水循环过程^[16]。

针对西藏林芝地区复杂的水文环境，如何精确模拟区域水文循环过程，已成为提升区域水循环机理认识的关键科学问题。已有研究已建立基于Navier-Stokes方程的径流汇流模拟模型，在林芝地区泥石流灾害预警中取得了良好的应用效果^[17,18]。本文以该模型为基础，进一步完善了冰川消融和闸坝调蓄模拟环节，发现完善后的模型更适合于林芝地区复杂水文环境下的径流模拟研究。本研究可为复杂水文环境下的径流模拟研究提供新方法，也为林芝地区水文管理和洪峰预警提供基础数据。

1 研究区概况与试验方案

1.1 研究区概况

林芝地处西藏东南部、雅鲁藏布江中下游，平均海拔3 100 m，地形呈西北高、东南低的走势。受孟加拉湾暖湿气流的影响，气候温暖湿润，夏季雨量充沛。区内广泛分布冰川与水库，其水文环境非常复杂。作为林芝地区的典型代表流域，尼洋河是雅鲁藏布江北侧的最大支流，河流全长约307.5 km，天然落差达2 080 m，平均纵坡为7.25‰，流域面积17 535 km²，日平均每秒流量为538 m³/s，年径流量为220亿m³。流域内冰川和常年积雪面积为6 728.2 km²，占土地总面积的38.45%，规模较大的冰川约有40处。此外，流域内有多布电站和老虎嘴电站等闸坝调蓄。尼洋河发源于工布江达县，在巴宜区附近汇入雅鲁藏布江，干流自汇入口以下流经米林县，沿途支流众多。其中，巴河是尼洋河最大且最为著名的支流，全长约89 km，流域面积4 229 km²，平均日流量为184 m³/s，海拔落差为768 m。本研究采用的主要数据及其来源见表1。

表1　研究数据及其来源

1.2 径流模拟方案

选择工布江达县和米林县两个研究区用于简单水文环境下的径流模型，以验证模型的适用性。在此基础上，选择巴河区，在模型中增加闸坝调蓄，验证模型的适应性；最后选择巴宜区，综合考虑冰川消融、闸坝调蓄，验证模型的适应性。各流域分区如图1所示，水文环境特征如表2所示。

2 模拟方法

2.1 简单水文环境径流模拟

借鉴陈军等^[17]提出的栅格模式下的降雨径流汇流模型，综合考虑了地形、降雨、气温、水位等因素对径流过程中降雨、蒸发、下渗、扩散等环节的影响^[19,20]，基于纳维-斯托克斯方程（N-S方程），在栅格框架下将每个栅格单元设定为独立的水文计算单元，统一描述坡面汇流与河道汇流为栅格尺度上的水流运动^[21]。每个栅格上的水体则被假定为一个可由水深表征的立方体水柱。在时间维度上，将完整的径流汇流过程被分解为若干相等的时间片，每一个时间片的各栅格单元负责完成平流计算和水动力模拟。对于时间片t_{i+\mathrm{1}}，栅格单元c的水量平衡方程可表示为：

W_c^{i+\mathrm{1}}=W_c^i+R_c^i+I{n}_c^i-I_c^i-E_c^i-Ou{t}_c^i（1）

式中，W_c^i和W_c^{i+\mathrm{1}}分别表示时间片i和i+1上栅格单元c的水深，单位m；R_c^i、I{n}_c^i、I_c^i、E_c^i、Ou{t}_c^i分别表示栅格单元上的降水量、流入量、下渗量、蒸发量和流出量，上述参数均换算为栅格水深，单位m。其中，下渗量I_c^i和蒸发量E_c^i依据下垫面性质设定下渗率和蒸发率，并与W_c^i相乘得到；I{n}_c^i和Ou{t}_c^i则通过平流计算进行模拟。

采用多流向法模拟平流项。将中心栅格c的8邻域像元设为b，其水流速度矢量为(v_b^x,v_b^y)，则流入中心栅格c的水量比例S_{b\to c}的计算公式如下：

S_{b\to c}=\mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{1}-\left|\mathrm{\Delta }x+v_b^x/v_{\mathrm{m}\mathrm{a}\mathrm{x}}\right|,\mathrm{0})\times \mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{1}-\left|\mathrm{\Delta }y+v_b^y/v_{\mathrm{m}\mathrm{a}\mathrm{x}}\right|,\mathrm{0})（2）

式中，(\mathrm{\Delta }x、\mathrm{\Delta }y)为邻域栅格坐标相对中心栅格坐标的偏移，单位为像素；v_{\mathrm{m}\mathrm{a}\mathrm{x}}表示时间片上允许的最大速度标量。

将式（2）中的b设置为c，此时\mathrm{\Delta }x和\mathrm{\Delta }y均为0，则计算结果为中心栅格剩余的水量。

2.2 冰雪融水模拟

在林芝地区冰川消融对水资源的贡献显著，因此在径流模拟中，有必要将冰川消融作为水分收入源加以考虑。本文根据度日模型原理，设定冰川消融与气温之间存在线性关系，即冰川消融量与超过某一阈值（通常为0°C）的累积正气温成正比。度日模型数学表达式如下^[22]：

M=\left\{\begin{array}{l}MF(T-T_{\mathrm{0}})\begin{array}{cc}& T>T_{\mathrm{0}}\end{array}\\ \begin{array}{cc}& \end{array}\mathrm{0}\begin{array}{cc}\begin{array}{cc}& \end{array}& \begin{array}{cc}& T\le T_{\mathrm{0}}\end{array}\end{array}\end{array}\right.（3）

式中，M为日消融量（mm），T为日平均气温（℃），T₀为冰川消融的温度阈值（℃），MF为度日因子（mm·℃^-1·d^-1）。鉴于消融量与气温密切相关，以及林芝地区复杂地形导致的气温空间异质性，本研究为确保模拟的准确性，首先基于数字高程模拟（DEM）与气温递减率，对空间插值后的气温进行修正。随后，将校正后的气温数据插值后代入公式（3），计算冰川所在栅格像元的日消融量。最后，将栅格单元的降水量与计算得到的冰川消融量相加，共同构成模型的水量输入。增加冰川消融量后模型的水量平衡公式如下^[23]：

W_c^{\mathrm{i}+\mathrm{1}}=W_c^i+R_c^i+I{n}_c^i-I_c^i-E_c^n-Ou{t}_c^i+W_c^i（4）

式中，W_{\mathrm{c}}^{\mathrm{i}}表示中心栅格时间片上的冰川消融量。

2.3 闸坝水文环境径流模拟

在林芝地区，部分河道建有水电站，水电站的闸坝对于径流具有蓄水与节流的调节作用。因此，为准确模拟受闸坝调蓄影响的径流过程，需要构建相应的径流模拟模型。由于闸坝的调蓄作用会阻碍河道水流向下游流动，形成截流效应。本研究为模拟闸坝调蓄，首先在闸坝点的河流两岸人工勾绘出闸坝线，然后将闸坝沿线栅格分为三类，即坝前栅格、坝上栅格和坝后栅格，如图2所示。其中，坝前栅格位于水库库区内，坝后栅格位于水库下游，两者均与坝上栅格具有8邻域邻接关系。通过分析坝上栅格和水流方向，即可有效识别并区分坝上栅格8邻域内的坝前栅格与坝后栅格。

闸坝调蓄模拟的关键在于如何实现坝前栅格的下泄水流向坝后栅格进行有效迁移。为此，首先将水库的下泄流量换算至时间片上的下泄水量Q，并建立其与泄洪水位变化的关系，计算公式如下：

Q={\displaystyle \sum _{i=\mathrm{1}}^m}Ou{t}_i{C}^{\mathrm{2}}（5）

式中，m为坝前栅格总数；Ou{t}_i为坝前栅格i的下泄流量折算的水位降低量（m）；C为地形栅格的分辨率（m）。若在式（5）中的Q为已知量，则Ou{t}_i为待求量。假设在单个时间片内各坝前栅格的水位在泄洪前后变化相同，则可采用二分法逼近求解。具体计算过程如下：

（1）将坝前栅格的最低高程和最高高程作为二分法求解的上下界，令初始二分法求解的水位下界{\mathrm{L}}_f为最低高程、上界L_t为最高高程；

（2）设当前水位L_{cur}=(L_f+L_t)/2，根据坝前栅格的地形高度和L_{cur}计算出Ou{t}_i，公式如下:

Ou{t}_i=\left\{\begin{array}{l}{L}_{cur}-H_i{L}_{cur}-H_i>\mathrm{0}\\ \mathrm{0}{L}_{cur}-H_i\le \mathrm{0}\end{array}\right.（6）

（3）将式（6）代入式（5），求出L_{cur}水位对应的Q_{cur}；

（4）若Q_{cur}与Q在一定的误差范围内，则算法结束；

（5）若Q_{cur}<Q，则令L_t=L_{cur}，转（2）；

（6）若Q_{cur}>Q，则令L_f=L_{cur}，转（2）。

在时间片的迭代过程中，首先，依据水流方向将坝前栅格的Ou{t}_i迁移至坝上栅格，并记录迁移时的水速方向和大小；之后，在下一个时间片迭代时，将水库的水量依据水速方向一次性迁移至坝后栅格，即完成了水库部分溃坝的模拟。

2.4 复杂水文环境下的径流模拟

首先，针对工布江达和米林研究区，应用简单径流模型进行降水、蒸发、扩散、下渗过程的模拟。同时，采用基于度日模型的融雪径流模型开展冰川融水模拟，并利用集成了闸坝调蓄功能的闸坝径流模型对部分建有水电站的河道进行模拟。模型构建思路见表3。

3 模型验证与分析

采集林芝地区60 m分辨率的地形数据、水库水文数据，以及对应日期的气温、降水、雪水当量等数据，开展林芝地区复杂水文环境径流模拟研究。

3.1 模拟验证方法

由于本文模型模拟结果为栅格单元上的平均水深（式1），而实测数据为河道断面水位，二者无法直接进行对比。此外，本模型为基于栅格进行的径流模拟，其核心思想是将每一个栅格上的水体假定为一个立方体水柱，即水体均匀“平铺”于每一个栅格单元上，这与实际地面差距甚大，导致模型模拟水深与实际水深难以进行精确对照。鉴于此，本文采用归一化的模拟水位曲线与归一化的真实水位曲线的相关系数作为模型精度验证指标。该相关系数越大，表明模型对水位随时间变化的模拟效果越好，越接近实际情况。

归一化模拟水位曲线的计算步骤如下：首先，获取模拟时段内模拟水深的最大值和最小值；其次，将每一时间片上的模拟水深减去最小值，再除以最大值与最小值的差，得到该时间片上的归一化模拟水位；最后，将各时间片上的归一化模拟水位按时间先后顺序连接，即可得到归一化模拟水位曲线。

归一化真实水位曲线的计算采用相同方法。即在模拟时段内，求出模拟水位的最大值与最小值，并将每一个时间片上的模拟水位减去最小值，再除以最大值与最小值的差，即可得到该时间片上的归一化真实水位。

3.2 简单水文环境模拟适应性分析

本研究将构建的简单水文环境径流模型，分别应用于林芝尼洋河流域上、下游的工布江达和米林两个研究区，并基于小时降水和气温数据进行了模拟。通过调整加速系数、扩散系数等参数，使用控制变量法、最小二乘法拟合等方法进行多次试验，率定模型参数，得到了较好的模拟水深效果（图3），线性拟合的相关系数R²分别达到0.74和0.94，充分验证了模型在复杂地形条件下的模拟能力。

分析发现，率定参数的简单水文环境径流模拟模型的归一化模拟水位与归一化观测水位在整体趋势上保持一致。模拟水深的波动幅度应与观测水位的波动幅度相匹配。模型模拟的峰值和谷值出现的时间点应与观测数据中的时间点相吻合。试验证明，该模型对当地水文环境有一定的适应性，模型在模拟水深变化方面表现良好，具有一定的解释力和相关性。

为了更好地验证该模型是否能够持续保持高度适应性，利用上述率定模型对工布江达和米林试验区的水深模拟适应性进行验证。从图4的模拟结果来看，模拟水深变化的大体趋势与真实水位相一致，模拟水深和真实水位之间存在较强的正线性相关性，线性拟合相关系数R²分别为0.82和0.86。从试验结果来看，工布江达试验区模拟效果明显优于米林试验区，如图4所示。这表明，即使在同一地区的相同水文参数，模型在不同时间段的表现也存在差异。这种差异可能源于部分水文要素在模型构建时未能被充分考虑，比如冰川融雪、地下水补给等影响。

3.3 闸坝水文环境模拟适应性分析

3.3.1 巴河试验区适应性分析

闸坝径流模型可以包容有融雪水源的情况。以巴河试验区为例，在闸坝水文环境下进行径流模拟试验，需为巴河试验区的老虎嘴水电站设置闸坝线，并需要将该时间段内的闸坝调蓄数据输入至模型中。参数率定试验得到适用于巴河试验区的模型参数，具体见表4。

为有效评估闸坝调蓄对径流模拟的影响，对比分析了未设置闸坝线、未加载闸坝调蓄流量数据的模型水深模拟结果与设置老虎嘴闸坝线、模型导入老虎嘴调蓄流量数据的模型水深模拟结果，具体如图5所示。结果显示，融雪径流汇流模型对水深的趋势模拟较差，模拟水深的变化相对平缓，且完全没有体现闸坝调蓄的影响。相比之下，闸坝径流模型模拟水深与观测水深的起伏较为接近，在某些时间点的模拟水深与观测水深较为吻合，可以体现闸坝调蓄的影响，线性拟合相关系数R²分别为-0.10（融雪径流汇流模型）和0.82（闸坝径流模型），具体如图6所示。

闸坝径流模型模拟结果的皮尔逊相关系数（0.82）远高于融雪径流模型的模拟结果（-0.10），这表明闸坝径流模型在模拟水深方面与实际观测值的相关性更强。闸坝径流模型模拟结果的决定系数（0.66）也显著高于融雪径流汇流模型（0.01），说明了闸坝径流模型在解释水位变异性方面更为有效，模型拟合得更好。通过对两个不同模型模拟结果的综合对比，并结合前文对融雪模型径流模型的分析，可以明确闸坝调蓄是影响径流模拟精度的关键因素，尤其在试验区水文站与水电站距离较近的情形下，对径流过程的扰动作用进一步凸显。

3.3.2 巴宜试验区适应性验证

采用闸坝径流模型在巴宜试验区闸坝水文环境下进行适应性验证。参数率定试验得到了适用于巴宜试验区的模型参数，具体见表5。

巴宜试验区案例模拟结果如图7所示。从图7来看，该区的适应性验证结果较优秀，相关系数在0.9以上，但是局部水位变化的误差较大，模型作为简化的物理模型，仍然未考虑水文要素，如地下水补给、极端气候事件或特殊地理条件等。在未来的研究中，可以考虑量化更多的水文要素以进一步提高模型模拟能力。巴河试验区和巴宜试验区模型模拟效果优良，证明该径流模型可以相对准确地模拟林芝地区闸坝水位环境试验区的径流过程。

4 结 论

本研究通过在径流汇流模型中增加冰川消融和闸坝调蓄模块，系统分析了不同水文要素如降雨、气温、地形、冰川融雪和闸坝调蓄对径流的影响，并评估了模型在不同流域与水文条件下的适应性和模拟能力。结果表明，经精确的参数校准和适应性验证，模型成功构建了适用于林芝地区复杂水文环境的径流模型，为林芝流域提供可靠的径流水位模拟。

（1）通过对林芝简单水文环境试验区模拟适应性进行分析，验证了简单径流模型在林芝地区简单水文环境下的适应性和稳定性，构建的径流模型为研究林芝地区更为复杂的水文环境模拟提供了基础。

（2）基于融雪径流模型，针对闸坝水文环境的特点，创新性地在模型中增加了闸坝调蓄模拟方法，将其迭代成可以模拟林芝流域闸坝调蓄的闸坝径流模型，验证了闸坝径流模型在林芝地区流域受到冰川融雪和闸坝调蓄双重影响下仍能保持良好的稳定性与适应性。

参考文献

[1] 刘化芳,张成福,冯霜,等.水文模型研究进展[J].绿色科技,2023,25(18):49-53+61.

[2] 封志明,杨艳昭,游珍.中国人口分布的水资源限制性与限制度研究[J].自然资源学报,2014,29(10):1637-1648.

[3] Sabzipour B, Arsenault R, Troin M, et al. Comparing a Long Short-Term Memory (LSTM) neural network with a physically-based hydrological model for streamflow forecasting over a Canadian catchment[J]. Journal of Hydrology,2023,627:130380.

[4] Soriano J M A, Reed M. Spatial aggregation effects on the performance of machine learning metamodels for predicting transit time to baseflow[J]. Environmental Research Communications, 2023,5(11):115002.

[5] Liang W T, Chen Y N, Fang G H, et al. Machine learning method is an alternative for the hydrological model in an alpine catchment in the Tianshan region, Central Asia[J]. Journal of Hydrology: Regional Studies, 2023,49:101492.

[6] 王也,龚定,李致家,等.基于不同河道汇流方法的新安江模型研究[J].水力发电,2024,50(1):8-14.

[7] Xi Y, Peng S S, Ducharne A, et al. Author correction: Gridded maps of wetlands dynamics over mid-low latitudes for 1980–2020 based on TOPMODEL[J]. Scientific Data, 2022,9(1):463.

[8] Qiu A N, Zhang Y J, Wang G X, et al. A modified TOPMODEL introducing the bedrock surface topographic index in Huangbengliu watershed, China[J]. Journal of Mountain Science,2022,19(12):3517-3532.

[9] Bahmani F, Fattahi M H, Sabzevari T, et al. Spatial–temporal analysis of landslides in complex hillslopes of catchments using Dynamic Topmodel[J]. Acta Geophysica, 2022,70(3):1417-1432.

[10] Singh V, Muñoz-Arriola F. Improvements in sub-catchment fractional snowpack and snowmelt parameterizations and hydrologic modeling for climate change assessments in the Western Himalayas[J]. Hydrology, 2021,8(4):179.

[11] 班春广,徐宗学,左德鹏,等.高寒流域水文模拟与径流水源解析——以雅鲁藏布江帕隆藏布上游流域为例[J].北京师范大学学报(自然科学版),2023,59(1):85-93.

[12] Pradhananga D, Pomeroy J W. Diagnosing changes in glacier hydrology from physical principles using a hydrological model with snow redistribution, sublimation, firnification and energy balance ablation algorithms[J].Journal of Hydrology, 2022,608: 127545.

[13] Ju Q, Liu X, Zhang D, et al.Application of distributed Xin'anjiang model of melting ice and snow in Bahe River basin[J].Journal of Hydrology: Regional Studies, 2024, 51(c):101638.

[14] Wang Y C, Li X S, Han Z Y. Determining the cooling magnitude in the east of China during the Last Glacial Maximum using a degree-day model[J]. Elsevier Besloten Vennootschap, 2023,664:13-19.

[15] Bulti A T. The influence of Dam construction on the catchment hydrologic behavior and its effects on a discharge forecast in hydrological models[J].Water Resources Management, 2021,35(6):2023-2037.

[16] 于峰,史正涛,李滨勇,等.SWAT模型及其应用研究[J].水科学与工程技术,2008(5):4-9.

[17] 陈军,李婷,禹政阳.一种基于通用计算的流域水质实时模拟方法[J].中国农村水利水电,2021(2):24-29+36.

[18] 陈宫燕,李婷,陈军,等.基于栅格径流汇流模拟的西藏林芝市泥石流灾害预警模型初探[J].中国地质灾害与防治学报,2023,34(1):110-120.

[19] 骆卓然,刘家宏,张尚弘,等.气候变化对青海湖水位的影响及机理分析[J].水力发电学报,2024,43(3):1-13.

[20] 曾方明,董博,张西营.青藏高原东北部60多年来的气候暖湿化：来自德令哈水文气候的记录[J].盐湖研究,2024,32(4):32-39.

[21] 冯熠.不可压流体力学方程弱解的能量等式[D].扬州:扬州大学,2025.

[22] 尹振良,冯起,刘时银,等.水文模型在估算冰川径流研究中的应用现状[J].冰川冻土,2016,38(1):248-258.

[23] 陈军,刘意,李婷.一种基于稳定流场的流域径流汇流模拟方法[J].中国农村水利水电,2020(4):76-81.

基金资助

西藏自治区自然科学基金项目(XZ202101ZR0042G)

西藏自治区自然科学基金项目(XZ202101ZR0042G)