1 引言 径流是水循环中的重要变量。准确可靠的长序列径流数据，是水资源调查评价的基础。全球径流观测站点分布不均，许多区域尤其是高寒地区缺乏径流观测数据。水文研究中常利用水文模型模拟径流，但通常需要实测流量数据率定模型参数，在无测站、缺资料流域的应用受到限制。为解决上述问题，清华大学水利系遥感水文团队，围绕国际水文科学协会(IAHS)十年计划(2003?2012)：Prediction in Ungauged Basins (PUB)相关科学问题，开展了长期系统性研究。近期，团队李雪莹博士生、龙笛教授及合作者，以土壤水动态变化信息约束水文通量，基于高维参数空间优化的数学模型思想，构建了一种不依赖实测流量的径流求解模型(Soil Moisture to Runoff, SM2R)。模型在青藏高原江河源区20个测站进行了初步验证，表现了较好的精度和稳健性。
1
引言
径流是水循环中的重要变量。准确可靠的长序列径流数据，是水资源调查评价的基础。全球径流观测站点分布不均，许多区域尤其是高寒地区缺乏径流观测数据。水文研究中常利用水文模型模拟径流，但通常需要实测流量数据率定模型参数，在无测站、缺资料流域的应用受到限制。为解决上述问题，清华大学水利系遥感水文团队，围绕国际水文科学协会(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 。
2
研究区与数据
本研究以青藏高原7个外流区(包括印度河、恒河、雅鲁藏布江、怒江、澜沧江、长江、黄河)，共20个子流域为测试流域(图1)，利用SM2R模型模拟过去40年(1981?2020年)月尺度径流。通过对比实测径流、跨部门影响模型比较计划(ISIMIP)中16种模式输出径流均值、再分析径流数据、基于实测资料率定参数的VIC模型模拟的径流，评估了SM2R模型在青藏高原不同气候地理特征的流域的模拟效果。
图1 青藏高原地理位置，径流观测站、河流、冰川、湖泊及流域分布
SM2R模型的驱动数据包括降雨、降雪、潜在蒸散发、土壤体积含水率(后文中简称 “土壤含水率”)，上述变量均由全球再分析数据集ERA5-Land(ERA5L)获取。冰川广泛分布的流域需考虑冰川净质量变化，可由Hugonnet等 [5] 研制的全球冰川质量变化数据获取。此外，SM2R模型利用气温、积雪覆盖面积、土壤性质数据作为辅助数据，用于确定流域特征和进行参数初始化。其中，气温由ERA5L数据获取；积雪覆盖面积由Immerzeel等 [6] 基于MODIS遥感产品，针对全球78个水塔(包括青藏高原及天山区域)研制的多年平均积雪覆盖面积提取；土壤性质数据由全球土壤数据库(HWSD)获取。
用于对比SM2R模拟效果的径流数据集如下：(1)ISIMIP 中H08、 WaterGAP2、MATSIRO和LPJml的多模型平均结果，其中，每个模型分别由GFDL、Hadgem、IPSL和Miroc5四种气候模式的输出结果驱动，ISIMIP可提供1981?2005年历史时段的月径流模拟值；(2)ERA5L再分析径流数据；(3)Gou等 [7] 基于VIC模型并利用国内约200个径流站的实测径流率定参数，模拟研制的1961?2018年中国天然径流数据集(CNRD)。
3
研究方法
SM2R模型的框架如图2所示，关键模块包括：(1)考虑冰川和积雪模块的土壤水量平衡，(2)建立水文通量(蒸散、下渗、及径流)与土壤体积含水率的函数约束。以系统状态变量达到最优为指标，即水量平衡求解的土壤含水率变化和ERA5L估算的土壤含水率变化最接近，利用基于梯度下降的高维参数空间优化方法，求解各参数并估算径流。
图2 SM2R模型框架图
3.1 建立土壤水量平衡
以一定深度的土层为研究对象，土壤水量平衡公式如下：
式中， z 表示土壤深度(mm)，在SM2R模型作为待求参数之一，是满足水量平衡闭合的参数； θ 表示土壤含水率(m 3 /m 3 )； P Eq 表示土壤包气带的等效输入水量(mm/mon)； L ( θ )表示水量输出通量(mm/mon)，包括蒸散( ET ( θ ))、下渗( D ( θ ))和径流( R ( θ ))(式(2))，本研究中径流由快速流 R fast 和慢速流 R slow 两部分组成(式(3))。
在青藏高原，降水相态和冰川质量变化等对径流模拟有重要影响。SM2R模型在计算等效输入水量 P Eq (式(1))时，将总降水划分为液态降水(降雨)和固态降水(降雪)，冰川积累/消融对土壤包气带等效输入水量的影响为：(1)冰川积累期，固态降水(降雪)首先储存在冰川上，其余形成土壤层的有效输入(式(4))；(2)冰川消融期，冰川融化为液态水，同液态降水(降雨)共同作为土壤层的有效输入(式(5))。
式中， dh g 表示冰川高程变化(mm/mon)； snowfall g 和 rainfall g 分别代表考虑冰川质量变化后的等效降雪和等效降雨； dh 表示由冰川净质量变化引起的流域平均水深变化，由式(6)计算：
式中， Area g 指流域内冰川面积(km 2 )；Area指流域集水面积(km 2 )； ρ g 指冰川密度(取850 kg/m 3 )，ρw指水密度(取1000 kg/m 3 )。
此外，高海拔地区(例如喜马拉雅山脉)的实测径流可能滞后于降水数月，且径流与降水出现峰值的时间不同，这是由于冬春季的固态降水通常在冬春季冻结，当温度升高后再融化成液态水参与到水循环中。本研究将具有“储存冬春季降雪并补偿夏季降雨”的流域定义为“降雪主导”流域，提出积雪指数(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模型的目标函数为：由水量平衡通量计算的土壤含水率变化，与由ERA5L估算的土壤含水率变化的均方根误差(式(19))：
式中， n 表示样本数目， X 表示待求的参数空间，即 X =[ a , b , c , e , f , g , z ]。SM2R模型设置优化目标变量(土壤含水率变化， dθ / dt )量级的1‰为优化阈值(式(20))。当两次优化的效果提升小于优化阈值时(式(21))，认为继续优化对模拟结果提升有限，即可终止迭代过程。
采用基于梯度下降的高维参数空间优化方法求解各参数，具体步骤如下：
4
主要研究结果
本研究利用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的评估结果)
分析SM2R模拟径流在不同气候地理特征流域的模拟效果发现(图4)：印度河(站点1?4)、恒河(站点5?7)、雅鲁藏布江(站点8?9)流域，SM2R模拟径流整体表现较好。在印度河流域，由于模型考虑了高山区径流和降水的滞后时间，因此，有效提高了模拟径流和实测径流的相关性( r ≥ 0.94 (SM2R)； r = 0.67?0.71 (ERA5L))。恒河和雅鲁藏布江流域受南亚季风影响显著，降雨产流是总径流的主要贡献，SM2R模型能反映季风主导影响下较稳定的降雨径流关系。
怒江(站点13?14)和澜沧江(站点15?16)同样受南亚季风影响，但这两个流域下游地区的河道狭长(特别是贡山站的集水流域)。驱动数据ERA5L(空间分辨率为0.1°(~10 km))在狭长的集水流域不确定性较大，导致SM2R径流模拟精度降低。SM2R模拟的径流在怒江和澜沧江流域总体表现较好，但在贡山站(站点13)的统计指标较低。
长江下游泸宁站(站点18)受东亚季风影响显著，SM2R模拟径流的各指标表现良好。长江源(直门达水文站集水流域，站点17)和黄河源(唐乃亥和玛曲水文站集水流域，站点19和20)位于青藏高原腹地，多年冻土和季节性冻土分布广泛。气候变化下冻土退化，引起土壤含水率、土壤性质和水文过程的一系列改变，对产流过程的影响复杂，SM2R模拟精度较低。
图4 SM2R模拟径流在不同站点的相关系数(a)和纳什效率系数(b)评估指标
对比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估算径流在直门达、唐乃亥及玛曲站的评估结果)。
5
讨论
本研究构建的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模型中各水文通量独立计算，即每种输入数据的不确定性，仅通过相关的数学函数影响特定分量和参数，可减小误差传递，提高模型鲁棒性。
图6 ERA5L和陆面模式估算的2003?2020年土壤体积含水率(左坐标轴)和土壤体积含水率变化(右坐标轴)
6
总结
本研究基于高维参数空间优化的数学思想，构建了以土壤水动态变化约束径流求解的新方法(SM2R)，实现了在不依赖实测径流资料条件下，较为准确地估算了青藏高原不同气候和地理特征流域的长时序径流。SM2R模型的输入数据和模型参数较少，径流模拟效果显著优于ERA5L和ISIMIP的径流估算，与利用实测径流率定参数的VIC模型表现接近。本研究提出的SM2R模型提供了“由状态量推求通量”的径流模拟新思路，有潜力在其他无测站、缺资料流域的径流估算和分析中提供参考借鉴。
本研究受到第二次青藏高原综合科学考察(2019QZKK0105)、基金委“西南径流重大研究计划”(92047301、92047203)等项目支持，谨致谢忱。
