1 引言
青藏高原约占中国陆地面积的四分之一， 平均海拔超过4000 m， 素有“世界屋脊”和“第三极”之称（Qiu， 2008）， 其巨大的地形和高海拔通过热力和动力强迫作用， 使得青藏高原对东亚乃至全球的天气及气候产生重大影响（Duan et al， 2012； Wu et al， 2015）。同时， 青藏高原广泛分布着多年冻土和冰川等， 是包括长江、 黄河在内的众多大江大河的发源地， 也被称为“亚洲水塔”（Immerzeel et al， 2010）， 在区域水文过程中发挥着重要作用。在全球气候变暖背景下， 青藏高原正经历着显著变暖， 尤其是20世纪90年代以来， 其升温速率要远快于东部季风区和西北干旱区（Kuang and Jiao， 2016； 杜勤勤等， 2018）。据研究， 青藏高原年平均气温（1970 -2014年）的升温速率约为0.35 ℃·（10a）-1（Yao et al， 2019）， 而年平均地表温度（1979 -2018年）升温速率可达0.46 ℃·（10a）-1（Yan et al， 2020）， 均显著高于北半球陆地平均升温速率。快速升温引起高原气候生态环境尤其是冰冻圈要素发生了一系列显著的变化， 如冰川萎缩（Yao et al， 2012； Yang et al， 2019）、 冻土退化（Cheng and Wu， 2007； Zhao et al， 2020）、 积雪面积减少（车涛等， 2019）等， 这必将会影响高原土壤-植被-大气系统的水热平衡， 进而对该区域水文（Immerzeel et al， 2010）、 生态（Shen et al， 2015）和气候（Ding et al， 2020）等产生深刻的影响。
蒸散发是水循环过程和地表能量平衡的重要组成部分， 地表蒸散发既是联系土壤-植被-大气连续体中水循环的纽带， 也是地表物质交换与能量平衡的关键环节（Wang and Dickinson， 2012）。过去几十年里， 围绕蒸散发的观测、 模拟等研究已有了很大的进展。观测方面， 主要方法包括蒸渗仪、 涡动相关法及闪烁仪等， 空间范围从几米到几十公里不等， 这些方法在点尺度或小范围空间上得到了广泛的应用（Liu et al， 2018）， 但难以获取大尺度空间的蒸散发量。模拟方面， 与蒸散发相关的模型众多， 大致包括遥感蒸散发模型， 比较有代表性的模型如： SEBS（Surface Energy Balance System）模型、 SEBAL（Surface Energy Balance Algorithm for Land）模型以及MOD16产品等； 陆面模型， 常使用的模型有Noah， CLM模型等； 以及侧重生态水文过程对蒸散发影响的模型， 如PML（Penman-Monteith-Leuning）， BEPS（Boreal Ecosystem Productivity Simulator）模型等。此外， 还有一些时空升尺度技术， 如： 机器学习和数理统计方法等， 以上的模型和时空升尺度技术在一定程度上突破了传统观测方法在时空方面的限制， 使得获取长时间， 大尺度空间的蒸散发变化成为可能， 也是目前区域蒸散发研究的主要手段（Mu et al， 2007； Jung et al， 2010； Ma et al， 2014a； Ma et al， 2017a； Zhang et al， 2019）。然而， 上述模型大都存在模型结构复杂， 参数众多且难以获取的问题， 升尺度技术存在机理欠缺的不足， 因而在应用时不同方法计算结果的差异性很大， 这也是当前估算区域蒸散发的一大瓶颈。在青藏高原地区， 过去在高原东部和南部陆续建立了一批综合观测点（Ma et al， 2008， 2017b）， 这些观测主要集中在季节冻土区， 蒸散发相关模拟的改进和验证工作也主要是针对季节冻土区展开（Ma et al， 2014b； Song et al， 2017； Zou et al， 2018）， 而对于广泛分布的多年冻土区， 蒸散发相关的观测或是模拟研究极其有限（Chang et al， 2018； 王利辉等， 2019； Wang et al， 2020a）。已有研究指出， 多年冻土区和季节冻土区的地表能量分配存在明显不同（Gu et al， 2015）， 因而有必要专门针对多年冻土区的蒸散发开展研究。
考虑到高原多年冻土区蒸散发相关的观测资料非常缺乏， 本研究中采用基于互补原理的方法模拟地表实际蒸散发， 该方法最大的优势是不需要植被， 土壤水分等数据， 需要输入的模型参数相对较少， 只需要常规的气象资料即可。互补理论最早由Bouchet（1963）提出， 现已被应用于解释“蒸发悖论”现象以及估算不同时空尺度实际蒸散发， 近些年来受到了越来越多的关注（Burtsaert and Parlange， 1998； Brutsaert， 2013； Brutsaert et al， 2017， 2020； Ma et al， 2019）。在Brutsaert（2015）提出广义非线性互补模型前， 比较有影响的互补蒸散模型有AA（Advection-Aridity）模型（Brutsaert and Stricker， 1979）， CRAE（complementary relationship areal evapotranspiration）模型（Morton， 1983）以及GG模型（Granger and Gray， 1989）， 这些模型均基于传统互补蒸散理论， 其核心思想是假设在最小平流影响下， 区域地表均质且足够大， 地表能量的输入一定时， 潜在蒸散发会随着实际蒸散发的变化而产生负反馈作用， 即在地表水分充足不受限制的情况下， 实际蒸散发与表观潜在蒸散发相等， 同时也等于湿润环境蒸散发。随着地表变干， 实际蒸散发的减少量与表观潜在蒸散发的增加量相等， 呈现出对称互补的关系， 因而传统互补模型是线性对称的。此后， 一些研究（Kahler and Brutsaert， 2006； Szilagyi， 2007）发现互补关系并不是对称的， 进而提出了非对称线性互补模型。然而， Han et al（2011）的研究表明无论是传统互补模型还是非对称线性互补模型， 其普适性都较差， 只在一些气候过渡区适用。为了解决该问题， Han et al（2012）首先将Budyko水热平衡理论中的边界约束条件引入到互补模型中， 并推导了一个新的互补模型解析式， 在干旱区验证时取得了较好的结果； 在此基础上， Brutsaert通过进一步研究提出了广义非线性互补模型（Brutsaert， 2015）， 随后， Crago et al（2016， 2018）和Szilagyi et al（2017）完善了广义非线性互补模型边界约束条件的不足； Han and Tian（2018）也在此前的基础上， 提出了一个具有物理基础的S型函数。最近， Gao and Xu（2021）提出了一个指数型函数， 结果表明该函数较前两者能更好地模拟出蒸散发变化。总体而言， 互补理论模型在不同植被类型下垫面和不同气候区已有了不少研究（Liu et al， 2016， 2018； Zhang et al， 2017； Hu et al， 2018； 郭晨露等， 2017）， 但在多年冻土区的研究还非常少， 对于互补模型在多年冻土区的适用性缺少系统性的评估工作， 这限制了互补理论模型在多年冻土区的发展和应用。
基于此， 本文利用位于青藏高原多年冻土区腹地唐古拉观测场的气象和通量等观测资料， 计算了实际蒸散发量和潜在蒸散发量等， 分析了三种广义非线性互补函数和S型函数模拟日尺度地表实际蒸散发的能力， 最后选择参数较多的两种模型进行了参数的敏感性分析。本文的研究结果可以为互补理论在高原多年冻土区的适用性提供一定的参考， 同时能为该地区的地表实际蒸散发估算提供一种新的可行的选择。
2 观测站点和数据
2.1 观测站点概况
唐古拉观测场（33°04′N， 91°56′E， 图1 ）位于青藏高原多年冻土区腹地， 海拔5100 m， 属于连续多年冻土区， 活动层厚度约为3.3 m， 下垫面为退化的高寒草甸。该地区年平均气温约为-4.9 ℃， 年降水量约为436.7 mm， 降水主要集中在5 -9月（Sun et al， 2020）。本文使用的唐古拉综合观测场资料来源于青海藏北高原冰冻圈特殊环境与灾害国家野外科学观测研究站。
2.2 观测数据
观测资料包括近地面空气温湿度（2， 5， 10 m）、 风速（2， 5， 10 m）、 辐射四分量、 降水以及浅层5 cm处土壤温湿度和土壤热通量数据。此外， 还使用了涡动相关系统观测得到的感热和潜热通量数据， 该系统架设高度为3 m， 由一套三维超声风速仪（CSAT3， Campbell， USA）和开路式CO2/H2O气体分析仪（Li-7500， Li-Cor， USA）组成， 原始10 Hz高频数据通过CR5000数据采集器采集。本研究中使用的数据时间段为2010年1月1日至2011年12月31日， 经处理后时间分辨率均为30 min， 最后转换成逐日数据。研究时段内数据整体质量较好， 缺失数据较少。
3 数据处理和研究方法
3.1 通量数据处理
本文使用EddyPro软件处理涡动相关原始10 Hz高频数据， 主要处理步骤包括野点值剔除， 延迟时间校正， 坐标旋转， 超声虚温修正和密度（WPL）修正。对输出后的30 min数据进行质量控制， 剔除仪器误差数据和摩擦速度小于0.1 m·s-1的数据， 以及降水时段的数据， 剔除的感热和潜热通量数据分别占到使用时段所有数据的21.76%和26.86%。对于质量控制后缺失的数据， 采用边缘分布采样法和空气动力学等方法填补。能量闭合率是评价涡动数据质量的重要指标， 通过对研究时段内观测点能量通量数据进行闭合率分析， 得到唐古拉观测场地区的能量闭合率约为0.82， 通常情况下能量闭合率的可接受范围在0.7~0.9， 说明本研究中使用的观测场涡动数据质量良好。
需要指出， 当地表有积雪覆盖时， 此时的蒸散发形式主要是积雪升华， 由于积雪观测数据缺失， 本文通过反照率albedo来判断地表是否有积雪覆盖：
式中： USR为观测的向上短波辐射； DSR为观测的向下短波辐射。为了排除太阳高度角的影响， 日均反照率采用10:00（当地时， 下同） -14:00的计算结果。从图2 中可以看出， 地表无积雪时的反照率较平稳， 当降雪时， 反照率会迅速增大， 随着积雪升华或融化， 反照率逐渐减小。由图2 还可发现， 0.3大致可以作为区分地表是否有积雪覆盖的阈值， 因此， 当albedo>0.3时， 认为地表有积雪覆盖； 当albedo<0.3时， 地表无积雪覆盖。此外， 我们发现研究时段内唐古拉站夏季仍有降雪事件发生， 地表会覆盖有较薄的积雪， 此时升华量很小， 因此， 本文只计算冷季（10月至次年4月）的积雪升华量。最终， 地表蒸散发（非积雪覆盖日）与升华（冷季积雪覆盖日）分别利用蒸发潜热λv（单位： J·kg-1）和升华潜热λs（单位： J·kg-1）结合涡动观测及插补得到的潜热通量折算成蒸散发量和升华量（Zhu et al， 2014）：
式中： Ta 为2 m高度处气温（单位： ℃）。
此外， 地表土壤热通量G使用浅层5 cm土壤热通量和土壤温湿度通过一维热传导方程（TDE法）（阳坤和王介民， 2008； 杨成等， 2020）计算， 具体计算过程见Yao et al（2011）的研究。
3.2 互补蒸散发模型
传统互补模型（Bouchet， 1963）认为在大且均匀、 受平流影响较小的下垫面， 可利用能量一定时， 在地表充分供水条件下， ， 其中E为实际蒸散发， 为表观潜在蒸散发， 表示在非潜在大气条件下从非饱和陆地表面内部的一个小的饱和面发生的蒸散发， 为湿润环境蒸散发， 表示发生在饱和的陆地表面的蒸散发。随着地表供水量的降低， E的减少量等于 的增加量， 其形式如下：
式中： α是Prisetley-Taylor系数， 默认值为1.26； Δ是当前气温对应的饱和水汽压曲线斜率（单位： kPa·℃-1）； γ是干湿表常数（单位： kPa·℃-1）； 和G分别表示净辐射和地表土壤热通量（单位： mm·d-1）； e*和 分别为气温对应的饱和水汽压和实际水汽压（单位： kPa）； 为风函数（单位： mm·d-1·kPa-1）， 本文中使用Rome风函数公式计算， 公式如下：
需要指出的是， Szilagyi and Jozsa（2008）认为在计算 时， 饱和水汽压曲线斜率Δ应该对应湿润环境的温度而不是测量到的空气温度， 并提出采用湿润环境空气温度下的饱和水汽压斜率ΔT wea（单位： kPa·℃-1）。因T wea在实际情况下通常无法获得， Szilagyi and Jozsa（2008）提出采用湿润环境表面温度T wes来近似， T wea可以用波文比 迭代求解：
其中： 是湿润环境表面的饱和水汽压， 迭代求解的 可能会大于 ， 此时须令 ， 修正后的公式为：
3.2.1 B2015模型
Brutsaert（2015）在传统互补模型的基础上， 基于严格的物理边界条件重新定义了互补理论中蒸散发相关变量的概念， 描述了极端干湿的物理边界条件， 提出了多项式广义互补蒸散模型， 其物理边界形式为：
基于上述的物理边界条件， Brutsaert（2015）给出了简化的多项式函数形式：
本文中简称为B2015基础式。此外， Brutsaert（2015）还提出了适用性更为广泛的一般多项式函数形式， 本文中简称B2015一般式， 公式如下：
式中： c为需要率定的参数， ； 其中 。
3.2.2 C2016模型
Crago et al（2016）认为B2015模型物理边界条件中y=0， x=0存在着缺陷， 因为PT公式计算的 不会无限小， 实际情况无法满足x=0的假设。因此， Crago et al（2016）引入最大潜在蒸散发E pmax的概念来计算x min， 重新定义Brutsaert（2015）提出的物理边界， 即y=0时， x=x min， 其中：
根据上述调整的物理边界条件， Cargo et al（2016）重新定义了广义非线性互补函数形式， 本文中简称C2016， 其解析形式为：
3.2.3 S2017模型
Szilagyi et al（2017）也同样认为B2015模型物理边界条件中y=0， x=0存在着缺陷， 并将E pmax定义为E=0时的E pa， 使用干燥表面空气温度T dry采用Penman公式计算得到：
T dry计算方法见参考文献（Szilagyi et al， 2017）。在Brutsaert（2015）提出的简化的多项式函数形式上， Szilagyi et al（2017）重新定义了函数中的自变量， 本文中简称S2017， 其解析形式为：
C2016和S2017模型的解析式中简化了未知参数， 不再包含参数c， 在使用时只需要率定参数α即可。
3.2.4 H2018模型
Han and Tian（2018）在前期研究（Han et al， 2011， 2012）的基础上， 重新定义互补模型物理边界的约束范围， 提出了S型广义互补蒸散模型， 本文中简称H2018， 其解析形式为：
式中： ， 的计算同 ， E rad的计算见公式（5） 。H2018的物理边界条件为：
式中： m和n可以通过AA模型参数α和b计算：
3.3 模型的评价指标
本研究选用的评价指标有均方根误差（Root Mean Square Error， RMSE）、 纳什效率系数（Nash-Sutcliffe efficiency coefficient， NSE）以及King-Gupta效率系数（King-Gupta efficiency coefficient， KGE）（Gupta et al， 2009； Kling et al， 2012）。各指标的计算公式如下：
其中，
式中： N为样本数； μ为平均值； σ为标准差； 下标sim为模拟结果。
4 结果分析
4.1 唐古拉地表能量和气候特征
从2010 -2011年唐古拉地表蒸散发及升华量的逐日变化（图3 ）中可以看出， 根据近地表5 cm处土壤温度T soil5的日变化， 可将冻融过程划分为完全冻结（T soil5_max<0 ℃）、 日冻融循环（T soil5_max>0 ℃且T soil5_min<0 ℃）以及完全融化（T soil5_min>0 ℃）三个阶段。在完全冻结阶段， 土壤处于完全冻结状态或地表有积雪覆盖， 此时实际蒸散发以冰雪升华形式发生， 地表日均蒸散发量约为0.2 mm·d-1， 2010年和2011年此阶段地表蒸散发量分别为29.6 mm和21.3 mm； 积雪日时， 升华量介于0.03~1.4 mm·d-1， 2010年和2011年此阶段地表积雪升华量分别为4.7 mm和11.6 mm。4月中旬， 随着气温的升高， 土壤出现昼融夜冻现象， 此时地表实际蒸散发速率迅速增加， 日均蒸散发量可达1 mm·d-1以上。5月中旬至10月中旬， 土壤由地表向下逐渐融化， 直到达到最大融化深度， 土壤中液态含水量增加， 同时伴随雨季产生大量降水， 该时期实际蒸散发以液态水蒸发为主， 日均蒸散发量约为2.4 mm·d-1， 2010年和2011年此阶段地表蒸散发量分别为350.5 mm和363.9 mm， 约占全年蒸散发总量的82.9%和85.4%。10月中旬至月底， 存在一个短暂的日冻融循环期， 此时土壤液态含水量迅速减少， 加之雨季结束， 地表蒸散发量明显减少。此后， 随着气温进一步降低， 土壤又进入完全冻结阶段， 地表蒸散发随之减小， 并降到最低值。
将研究时段内蒸散发及升华量和降水量对比后发现， 冷季时， 降水量整体上小于实际蒸散发量， 这可能是由于土壤中存在的未冻水， 为地表蒸散发提供了一定的水分； 暖季时， 降水量略大于实际蒸散发量， 因唐古拉观测场整体较为平坦， 降水难以产生地表径流， 因此该地区降水主要以地表蒸散发的方式返回大气。年尺度上， 2010年和2011年的累积地表蒸散发和升华量总和分别为433 mm和446.5 mm， 同期降水量均为394 mm。上述地表蒸散发和升华总量大于降水总量的可能原因是受高寒恶劣环境的影响， 观测场涡动观测的蒸散发数据存在较多的缺失， 本文通过插补的方法填充会给最终的结果带来一定的偏差， 尤其是冷季， 偏差更为明显。此外， 高原上降水观测也存在着很大的不确定性。因此， 将本文的观测结果与王利辉等（2019）在距离观测场约20 km处的唐古拉山冬克玛底河观测点的实测结果进行了比较（表1 ）， 结果表明， 无论是实际蒸散发量还是降水量， 两者在植被生长季的结果都是比较接近的， 尽管与观测场真实情况会存在着一定的偏差， 但本研究的实际蒸散发和降水数据仍是可靠的。
表1 本研究观测数据与其他观测结果的对比Table 1 Comparsion of observation data of this study and other observation results |
| 观测点 | 时间段 | 观测值 | 来源 |
|---|---|---|---|
| 唐古拉山冬克玛底河流域 （33°03′N， 92°02′E） | 2010年6月17日至9月24日 | 实际蒸散发量： 280.3 mm 降水量： 342.9 mm | 王利辉等（2019）， 小型蒸渗仪观测 |
| 2011年6月17日至9月24日 | 降水量： 237.2 mm | ||
| 唐古拉综合观测场 （33°04′N， 91°56′E） | 2010年6月17日至9月24日 | 实际蒸散发量： 265.3 mm 降水量： 290.9 mm | 本研究， 涡动系统观测 |
| 2011年6月17日至9月24日 | 降水量： 275.9 mm |
4.2 不同互补模型模拟结果评价
本文使用四种不同的互补模型分别采用默认参数以及率定参数模拟地表实际蒸散发， 比较了四种互补模型模拟实际蒸散发的能力， 旨在探讨互补理论在青藏高原多年冻土区的适用性， 以及认识该地区的互补模态是如何的。本研究中将2010年的逐日蒸散发数据作为率定数据， 然后利用2011年的日蒸散发数据作为验证数据， 着重分析不同互补模型在青藏高原多年冻土区模拟日地表实际蒸散发的能力和不确定性。
4.2.1 利用默认参数时模型模拟结果评价
本文首先使用默认的参数来模拟唐古拉站2011年的地表实际蒸散发， B2015模型默认参数分别为α=1.26， c=0； C2016和S2017模型默认参数为α=1.26； H2018模型默认参数分别为α=1.26， b=1， x max=1和x min=0。图4 展示了不同互补模型定义的自变量x（或X）与因变量y的关系， 从图4 中可以看出， 不同互补函数使用默认的参数值模拟的结果存在着较大差异。B2015模型模拟的结果［图4 （a）中实线］远远偏离了实测资料计算的结果［图4 （a）中散点］， 根据Brutsaert广义非线性互补模型的前提假设， 即当下垫面接近饱和时， ， 然而使用默认参数时， 和 均明显大于E。C2016和S2017及H2018模型模拟的结果整体均高于实测资料计算的结果， 但偏差较B2015模型模拟结果要小， 以上结果表明使用四种模型默认的参数值并不适合多年冻土区实际蒸散发的模拟。
将模拟的2011年日地表实际蒸散发量和实测的2011年日实际蒸散发量点绘于图（图5 ）可以看到， 四种模型模拟的实际蒸散发量均集中在1∶1线上方， 说明模拟的蒸散发量高于实测蒸散发量。已有研究（杨汉波等， 2008， 2009）指出， 参数α在互补模型中对模拟精度有重要的影响， α值的大小与纬度和离海洋的距离有关， 一般而言， α值随着纬度的增加而增加， 随着离海洋距离的增大而减小。结合图4 来看， 默认的参数α值（1.26）在中低纬且居青藏高原腹地的多年冻土区偏大， 进而造成模型计算的湿润环境蒸散发量 偏大， 使得最终的实际蒸散发结果高估。根据KGE值来看， 模拟精度为S2017（0.85）>H2018（0.79）>C2016（0.7）>B2015（0.27）。C2016、 S2017和H2018模型均是在B2015模型基础上改进的， 从以上结果可以看出， 在使用相同默认参数的情况下， C2016、 S2017和H2018模型模拟日地表实际蒸散发的能力要远远好于B2015模型， 这表明过去一系列针对物理边界约束条件和函数形式的改进， 使得互补理论模拟实际蒸散发的能力确实得到提高， 本文的结果也印证了此点。另外， 值得注意的是， C2016和S2017模型中均只有一个未知参数， 需要率定的参数少于B2015模型和H2018模型， 而模拟精度却较高， 这对于实测数据相较缺乏的青藏高原多年冻土区而言， 在模拟地表实际蒸散发时具有一定的优势。
图5 使用不同互补函数形式默认参数值模拟得到的2011年逐日实际蒸散发量与观测的实际蒸散发量散点分布 虚线为线性拟合结果, 实线为1∶1线, N表示使用的样本数Fig.5 Scatter plots of simulated daily actual evapotranspiration used by different complementary functions with default parameters and observed daily actual evapotranspiration in 2011.The dash lines show linear fitting results， the solid lines show one-to-one relationship， N represents used sample numbers |
4.2.2 利用校正后的参数时模型模拟结果评价
根据前面的研究， 得出在应用互补模型模拟地表蒸散发前， 首先需要对互补模型中的参数进行率定， 本文通过计算模拟的2010年日实际蒸散发E sim值与观测的日实际蒸散发E值之间的最小均方根误差（RMSE）来确定模型中参数的最优解。同时， 参考前人（Liu et al， 2016； Wang et al， 2020b）的研究， 将参数α、 c和b的范围分别设置为［1， 1.5］， ［-1， 4］和［1， 4］。图6 展示了四种不同互补函数率定参数后模拟（实线）与实际观测（散点）的结果。从图6 中可以看出， 率定参数后模拟的实际蒸散发值相较未率定参数模拟的结果有了很大的提高， 自变量x（或X）和y间的关系基本能够模拟出来。尽管各模型率定后的参数值不一致， 但分析后发现仍存在一定的规律。以参数α为例， 率定后B2015模型中为1.02， C2016和S2017模型中为1.13， H2018模型中为1.06， 这与Ma et al（2015）、 Yang et al（2020）在高原冻土区计算的结果十分接近， 其α值分别为1.13和1.19。此外， H2018模型中参数b值（1.62）也与Wang et al（2020b）给出的相似下垫面类型的b值（0.98±0.48）较接近。因此， 互补理论中的参数值在相似下垫面空间尺度上可能具有一定的可移植性， 这可以为缺资料地区估算地表蒸散发时的参数值提供一定的参考。另外， 有研究（杨汉波等， 2008， 2009； Wang et al， 2020b）指出， 参数α和b值还存在着时间变异性， 即在不同时间尺度上参数值的差异大， 且该差异要大于不同年份之间的差异， 因此， 在模拟不同尺度（日、 月、 年）的时候还需要注意区分时间尺度来获取可靠的参数值。
为了识别率定后的参数是否能够模拟实际蒸散量及其变化， 以及和使用默认参数相比， 结果如何变化两个问题， 本文首先计算2011年逐日潜在蒸散发 和湿润环境蒸散发 ， 再和各率定后的参数值一同代入互补模型函数式中， 即可估算2011年逐日的实际蒸散发量， 最后和未率定参数模拟的结果（图4 ， 图5 ）进行对比。图7 展示的是四种不同互补函数形式采用率定的参数值模拟的2011年逐日蒸散发和实测的逐日蒸散发值散点分布。从图7 中可以看出， 各模型模拟的结果和观测值均接近于1∶1线， B2015、 C2016、 S2017和H2018模型模拟实际蒸散发的KGE值（表2 ）分别为0.82、 0.89、 0.87和0.84， 各模型模拟结果间的差异较小。进一步分析发现， 除S2017模型外， 其余各模型模拟的结果均显示出在日蒸散发量E>2 mm时， 数据较为集中地分布在1∶1线附近， 而日蒸散发量E<2 mm时， 数据分布很零散， 偏离1∶1线， 该时期对应着在土壤完全冻结期和日冻融循环期， 以上结果说明完全冻结期和日冻融循环期模拟结果的不确定性要大于完全融化期， 精度上低于后者。
图7 不同互补函数形式采用率定参数值模拟得到的2011年逐日实际蒸散发量与观测的实际蒸散发量散点分布 虚线为线性拟合结果， 实线为1∶1线， N代表样本数Fig.7 Scatter plots of simulated daily actual evapotranspiration used by different complementary functions with calibrated parameters and observed daily actual evapotranspiration in 2011.The dash lines show linear fitting results， the solid lines show one-to-one relationship， N represents used sample numbers |
表2 不同互补模型在冻融各阶段模拟的地表实际蒸散发量与实测的对比Table 2 Comparsion of simulated E by different complementary models and observed E under each freezing-thawing processes |
| 冻融阶段 | 率定前 | 率定后 | ||||||
|---|---|---|---|---|---|---|---|---|
| B2015 | C2016 | S2017 | H2018 | B2015 | C2016 | S2017 | H2018 | |
| 完全冻结期 | 98.1（360.6%） | 51.3（140.8%） | 26.5（24.4%） | 43.6（104.7%） | 46（116%） | 42.7（100.5%） | 19.2（-9.9%） | 47.4（122.5%） |
| 日冻融循环 | 99.5（143.3%） | 76.7（87.5%） | 63.7（55.7%） | 74.7（82.6%） | 62.5（52.8%） | 62.2（52.1%） | 47.2（15.4%） | 62.6（53.1%） |
| 完全融化期 | 457.1（25.6%） | 447.1（22.9%） | 431.2（18.5%） | 402.8（10.7%） | 330.6（-9.2%） | 362.6（-0.36%） | 339.5（-6.7%） | 335.9（-7.7%） |
2011年1月1 -5日实测资料缺失， 未统计在内。完全冻结期： 2011年1月6日至4月8日、 2011年10月24日至12月31日； 日冻融循环期： 2011年4月9日至5月16日、 2011年10月12 -23日； 完全融化期： 2011年5月17日至10月11日； 单位： mm； 括号内为相对偏差。 |
图9 不同互补模型参数率定前、 率定后模拟的2011年地表实际蒸散发与实测结果对比的箱型图Fig.9 Boxplot of comparsion of simulated E by different complementary models before calibration， after calibration and observed E |
表3 不同互补模型率定参数前后模拟的地表实际蒸散发结果与实测的误差情况Table 3 Error of simulated E by different complementary models before calibration， after calibration and observed E |
| 参数 | 模型 | α | c | b | RMSE/（mm·d-1） | NSE | KGE | N |
|---|---|---|---|---|---|---|---|---|
| 默认参数 | B2015 | 1.26 | 0 | - | 0.77 | 0.46 | 0.27 | 267 |
| C2016 | 1.26 | - | - | 0.47 | 0.8 | 0.7 | 274 | |
| S2017 | 1.26 | - | - | 0.44 | 0.82 | 0.85 | 274 | |
| H2018 | 1.26 | - | 1 | 0.46 | 0.84 | 0.79 | 341 | |
| 校正参数 | B2015 | 1.02 | 1.75 | - | 0.31 | 0.92 | 0.82 | 272 |
| C2016 | 1.13 | - | - | 0.31 | 0.93 | 0.89 | 318 | |
| S2017 | 1.13 | - | - | 0.37 | 0.9 | 0.87 | 318 | |
| H2018 | 1.06 | - | 1.62 | 0.31 | 0.92 | 0.84 | 293 |
加粗字体是指各指标统计地模拟实际蒸散发精度最优的模型， -表示无数据 |
4.3 不同互补模型中参数的敏感性分析
图10 B2015模型（a）和H2018模型（b）中一个参数变化±10%、 ±20%、 ±30%、 ±40%、 ±50%而引起的模拟的地表实际蒸散发的变化虚线表示模拟地表实际蒸散发变化的平均值， 竖实线上下端分别表示模拟地表蒸散发变化的最大值和最小值 Fig.10 The maximum， minimum， and mean changes in the simulated actual evapotranspiration （ ） values with perturbations （within the range of ±10% to ±50%） of α and c for B2015 model （a）， α and b for H2018 model （b） around their calibrated values， respectively.The dash line represents the average value of ， the upper and lower side of vertically solid lines represent the maximum and minimum values of ， respectively |
5 结论与讨论
通过青藏高原多年冻土区唐古拉观测场实测的气象和涡动数据资料， 评估了四种不同互补模型在日尺度上模拟地表实际蒸散发的能力， 并对模型中参数的敏感性进行了分析， 得出以下主要结论：
（1） 利用默认参数模拟地表实际蒸散发时， S2017模型模拟精度最高， H2018和C2016模型次之， B2015模型模拟精度最低， RMSE值分别为0.44 mm·d-1、 0.46 mm·d-1、 0.47 mm·d-1、 0.77 mm·d-1； 四种互补模型模拟的地表蒸散发量整体均高于实测值， 相较于完全融化期， 完全冻结期和日冻融循环期模拟的地表实际蒸散发量与实测值偏差更大， 这可能与参数的不合理取值有关。
（2） 校正参数后， 各模型模拟精度为C2016>S2017>H2018>B2015， 模拟的年内日地表实际蒸散发量与实测值接近， RMSE值均在0.3~0.4 mm·d-1， NSE值均在0.9以上， KGE值在0.8~0.9。日冻融循环期和完全融化期模拟的结果相较率定参数前有明显提高， 而完全冻结期模拟的结果则较差。表明各互补模型整体上能够模拟出青藏高原多年冻土区日地表蒸散发量及变化， 但在年内不同时期的模拟精度有较大差异。
（3） B2015和H2018模型中参数α的值对模拟结果很敏感， 当α的取值越大于合理值， 模拟的地表实际蒸散发量会越大， 而B2015模型中的参数c和H2018模型中的参数b对模拟结果的影响较小。因此， 合理的α取值是十分重要的。
本文对不同互补模型模拟青藏高原多年冻土区地表实际蒸散发的能力进行了评估， 结果可以为互补理论应用于该地区蒸散发的模拟提供一定的参考。然而， 受限于观测数据的时间跨度， 文中仅研究了日尺度实际蒸散发的模拟效果， 缺乏对不同时间尺度的研究。此外， 本文缺乏对模型中的参数进行详细分析， 已有研究表明， 互补模型中的α值与下垫面植被类型、 干燥指数（E pa /P）、 地理位置等因子具有紧密联系， 如何建立参数α与环境因子之间的关系， 对于合理的α取值具有重要的意义， 同时这也是互补理论应用的重要研究方向， 而在多年冻土区， 近地表冻融过程对参数α的影响如何， 还有待于探索。下一步研究将收集更多的资料， 完善以上研究过程中的不足。