青藏高原多年冻土退化对蒸散发的影响

季 芳，范林峰，匡星星，邹一光，郑春苗，3

(1. 哈尔滨工业大学环境学院，黑龙江 哈尔滨 150090；2. 南方科技大学环境科学与工程学院，广东 深圳 518055；3. 东方理工高等研究院，浙江 宁波 315200)

在过去的约170 a(1850—2019年)，全球年平均地表气温升高了约1.05 ℃[1]。根据CMIP6模式(第6次国际耦合模式比较计划)的不同预估情景，到2100年，预计全球年平均地表气温将继续增加2.23～4.95 ℃[1]。气候变化加剧导致全球水文循环和碳循环增强[2-5]，改变了水文-能量平衡，使得干旱区更加干旱、湿润区更加湿润[6]。蒸散发是水圈、大气圈、生物圈水能交换过程的关键变量，与降水一同对地表水的可用性和分布起着重要的调节作用[7]。青藏高原是长江、黄河、雅鲁藏布江等亚洲重要河流的发源地，被称为“亚洲水塔”。蒸散发、降水和径流量的平衡对下游的水资源供应有很大的影响，在亚洲水能循环中发挥着重要作用，对全球气候变化极其敏感[8-10]。由于青藏高原特殊的地理单元和较低的年平均气温，青藏高原约40%的地区上覆多年冻土[11]。由于气候变暖的加剧，多年冻土正在经历着广泛的退化[10-12]。冻土的存在阻碍了地表与大气间的物质、能量和水分交换，增加了地表水文过程的复杂性[13]。因此，准确量化地表蒸散发的时空动态变化有助于更好地理解多年冻土地区水文过程对气候变化的响应。

多年冻土退化的表现之一是活动层厚度的增加。活动层，即多年冻土区冬季冻结、夏季融化的地表冻土层[14]。活动层的夏季融化加之降水的入渗导致表层土壤含水量增加，使得蒸散发作用增强；但随着表层冻土的不断融化，活动层的不断加深，土壤水也可能通过未冻结的土层流向深部，使得地表土壤蓄水量减小，导致蒸散发降低[12-13]。此外，冰的融化和升华作用需要消耗更多的热量，冻土的冻融过程改变了原有的能量平衡，因此，在其他条件相同时，相比正常土层，理论上多年冻土区的蒸散发更小。如Wang等[15]结果表明，若忽略冰的冻融过程所需的热量，青藏高原的平均年蒸散发将被高估8.88%。此外，活动层通过植物的根系为植物提供水分和养分，气候变暖导致的活动层厚度增加可能导致植物根区的加深，进而导致植被蒸腾作用增强[16-17]。因此，气候变化导致的多年冻土活动层加深对蒸散发的影响非常复杂。目前，大部分相关研究尤其是冻土退化导致的植物根区的加深对蒸散发的影响机制仅停留在理论分析阶段，缺乏模型模拟的量化分析。

Budyko模型[18]是一种基于水分与能量耦合平衡方程的经验模型，傅抱璞[19]在Budyko模型的基础上发展的Budyko-Fu模型能够表征气候、植被和下垫面特征与水文过程的相互关系。仅需1个参数的Budyko-Fu模型在适应不同区域地表特征方面具有更强的灵活性，被广泛运用于流域尺度、长时间序列平衡态陆面水热过程模拟和机理分析[20-23]。随着学者们对模型参数的不断优化，Budyko-Fu模型的适用性变得更加广泛，越来越多地被应用于高空间分辨率、年尺度的水文过程分析[6，24]。定量表示冻土地区活动层变化所导致的下垫面条件的改变，并将其耦合至Budyko-Fu模型中，可以定量分析冻土退化导致的植物根区变化对蒸散发的影响。

本文基于Budyko-Fu模型，在充分考虑气候、植被的基础上，构建冻土退化所致下垫面条件改变的水热耦合模型，分析模型参数与活动层厚度的统计关系并建立适用于青藏高原多年冻土区的经验公式，探究青藏高原多年冻土地区蒸散发的时空动态演变趋势，并量化青藏高原蒸散发对冻土退化的响应。研究结果可以提高对变化环境下水循环演变机理的认识，为高寒区水能循环过程的研究提供支撑。

1 研究区概况

青藏高原位于中国西南部，平均高程大于4 000 m[25]，被称为“世界屋脊”(图1(a))。青藏高原多年冻土区土壤类型以砂壤土和壤土为主(图1(b))，其中砂壤土占比74.3%，壤土占比20.9%。青藏高原幅员辽阔，多年冻土区植被类型多种多样，西部植被类型以草原为主，其中西北部高山地带分布有高山植被；东部地区以草甸和灌丛为主；东南部主要为森林(图1(c))。1982—2018年间，青藏高原多年冻土区的年平均气温为-5.5 ℃，年平均降水量约为336.6 mm。气温和降水量空间分布高度不均，年平均气温由东南地区的20 ℃降低到西北地区的-15℃，年平均降水量由东南地区的1 200 mm减小到西北地区的120 mm(图1(d)和图1(e))。青藏高原多年冻土区1982—2018年的平均活动层厚度为2.35 m，分布规律大致与气温相同(图1(f))。图1所示数据来源详见2.3节，文中所有分布图的投影坐标系均为Albers Conical Equal Area坐标系，分布图中白色区域为季节性冻土区，不列为研究对象。

图1 青藏高原地理及气象要素分布Fig.1 Spatial distribution of geographical and meteorological elements on the Qinghai-Tibet Plateau

2 研究方法与数据

收集青藏高原1982—2018年水文气象及归一化植被指数(Normalized Differential Vegetation Index，NDVI)数据，利用Penman-Monteith公式计算潜在蒸散发(ET，p)，基于Budyko-Fu模型，建立考虑活动层厚度的青藏高原多年冻土区水热耦合模型，通过最小二乘法确定模型参数并进行校正；将校正的模型应用于青藏高原多年冻土区，与现有蒸散发数据集进行对比，并探究该地区蒸散发的时空分布与演变特征；通过设置活动层厚度不变情景，与实际情况对比，研究活动层厚度与蒸散发的关系，定量分析冻土退化对青藏高原多年冻土区蒸散发的影响。

2.1 Budyko-Fu模型

实际蒸散发(ET)是大气降水(P)和潜在蒸散发相互控制的结果，下垫面特征或自然景观会影响二者的控制平衡，傅抱璞[19]基于原始的Budyko曲线[18]利用量纲分析方法和微分方程理论讨论了下垫面特性或自然景观对蒸散发的影响程度，得到Budyko-Fu模型：

(1)

为了提高计算的准确度， 本文先计算了研究区的日尺度ET， p， 将日尺度ET，p加和得到年尺度潜在蒸散发。 日尺度ET，p(mm/d)采用考虑能量平衡和空气动力学的Penman-Monteith公式[26]计算：

(2)

式中：Δ为饱和蒸汽压曲线斜率，kPa/ ℃，由气温计算得到；γ为干湿常数，kPa/ ℃，由降水计算得到；Rn为净辐射，MJ/(m2·d)，由地面净太阳辐射和地面净热辐射计算得到；G为土壤热通量，MJ/(m2·d)，由于本文中ET,p计算时间尺度为1 d，故可认为G=0；T为气温，℃；u2为地面以上2 m风速，m/s，由10 m风速转换而来；es为饱和蒸汽压，kPa，由日最高气温和最低气温计算得到；ea为实际蒸汽压，kPa，由比湿、降水和水蒸气分子质量与干燥空气分子质量的比值计算得到。式中各项气象要素的具体计算方法请参考世界粮农组织灌溉和排水文件[27]，计算所用数据来源见2.3节。

已有研究表明，植被功能、土壤类型和降水强度等因素对Budyko-Fu模型的参数ω有显著影响，如Donohue等[28]基于澳大利亚各流域情况建立了Budyko-Choudhury-Porporato模型：

(3)

式中：θ为土壤有效含水量， m3/m3， 定义为田间持水量与凋萎系数之差， 其值与土壤类型有关；Ze为植物有效生根深度， m， 与植被类型有关；α为研究期降水强度泊松分布的平均值， m。

另外，植被覆盖度也会影响ω的取值，Li等[20]基于全球不同流域拟合了ω与植被覆盖度的关系：

ω=2.36M+1.16

(4)

由于多年冻土区的植物根系吸水能力受活动层厚度的影响[16]，活动层加深可能导致植物根区向下延深，因此本研究结合式(3) 、式(4) 中土壤有效含水量、植物有效生根深度、降水强度和植被覆盖度与参数ω的相关关系，建立综合考虑植被功能、土壤类型、气候条件和多年冻土退化情况的ω参数化方案，使其可适用于多年冻土地区：

(5)

式中：Zt为多年冻土地区活动层厚度， m；a和b均为该半经验公式的回归参数。

2.2 趋势分析法

采用Theil-Sen slope方法[30]分析了青藏高原多年冻土区降水、潜在蒸散发和模拟的蒸散发的变化趋势，同时采用Mann-Kendall[31]非参数检验方法进行显著性检验。

2.3 数据来源

图1中青藏高原土壤类型数据来自世界土壤数据库(Regridded Harmonized World Soil Database，HWSD)，http:∥dx.doi.org/10.3334/ORNLDAAC/1247；植被类型数据来自中国科学院地理科学与资源研究所资源环境科学与数据中心，https:∥www.resdc.cn/data.aspx?DATAID=122;近地面(2 m)气温、降水、气压、比湿和10 m风速数据采用青藏高原数据中心(https:∥data.tpdc.ac.cn/zh-hans/data/8028b944-daaa-4511-8769-965612652c49/)的中国区域地面气象要素驱动数据集(China Meteorological Forcing Dataset，CMFD);地面净太阳辐射和地面净热辐射数据来自欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts，https:∥cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=overview)的第五代欧洲再分析陆地成分全球数据集(the global dataset for the land component of the fifth generation of European ReAnalysis，ERA5-Land)。各数据名称、单位、来源、时间分辨率及空间分辨率等信息详见表1。

表1 用于计算ET，p的数据具体信息

NDVI数据来自美国国家海洋与大气管理局(National Oceanic and Atmospheric Administration，NOAA)的AVHRR卫星数据集(https:∥www.ncei.noaa.gov/data/avhrr-land-normalized-difference-vegetation-index/access/)，时间分辨率为逐日，空间分辨率为0.05°，时间跨度为1982—2018年。模型中需要的植物可吸水含量数据与土壤类型有关，经查阅[32]，不同土壤类型中土壤有效含水量值如表2所示。

表2 不同土壤类型中土壤有效含水量

基于经典的Stefan模型[33]，通过考虑不同含水量对土壤热性质的影响，计算了1982—2018年青藏高原多年冻土活动层厚度，所得结果与青藏高原多年冻土区181个钻孔观测站点记录的1995—2018年间329个观测数据相比吻合良好，相关系数(R)为0.90，均方根误差(ERMS)为0.25 m，说明Stefan模型适用于青藏高原多年冻土区模拟，因此将此数据用于本模型输入[34]。青藏高原多年冻土区的划分使用Zou等[11]的模拟结果。由于模型校正需要，本文搜集了青藏高原多个实测蒸散发数据用于参数确定，具体站点、位置、年份、数据来源等信息详见表3，空间分布位置见图1(f)。另外，采用青藏高原数据中心的8 km分辨率蒸散发数据集(http:∥data.tpdc.ac.cn/zh-hans/data/46fcf412-d675-4137-a81a-a4570e45d3c7/)与模拟结果对比。

表3 用于模型校正的站点详细信息

2.4 活动层厚度不变情景

为了定量探究冻土退化对青藏高原蒸散发的影响，本文设置了假设情景，即活动层厚度不变情景：其他水文气象因素与实际情况相同，但活动层厚度取1982—1990年的平均值。将假设情景下模拟的蒸散发与真实情景(对照情景)进行比较，由于对照和假设情景下的其他所有变量都是相同的，因此对照和假设情景下蒸散发的差异揭示了冻土退化导致的植物根区深度变化对蒸散发的影响。

3 结果与讨论

3.1 模型校验

采用最小二乘法，对式(5)中的拟合参数a和b不断调整，使得模型计算的实际蒸散发与站点实测蒸散发的相关系数最大，而均方根误差尽可能小。由此确定，当a=1.25、b=1.56时，实际蒸散发与实测蒸散发的相关系数为0.84，均方根误差为87.1 mm(图2(a))。这表明，校正后的ω-Zt函数关系在青藏高原多年冻土区有较好的适应性，耦合了活动层厚度的Budyko-Fu模型具有较好的模拟能力。

图2 模型校正及与其他数据对比Fig.2 Model calibration and comparison with other dataset

将校正后的模型应用于青藏高原多年冻土区，得到1982—2018年期间的平均年蒸散发分布(图3(a))，并与已有数据集(青藏高原8 km蒸散发数据集)进行对比(图2(b))。结果表明，采用耦合了活动层厚度的Budyko-Fu模型模拟得到的平均年蒸散发与8 km数据集的ERMS=154.3 mm、R=0.74，模拟效果较好。

3.2 青藏高原多年冻土区蒸散发时空变化特征

本研究结果表明，1982—2018年间青藏高原多年冻土区平均年蒸散发为275.6 mm，空间分布不均(图3(a))，大部分地区的平均年蒸散发在100～600 mm之间，整体上分布规律与降水大致相同，由东南地区(接近1 500 mm)向西北地区(100 mm)递减。青藏高原多年冻土区年蒸散发的变化趋势时空分布如图3(b)所示。结果表明，青藏高原多年冻土区年蒸散发呈显著上升趋势(3.57 mm/a，p<0.05)，从1982年的227.8 mm增加到2018年的307.5 mm。在空间分布上，大部分地区(约80%面积)的年蒸散发呈显著的增加趋势(p<0.05)，东南部局部地区的年蒸散发呈降低趋势，这与Cui等[39]研究得到的空间变化趋势大致相同。

图3 青藏高原1982—2018年平均蒸散发及其变化趋势空间分布Fig.3 Spatial distributions of the annual average evapotranspiration and its changing rate during 1982—2018 on the Qinghai-Tibet Plateau

总体上，本模型在青藏高原多年冻土区的蒸散发及其变化趋势拟合效果较好，与其他相关研究结果相符。例如，Zhang等[40]计算得到青藏高原1966—2000年的平均年蒸散发为248.4 mm，研究期内年蒸散发以0.7 mm/a的速率显著增加；Yin等[41]研究结果表明，1981—2010年间青藏高原平均年蒸散发为255.8 mm；Song等[42]估算出青藏高原2000—2010年平均年蒸散发为350.3 mm；Wang等[15]估计了青藏高原蒸散发的时空格局，指出青藏高原1961—2014年平均年蒸散发为294.2 mm，年蒸散发增长速率为0.38 mm/a；Feng等[43]得出青藏高原2003—2014年间的平均年蒸散发为380.0 mm，年蒸散发以1.2 mm/a的速率增加。此外，这些研究都表明青藏高原的蒸散发空间分布不均，整体上从东南向西北递减。尽管与前人研究的区域和时期不完全一致，本研究模拟的蒸散发及其变化趋势介于前人的估算范围内。

为了探究引起青藏高原多年冻土区蒸散发增加且空间分布不均的控制因素，本文分析了Budyko-Fu模型中的2个控制因素(降水、潜在蒸散发)的变化趋势，如图4(a)和图4(b)所示。结果表明，青藏高原多年冻土区的降水和潜在蒸散发在1982—2018年间整体上均呈现明显的上升趋势，其中，降水的上升速率为4.26 mm/a，潜在蒸散发的上升速率为2.85 mm/a。在空间上，降水变化趋势的分布与蒸散发变化趋势的分布(图3(b))呈现高度的一致性，且在区域平均值的年际变化(图4(c))上步调一致；而潜在蒸散发变化趋势分布与蒸散发变化趋势分布之间的没有明显的相关关系。因此，相较潜在蒸散发(需求)而言，降水(供应)是青藏高原蒸散发过程的主要控制因素，这与青藏高原整体处于干旱、半干旱区的情况相吻合。

图4 研究区降水和潜在蒸散发变化趋势空间分布及区域平均降水和蒸散发的年际变化Fig.4 Spatial distributions of the precipitation and potential evapotranspiration，and interannual evolutions of regional average precipitation and evapotranspiration of the study area

3.3 活动层厚度变化对蒸散发的影响

为了探究冻土退化对高原蒸散发的影响，比较了对照情景和假设情景的2种不同情景下模拟的1982—2018年平均年蒸散发(图5(a)，图中蒸散发变化率为对照情景与假设情景的年蒸散发之差与对照情景年蒸散发的比值)，图5(b)是2种情景下多年平均活动层厚度之差的空间分布(对照情景减去假设情景)。结果表明，由气候变化引起的青藏高原广泛的活动层厚度加深(活动层厚度之差>0)会引起蒸散发的增大(蒸散发变化率>0)，在区域平均活动层厚度增大0.27 m的情景下，平均年蒸散发会增大2.2%(平均年蒸散发变化量为6.01 mm)。换言之，若不考虑气候变化导致的活动层厚度的加深(平均加深11.5%)，会导致模拟的青藏高原多年冻土区平均年蒸散发比真实值低估2.2%。虽然2.2%的低估在整个区域来看相对较小，但是在例如三江源等区域，蒸散发变化率可达6%，对局部水文循环过程模拟具有关键影响。

图5 2种情景下蒸散发变化率、活动层厚度之差及二者关系Fig.5 Changing rate of the annual evapotranspiration and change in active layer thickness between control and hypothetical scenario，and their relationship

在空间分布上，对比分析图5(a)和图5(b)，2种情景对比的蒸散发变化率并不随着活动层厚度之差的增大而增大。以青藏高原东南部为例，该地区的活动层厚度之差为0.3～0.6 m，对应的蒸散发变化率只有1%～2%，蒸散发对冻土退化的响应极小；而高原东部黄河源区活动层厚度之差只有0.1～0.2 m，对应的蒸散发变化率却达4%～5%，蒸散发对冻土退化的响应剧烈。为了探究引起这种情况的原因，进一步分析了青藏高原各点处蒸散发变化率与活动层厚度之差的关系，如图5(c)所示，不同颜色表示对照情景中不同的模型参数ω大小。结果表明：① 整体上活动层厚度的增大会导致更大的蒸散发，但两者变化幅度不同，与对照情景模型参数ω有关系。② 在相同的活动层厚度变化的情况下，对照情景模型参数ω越大，蒸散发的变化率越小，即蒸散发对活动层厚度变化的响应越小；反之，蒸散发对活动层厚度变化的响应越敏感。③ 在对照情景模型参数ω相近(颜色相同)情况下，蒸散发变化率随着活动层厚度变化的增大而增大，即蒸散发的响应随活动层厚度的增大而更剧烈。由于模型的参数ω越大，表示该处土壤有效含水量和植被覆盖度越高，越有利于水分和能量向大气传输，故可知，在土壤有效含水量和植被覆盖度较低的地区，蒸散发对活动层厚度变化的响应更加敏感。

3.4 模型局限性

本文通过改进Budyko-Fu模型的参数化方案，总体更为合理地分析了青藏高原多年冻土区蒸散发的时空变化特征，量化了青藏高原冻土退化对蒸散发过程的影响，但同时也存在一定的局限性：① 本研究在进行参数设置时，考虑到研究区植被类型以草甸、草原等草本植物为主，根系深度相对较浅，因此活动层的加深会促进植被根系向下延展，故假设研究区内植被根系随活动层厚度的增加而加深。但由于青藏高原环境复杂，缺乏实测的根系深度数据资料支撑该假设，加之森林等木本植物的根系较深，可能不会受到活动层加深的影响，因此该假设可能会对结果带来一定程度的不确定性。② 冻土退化的表现特征很多，如地温升高、活动层厚度加深、地下冰融化等[34， 44-45]，这些变化均会对蒸散发过程造成影响，如地温升高会导致模型中需求边界的提升，地下冰融化为液态水的能量消耗将会导致蒸散发被高估[15]，而Budyko-Fu模型无法考虑这部分能量消耗，只考虑了活动层厚度加深对蒸散发的影响，因此为模型带来不确定性。因此，如何更全面地考虑冻土退化对水文过程的影响仍是高寒区水文循环研究的难点，有待进一步研究。

4 结 论

本文基于Budyko-Fu模型，将多年冻土活动层厚度因素耦合入模型参数ω中，构建了适用于多年冻土地区的水热耦合模型，研究了青藏高原多年冻土地区1982—2018年间蒸散发的时空演变特征，并通过情景设置分析了活动层厚度变化对蒸散发的影响，主要结论如下：

(1) 提出的耦合活动层厚度的Budyko-Fu模型可以较好地模拟青藏高原多年冻土区的蒸散发(相关系数为0.84，均方根误差为87.1 mm)，探讨了活动层增厚对蒸散发的影响。

(2) 青藏高原多年冻土区1982—2018年间的平均年蒸散发为275.6 mm，且空间分布与降水分布相似，均为由西北向东南递增；研究期内年蒸散发整体上以3.57 mm/a(p<0.05)的速率显著增加。空间上，80%的地区蒸散发呈显著上升趋势，而仅在东南部降水下降的地区蒸散发呈下降趋势。

(3) 活动层厚度的增加会导致蒸散发的增大，不考虑冻土退化会使得蒸散发被低估2.2%。冻土退化对蒸散发的影响存在空间差异性，土壤有效含水量和植被覆盖度较低的地区对冻土退化的响应更敏感。

多年冻土退化对蒸散发的影响机制极为复杂，未来需建立综合考虑冻土退化引起的地表水文过程、能量平衡及植被功能变化的水热耦合模型以进一步揭示蒸散发对多年冻土退化的响应机制。