清华大学水利系遥感水文课题组在缺资料流域径流估算方面取得新进展

引言

径流是水循环中的重要变量。准确可靠的长序列径流数据，是水资源调查评价的基础。全球径流观测站点分布不均，许多区域尤其是高寒地区缺乏径流观测数据。水文研究中常利用水文模型模拟径流，但通常需要实测流量数据率定模型参数，在无测站、缺资料流域的应用受到限制。为解决上述问题，清华大学水利系遥感水文团队，围绕国际水文科学协会(IAHS)十年计划(2003?2012)：Prediction in Ungauged Basins (PUB)相关科学问题，开展了长期系统性研究。近期，团队李雪莹博士生、龙笛教授及合作者，以土壤水动态变化信息约束水文通量，基于高维参数空间优化的数学模型思想，构建了一种不依赖实测流量的径流求解模型(Soil Moisture to Runoff, SM2R)。模型在青藏高原江河源区20个测站进行了初步验证，表现了较好的精度和稳健性。

SM2R模型由SM2RAIN模型 ^[1-4] 发展而来，通过建立径流与土壤含水率、降水之间的函数约束，克服了SM2RAIN仅能在站点及网格尺度应用且忽略地表径流的局限，将“由状态变量推求通量”的思路，应用于流域径流量估算。此外，本研究针对青藏高原高寒流域，量化了冰川净质量变化、高海拔地区降水产流的滞后时间对土壤包气带等效输入水量的影响，一定程度扩大了SM2R模型的适用范围。研究成果近期以Soil Moisture to Runoff (SM2R): A Data-Driven Model for Runoff Estimation Across Poorly Gauged Asian Water Towers Based on Soil Moisture Dynamics为题，发表于水文水资源领域权威期刊 Water Resources Research 。

研究区与数据

研究方法

SM2R模型的框架如图2所示，关键模块包括：(1)考虑冰川和积雪模块的土壤水量平衡，(2)建立水文通量(蒸散、下渗、及径流)与土壤体积含水率的函数约束。以系统状态变量达到最优为指标，即水量平衡求解的土壤含水率变化和ERA5L估算的土壤含水率变化最接近，利用基于梯度下降的高维参数空间优化方法，求解各参数并估算径流。

3.1  建立土壤水量平衡

在青藏高原，降水相态和冰川质量变化等对径流模拟有重要影响。SM2R模型在计算等效输入水量 P _Eq (式(1))时，将总降水划分为液态降水(降雨)和固态降水(降雪)，冰川积累/消融对土壤包气带等效输入水量的影响为：(1)冰川积累期，固态降水(降雪)首先储存在冰川上，其余形成土壤层的有效输入(式(4))；(2)冰川消融期，冰川融化为液态水，同液态降水(降雨)共同作为土壤层的有效输入(式(5))。

式中， dh _g 表示冰川高程变化(mm/mon)； snowfall _g  和 rainfall _g 分别代表考虑冰川质量变化后的等效降雪和等效降雨； dh 表示由冰川净质量变化引起的流域平均水深变化，由式(6)计算：

式中， Area _g 指流域内冰川面积(km ² )；Area指流域集水面积(km ² )； ρ _g 指冰川密度(取850 kg/m ³ )，ρw指水密度(取1000 kg/m ³ )。

此外，高海拔地区(例如喜马拉雅山脉)的实测径流可能滞后于降水数月，且径流与降水出现峰值的时间不同，这是由于冬春季的固态降水通常在冬春季冻结，当温度升高后再融化成液态水参与到水循环中。本研究将具有“储存冬春季降雪并补偿夏季降雨”的流域定义为“降雪主导”流域，提出积雪指数(Snow Index ( SI )，式(7))表征流域受固态降水和积雪的影响。

式中， r _{snow / prep} 指多年平均降雪量占总降水量的比例； r _{persistent snow cover} 表示流域多年平均积雪覆盖比例。以 SI > 0.16的流域被认为是“降雪主导”流域，该比例阈值参考了聚类方法中常用的head/tail breaks原则。

针对“降雪主导”流域，月尺度液态降水( rainfall _g , 式(5))全部参与流域水循环，即补给土壤水并形成输出通量。月尺度固态降水( snowfall _g , 式(4))参与水循环的实际水量，由该月温度决定，假设每月固态水可融化为液态水的实际水量，与当月气温减去温度阈值的差值成正比：

式中， snowfall _a 代表每月可融化的实际水量(mm/mon)；k表示温度?融化量的比例因子； T ( t )指月平均气温(°C)， T _threshold 指温度阈值，设定为研究时段(1981?2020)的多年月平均气温减去1倍标准差。由于在年尺度，总的实际融化水量(Snowfalla, 式(8))应该与实际固态水量( snowfall _g , 式(4))相等，比例因子 k 可由式(9)计算：

在考虑冰川和积雪的共同影响后，“降雪主导”流域的等效输入水量由式(10)计算，非“降雪主导”流域的等效输入水量由式(11)计算。

3.2  确定各变量的函数约束

SM2R模型中建立各输出通量(蒸散、下渗和径流)与土壤含水率之间的数学函数约束。实际蒸散发( ET )由双曲正切函数约束(式(12))，表示当土壤含水率为零时，实际蒸散发为零，此时，实际蒸散发受水分控制；当土壤接近饱和含水量时，实际蒸散发接近潜在蒸散发( PET )，函数值逼近渐近线值域，此时，实际蒸散发受能量控制。式(12)中参数 a 表示蒸散发由水分控制向能量控制的过渡，可由sigmoid函数(式(13))拟合，sigmoid函数的值域为[0, 1]。

式(12)及后续式(15)?(17)中的变量 φ 表示土壤孔隙度，由式(14)计算，其中 F _sand and F _clay 分别指砂土和黏土的百分比。

下渗 D 表示土壤水垂向渗漏补给地下水的部分，可表达为土壤体积含水率的幂函数(式15)。式(15)中参数 b 指下渗量的最大值，参数 c 表示水量流失的速率，与土壤包气带的饱和度有关。

在SM2R模型中，快速流Rfast的数学函数形式可表达为等效输入水量和土壤含水率的幂函数(式(16))，参数 e 可表示降水产流的速率。慢速流 R _slow 可表示为土壤含水率的幂函数(式(17))，f可表示慢速流的最大值， g 表示慢速流产流速率。

3.3  参数优化和求解

在SM2R模型中，各变量由7维参数向量空间确定，即 X = [ a , b , c , e , f , g , z ]。本研究通过土壤性质等数据进行参数初始化，以提高优化收敛效率和参数的物理约束，但各参数最终由优化过程确定。由于参数 a (确定实际蒸散发)已被sigmoid函数约束，因此，对初始值无限制。参数 b ， c 用于确定下渗量，与土壤的水力学特性和下渗能力有关，Jahlvand 等 ^[8] 和Brocca等 ^[9] 基于实测土壤资料，根据土壤类别(可由HWSD数据集获取)提供参数 b 和 c 的推荐值。

       参数 e 表示降水产流的速率，由流域多年平均径流系数进行初始化。参数  f  和 g 用于确定慢速流，其中，参数f表示慢速流的最大值，与土壤持水能力和含水层给水度有关系，以式(18)对参数  f  进行初始化：

式中， AWC (mm/m)表示土壤包气带的持水能力，可基于HWSD数据集中针对不同土壤类别的推荐值确定； SY 表示含水层给水度，指土壤包气带可排出水量占土壤体积的百分比(%)； z’ 指慢速流产流的土壤深度，设置为100 cm，这与HWSD数据集估算土壤持水能力( AWC )时设置的参考深度相同。针对更深层土壤(>100 cm)，假设土壤水流失的主要过程为下渗，而非产生慢速流。

参数 g 确定慢速流的幂函数形式，以多年平均壤中流和多年平均总径流的比值，对其进行初始化(以ERA5L估算的壤中流和总径流设定初始值)。最后，参数z表示满足水量平衡闭合(式(1))的土壤层深度，利用ERA5L中的土层深度(2890 mm)对z进行初始化。

主要研究结果

本研究利用SM2R模型在青藏高原20个测试流域模拟1981?2020年月尺度径流，对比实测径流，从相关性( r )、归一化均方根误差( NRMSE )、纳什效率系数( NSE )、对数纳什效率系数( logNSE )、Kling-Gupta效率系数( KGE )等五个统计指标，评估SM2R模拟径流的表现(图3)。结果表明：20个站点的相关系数 r 均高于0.74，归一化均方根误差 NRMSE 均小于0.22；16个站点综合性评估指标表现良好( NSE > 0.56， logNSE > 0.43， KGE > 0.46)，表明在不依赖实测径流的情况下，SM2R模型在多数测试流域的模拟结果较为可靠。

尽管SM2R模型以ERA5L为驱动数据，但除了在恒河的Kali Khola站点和雅江的奴各沙站点，二者估算径流的表现相近外，其余18个站点ERA5L径流的评估指标整体低于SM2R模拟径流(图3)，特别是在冰川积雪广泛分布的印度河流域(SM2R模拟径流的NSE较ERA5L再分析径流提升134%?220%)和冻土广泛分布的长江源和黄河源(SM2R模拟径流的NSE较ERA5L提升101%?118%)。

图3 SM2R和ERA5L估算径流与实测径流对比评估指标雷达图。评估指标NRMSE在雷达图中展示为1-NRMSE，因此所有指标的理想值均为1。当评估指标<0时不在雷达图中显示，只标记具体值(例如ERA5L估算径流在Shatial Bridge、Skardu Kachura、Gongshan、Zhimenda、Tangnaihai及Maqu的评估结果)

对比1981?2018年SM2R模拟和Gou等 ^[7] 基于VIC模型研制的CNRD数据集发现(图5)，二者在青藏高原中国境内的13个站点的径流模拟效果接近。在SM2R模型评估指标较低的贡山、直门达、唐乃亥和玛曲站，CNRD数据集的评估指标较SM2R模拟结果更低，反映在土壤冻融过程复杂的江河源区，降水径流关系的鲁棒性较低，径流模拟的不确定性较大。ISIMIP可提供径流分析的常用模型，但其模拟效果在青藏高原表现较差。评估1981?2005年16种ISIMIP模式平均在青藏高原输出的径流结果发现，大多数站点(13个)的综合性评价指标 NSE < 0。

图5 SM2R和CNRD估算径流与实测径流对比的评估指标雷达图。评估指标NRMSE在雷达图中展示为1-NRMSE，因此所有指标的理想值均为1。当评估指标<0时不在雷达图中显示，只标记具体值(CNRD估算径流在直门达、唐乃亥及玛曲站的评估结果)。

讨论

本研究构建的SM2R模型，实现了利用土壤水动态变化信息，对流域径流较为可靠的模拟，其优势主要体现在：(1)在无测站、缺资料流域应用潜力较大。(2)对比水文模型或陆面模式，SM2R模型的输入数据(包括ERA5L获取的降水、潜在蒸散、土壤体积含水率三类)和模型参数(共7个)较少，降低了输入数据不一致性带来的不确定性和“异参同效”的风险。(3)对比Koster等 ^[10] 利用线性回归拟合土壤含水率和径流相关关系的研究，SM2R模型建立了土壤含水率、降水和径流之间的幂指数函数约束，有效提高了径流模拟效果。Koster等 ^[10] 利用实测径流率定线性回归参数，但在美国145个测试流域中，仅有22个站点对比实测径流的相关系数高于0.7，而本研究中所有流域的相关系数均高于0.74。

SM2R模型的局限主要体现在两方面：(1)由于SM2R模拟径流基于径流、土壤含水率和降水间的幂指数函数约束，这种函数关系在日尺度或小时尺度的鲁棒性较差，因此，SM2R模型更适用于长时序月径流分析。(2)SM2R模型在土壤冻融过程复杂的长江源和黄河源表现不佳。一方面，目前土壤水再分析和遥感产品在长江源/黄河源等冻土广泛分布地区的误差较大；另一方面，由于冻土退化后土壤孔隙度变大，直接影响各水文通量计算。因此，获取更精确的土壤水数据和土壤性质数据，有望进一步提升SM2R的模拟精度。

SM2R模型的鲁棒性可从以下三个方面分析：(1)模型优化的目标变量是土壤含水率变化( dθ / dt )。图6表示ERA5L和陆面模式(GLDAS NOAH和CLSM模式平均)在典型流域估算的土壤含水率和土壤含水率变化，研究发现：尽管不同数据源估算的土壤含水率(θ)的绝对值差异较大，但土壤含水率的变化(dθ/dt)一致性较高，这是SM2R模型优化参数和准确模拟径流的关键。(2)SM2R模型基于可靠的土壤含水率变化约束降水产流过程，通过参数优化和调整，降低降水数据不确定性对径流的影响，最终实现较可靠的径流估算。(3)SM2R模型中各水文通量独立计算，即每种输入数据的不确定性，仅通过相关的数学函数影响特定分量和参数，可减小误差传递，提高模型鲁棒性。

总结

本研究基于高维参数空间优化的数学思想，构建了以土壤水动态变化约束径流求解的新方法(SM2R)，实现了在不依赖实测径流资料条件下，较为准确地估算了青藏高原不同气候和地理特征流域的长时序径流。SM2R模型的输入数据和模型参数较少，径流模拟效果显著优于ERA5L和ISIMIP的径流估算，与利用实测径流率定参数的VIC模型表现接近。本研究提出的SM2R模型提供了“由状态量推求通量”的径流模拟新思路，有潜力在其他无测站、缺资料流域的径流估算和分析中提供参考借鉴。

本研究受到第二次青藏高原综合科学考察(2019QZKK0105)、基金委“西南径流重大研究计划”(92047301、92047203)等项目支持，谨致谢忱。