基于Budyko假设的祁连山排露沟流域蒸散发估算及其变化归因分析

Estimation and attribution analysis of evapotranspiration variation in Pailugou watershed of Qilian Mountains based on Budyko hypothesis

摘要

目的 探讨西北干旱区域水资源短缺和森林植被碳汇问题。 方法 以祁连山排露沟流域蒸发皿、气象站、水文站、固定样地等实测数据为基础，采用线性回归法、Mann-Kendall、滑动t、Pettitt、累积距平检验方法，Budyko 假设模型、材积源生物量等方法以及Matlab编程软件对实际蒸散发量进行估算，且对其变化趋势及归因进行研究。 结果 1994~2020 年排露沟流域蒸散发总体呈增加趋势，但在2001年有次突变，突变前后变异系数由0.13增大到0.22，属于较强变异；在2002~2020 年的突变期，蒸散发量对降水、径流、潜在蒸发和下垫面参数的弹性系数分别为0.69、-0.46、-0.47、0.62，且各因子对流域蒸散发变化的贡献率分别为71.65%、-85.57%、69.61%、44.31%，降水量变化对流域蒸散发最为敏感，流域下垫条件变化次之。2003~2021年排露沟流域青海云杉平均碳密度为65.8 mg/hm²，增长速率为1.99 mg/（hm²·a），与流域蒸散发的变化趋势基本相一致。祁连山排露沟流域蒸散发的变化不仅受气候变化的影响，同时在短时期内受到青海云杉生长变化的影响，其变化是气候与植被生长综合影响的结果。 结论 研究结果可对全球气候变暖下的水资源短缺与森林植被的碳汇功能降低的研究提供了理论依据。

Abstract

Objective The study aimed to investigate water resource scarcity and carbon sequestration within forest vegetation in the arid regions of northwest China. Method Utilizing data from potential evaporation measurements in evaporating dishes within the Pailougou watershed of the Qilian Mountains，along with data from meteorological stations，hydrological stations，and fixed sample plots of Picea crassifolia within the woodland，the actual evapotranspiration and its changing trends and attribution were determined by employed various methods，such as linear regression analysis，the Mann-Kendall test，the moving t test，the Pettitt test，the cumulative distance test，the Budyko hypothesis model，biomass calculations，MATLAB，Maple software programming，and other analytical tools. Result Between 1994 and 2020，evapotranspiration within the Pailugou watershed exhibited an overall increasing trend，with a notable abrupt change occurring in 2001.The coefficient of variation increased from 0.13 to 0.22 post the abrupt shift，indicating a significant variation.During the period of abrupt change from 2002 to 2020，the elastic coefficients of evapotranspiration in the Deeping Ditch Basin concerning precipitation，runoff，potential evaporation，and underlying surface parameters were 0.69，-0.46，-0.47，and 0.62，respectively.The contribution ratios of each factor to evapotranspiration in the basin were 71.65%，-85.57%，69.61%，and 44.31%，with precipitation exhibiting the highest sensitivity to evapotranspiration changes，followed by alterations in the basin's underlying conditions.The average carbon density of Picea crassifolia in the Pailugou Basin was 65.8 mg/hm² between 2003 and 2021，with an average annual growth rate of mg/(hm²·a)，aligning closely with the evapotranspiration trends in the watershed.Changes in evapotranspiration within the Pailugou watershed of the Qilian Mountains are not solely influenced by climate variations but also by rapid growth changes in Picea crassifolia.These changes reflect the combined impact of climate dynamics and vegetation growth. Conclusion The findings offered crucial insights and served as a foundational reference for further investigations into the issues of water resource scarcity and forest vegetation carbon sequestration under the influences of global climate change.

Graphical abstract

关键词

流域蒸散发 / 归因分析 / 气候变化 / 青海云杉林碳密度 / 祁连山排露沟流域

Key words

watershed evapotranspiration / attribution analysis / climate change / carbon density of picea crassifolia forest / Pailugou watershed of Qilian Mountains

Author summay

金铭，博士，研究员，主要从事生态水文方面的研究。E-mail：562988898@ qq.com

在全球气候变暖和人类活动的双重驱动下，蒸散发、水汽输送、降水、下渗、径流等水循环各个环节均发生了不同程度的改变^［1-2］，进而导致陆地与海洋间的大循环以及陆地、海洋各自内部的小循环过程随之发生着一系列的变化，由于蒸散发是水循环从液态变为气态的相变和潜热变化，对全球气候变暖影响指示更加敏感。另外，由于植物光合作用的碳汇过程与呼吸作用的蒸腾过程对诸多环境因子的响应往往具有一致的规律特征^［3-4］，蒸散发是探索碳循环与水循环联系的重要指标之一。因此，探索蒸散发的动态变化特征及成因机制，对于全球气候变暖和人类活动影响下水资源短缺与森林植被碳汇的水碳循环问题、水资源管理与“碳中和”任务权益平衡等问题的研究均具有重要的现实意义。

流域蒸散发是指一个流域内各种下垫面条件下的总散发，是流域水热循环过程的重要组成和直接环节，也是流域碳循环和水循环的重要“纽带”^［5］，但难以直接观测，通常利用潜在蒸发量（Potential evapotranspiration，PET）数据结合流域水量平衡原理推算流域蒸散发（Basin evapotranspiration，BET）。目前，国内外学者经常利用气象因子对PET的敏感性、通径、回归等进行研究^［6-8］，由于PET与气象因子是非线性关系，变化复杂且存在较大的不确定性。有一部分学者认为PET呈下降趋势^［9-10］，但另一部分学者认为近年来PET出现上升趋势^［11］，其实不同区域的PET对各气象要素的敏感度不同^［12］，海南岛PET对最低和最高温度最敏感^［13］，辽宁省PET对平均风速最敏感^［14］，印度干旱区PET对温度和净辐射极为敏感^［15］，伊朗北部的PET受风速和相对湿度的影响较大^［16］；洮河流域PET对净辐射敏感性最高，气温和净辐射的变化是该流域PET增加的主要原因^［17］。通过分析这些研究成果发现，目前利用气象要素或者基于Penman-Monteith （P-M） 模型方法对PET变化趋势和归因问题进行分析是比较主流的研究方法，并且已取得了很多相关研究成果，但利用Budyko 假设水热耦合平衡模型对基于实测数据估算BET和变化趋势以及归因问题的相关研究极少^［18-20］。基于此，本文中在敏感性分析的基础上，引入贡献率的概念，通过定量分析，进一步理解降水、径流、蒸散发能力、流域下垫面参数对流域蒸散发变化的贡献和响应。本文在前人的基础上，以潜在蒸发、径流、固定样地每木检尺等实测数据，利用M-K、滑动t、Pettitt、累积距平检验及敏感性分析法对祁连山排露沟流域蒸散发演变从时间尺度上进行研究，采用Budyko模型对流域蒸散发变化进行归因分析，并讨论了气候变化和青海云杉生长对流域蒸散发变化的影响。

1 材料与方法

1.1 研究区概况

排露沟流域（E 100°17'~100°18'，N 38°32'~38°33'）位于祁连山北麓中段的黑河上游高山区，属典型的大陆性半干旱高山气候，受西风环流和极地冷气团的影响有明显的干湿季之分。流域总面积2.73 km²，海拔2 640~3 800 m，青海云杉（Picea crassifolia）是该流域的主要建群树种。沿海拔梯度设置其固定样地（图1）。在流域出口处建有水文站（图1），利用雷达水位计和流速仪、浮标相互配合检测径流量。在流域下游不远处有国家三级标准的气象站，站内建有人工蒸发皿观测仪器（图1）。

1.2 数据来源

1994~2020年连续日潜在蒸发量数据（PET）来源于排露沟流域20 cm口径的蒸发皿，每日人工观测数据、降水（P）及气温数据（T）来源于祁连山排露流域自动气象站（Campbell-CR1000）；径流量数据（R）来源于流域出口处的水文监测站；青海云杉胸径和树高数据来源于2003~2021年排露沟流域每木检尺数据。对收集到的资料进行分类汇总，并运用excel、Spss13.0、Origin 10.0、matlab等软件进行统计分析。

1.3 研究方法

1.3.1 BET和流域下垫面参数估算

一般认为流域多年蓄水量变化可忽略不计，根据流域水量平衡理论和蒸散发量公式BET=P-R，利用实测径流深R（mm）、实测降水量P（mm）求得BET（mm）；然后基于 Budyko 假设理论和模型BET=\frac{P\times PET}{{\left(P^{\omega }+PE{T}^{\omega }\right)}^{\frac{\mathrm{1}}{\omega }}}，在MATLAB软件拟合流域下垫面参数\omega。

1.3.2 BET的变化趋势及突变检验

采用排露沟流域1994~2020年所求的BET数据，利用线性回归和5 a 滑动的方法进行BET变化趋势分析，利用Mann-Kendall、滑动、Pettitt、累积距平4种方法对BET突变检验，6种方法综合研判，寻找BET序列的突变时间点^［21-23］。

Mann-Kendall 检验法：在时间序列BET上，将第 i 个样本 BE{T}_i > BE{T}_j 的累计数统计成为：

BE{T}_k={\sum }_{i=\mathrm{1}}^{k}BE{T}_i(\mathrm{1}\le k\le n)

式中：BE{T}_i为前后两个样本BET的值分别取0 和 1。

1.3.2.1 滑动t检验法

假设BET存在一个突变点，突变点前后R容量分别为n_{\mathrm{1}}、n_{\mathrm{2}}

\mathrm{\sigma }=\frac{|{\overline{BET}}_{\mathrm{1}}-{\overline{BET}}_{\mathrm{2}}|}{\sqrt[]{\frac{n_{\mathrm{1}}{\mathrm{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+n_{\mathrm{2}}{\mathrm{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}{n_{\mathrm{1}}+n_{\mathrm{2}}-\mathrm{2}}(\frac{\mathrm{1}}{n_{\mathrm{1}}}+\frac{\mathrm{1}}{n_{\mathrm{2}}})}}

式中：\mathrm{\sigma }服从自由度为n_{\mathrm{1}}+n_{\mathrm{2}}-2的t分布；n_{\mathrm{1}}、n_{\mathrm{2}}分别为某一基准年前后两段BET值个数；{\overline{BET}}_{\mathrm{1}}、{\overline{BET}}_{\mathrm{2}} 分别为某一基准年前后两段BET相邻时段斜率平均值；{\mathrm{\sigma }}_{\mathrm{1}}、{\mathrm{\sigma }}_{\mathrm{2}}分别为前后两段样本标准差。

1.3.2.2 Pettitt法

假设BET存在突变点m，突变时刻t_p，采用 Mann-Whitney 的统计量U_{t_p}，n检验BET突变点m前后两个子样本｛［BET1］、（BET2）、…［BETtp］｝和｛［BETtp］+1、［BETtp］+2、…［BETn］｝二者间累积分布是否存在显著差异。统计量U_{t_{{\mathrm{p}}^{,n}}}计算公式为：

U_{t_{{\mathrm{p}}^{,n}}}=U_{t_{{\mathrm{p}}^{-\mathrm{1},n}}}+{\sum }_{\mathrm{i}=\mathrm{1}}^n\mathrm{s}\mathrm{g}\mathrm{n}\left(BE{T}_{t_p}-BE{T}_{\mathrm{i}}\right)({\mathrm{t}}_{\mathrm{p}}=\mathrm{2,3},\cdots ,n)

K={}_{\mathrm{1}\le {\mathrm{t}}_{\mathrm{p}}\le n}{}^{\mathrm{m}\mathrm{a}\mathrm{x}}|U_{t_p}|,P\cong {\mathrm{2}}_{\mathrm{e}\mathrm{x}\mathrm{p}}[-\mathrm{6}{\mathrm{k}}_{\mathrm{t}}^{\mathrm{2}}/(n^{\mathrm{3}}+n^{\mathrm{2}})]

式中：BE{T}_{t_p}、BE{T}_i分别为第t_{\mathrm{p}}、i 年BET斜率值；K为 t_{\mathrm{p}} 时刻出现 U_{t_p}，n最大统计量；若概率P≤0.05，则m 点为显著突变点。

1.3.2.3 累积距平法

差积曲线上的拐点即为BET的突变点。

BE{T}_{{\mathrm{t}}_a}={\sum }_{i=\mathrm{1}}^{{\mathrm{t}}_a}(BE{T}_i-\overline{BET})(\mathrm{1}\le {\mathrm{t}}_a\le n)，

式中：BE{T}_{{\mathrm{t}}_a}为BET在时间序列 {\mathrm{t}}_a 时段的累积距平值；BE{T}_i 为BET在 i时刻斜率值；\overline{BET}为BET序列平均值；n为BET的年际区间。

1.3.3 BET变化归因分析

1.3.3.1 BET敏感度分析法

基于 Budyko 假设水热耦合平衡原理^［20］，假设 P、R、PET、\omega为相互独立的变量x，用弹性系数 ε _{x表示变量x变化对BET的影响敏感度，ε x 含义为x变化1个单位所引起的BET变化量。利用MATLAB和Maple软件编程求得P、R、PET、\omega的偏导数和弹性系数ε x。}

1.3.3.2 BET贡献量分析法

dBE{T}^{\text{'}}=dBE{T}_p+dBE{T}_R+dBE{T}_n+dBE{T}_{\omega }

={\epsilon }_p\frac{BET}{P}dP+{\epsilon }_R\frac{BET}{R}dR+{\epsilon }_{PET}\frac{BET}{PET}dPET+dBE{T}_{\omega }

式中：dBET′为Budyko 假设求算的BET突变前后的变化量，dBE{T}_p、dBE{T}_R、dBE{T}_{PET}、dBE{T}_{\omega }为P、R、PET、\omega变化引起的BET变化量，dP、dR、dPET、d\omega分别为 P、R、PET、\omega突变前后多年平均值的差值。

{\mathrm{C}}_{\mathrm{p}}=\frac{dBE{T}_p}{dBE{T}^{\text{'}}}\times \mathrm{100}%

{\mathrm{C}}_{\mathrm{R}}=\frac{dBE{T}_R}{dBE{T}^{\text{'}}}\times \mathrm{100}%

{\mathrm{C}}_{\mathrm{P}\mathrm{E}\mathrm{T}}=\frac{dBE{T}_{PET}}{BdE{T}^{\text{'}}}\times \mathrm{100}%

C_{\omega }=\frac{dBE{T}_p}{dBE{T}^{\text{'}}}\times \mathrm{100}%

式中：{\mathrm{C}}_{\mathrm{p}}、{\mathrm{C}}_{\mathrm{R}}、{\mathrm{C}}_{\mathrm{P}\mathrm{E}\mathrm{T}}、C_{\omega }分别为P、R、PET、\omega变化引起BET变化的贡献率。

1.3.4 青海云杉林碳密度计算

1.3.4.1 固定样地每木检尺方法

采用常规的每木检尺方法，在2003~2021年期间对研究区16个固定样地的1 120棵青海云杉每隔5 a或4 a的6~8月份进行一次固定样地调查，调查内容较多，汲及本文的主要包括青海云杉的树高和胸径，固定样地特征因子如表1所示。

1.3.4.2 青海云杉林碳密度计算方法

本文主要采用林积源生物量法对青海云杉林碳密度进行估算，其表达式简述如下^［24-25］：

V_{ij}=a({D_{ij}}^{\mathrm{2}}{H}_{ij}{)}^b,B_{ij}=c{V}_{ij},C_{ij}=d{B}_{ij}

C_j={\sum }_{i=\mathrm{1}}^n{C}_i，C_{dj}=C_j/A_j

式中：V_{ij}、D_{ij}、H_{ij}、B_{ij}、C_{ij}分别为第j个固定样地内第i棵青海云杉单株蓄积量（m³）、单株胸径（cm）、单株树高（m）、单株生物量（kg）、单株碳储量（kg），a、b分别为青海云杉全株蓄积量经验系数，本文取值分别为0.431 2，0.742 5 ^［26-27］，c、d 分别为青海云杉生物量、碳密度经验系数，本文取值分别为0.593 9、0.462^［26-27］。C_j、n、C_{dj}、A_j分别第j个固定样地内青海云杉生物量（kg）、株数、碳密度（Mg/hm²）、样地面积（hm²）。

2 结果与分析

2.1 流域蒸散发变化趋势分析及突变检验

2.1.1 流域蒸散发变化趋势

利用排露沟流域降水量和径流量的实测数据，采用BET的计算方法，求出1994~2020年BET。如图2所示，BET多年平均为271.69 mm，2017年出现最大值（424.76 mm），2020年出现最小值（170.60 mm）。从线性拟合情况来看，BET变化总体呈波动增大趋势，变率为7.0 mm/10 a。从5 a 滑动平均曲线来看，BET围绕多年平均线呈上下波动的交替变化趋势。

2.1.2 流域蒸散发突变检验

图3-A 为 M-K 突变检验分析图，统计量UF（k）在 1994~2001 年间均为负值且处于 0.05 显著水平临界线（+1.96）内，表明BET在此时间段内呈现波动下降趋势，在2002~2020 年间UF（k）基本为正值且基本处于0.05 显著水平临界线（+1.96）内，表明BET在此时间段内呈波动上升趋势；UB（k）和UF（k）两条曲线有多次相交，但2001年的相交点较为理想；根据M-K突变检验规则可知，BET突变时间最有可能发生在2001年。图3-B为滑动t检验，图中统计量 t 在2001 年超出给定P=0.05显著水平临界线，表明BET在2001年发生了突变且较显著。图3-C为Pettitt检验，根据Pettitt 突变检验中T₀值与UT值的交点原则，表明BET可能在2013年前后发生了一次显著突变。图 3-D为累积距平曲线，由图可以看出BET累积距平曲线，根据曲线最低点原则，可能在2001年、2005、2013年BET发生过突变。

综上分析，虽然M-K、滑动t、Pettitt、累积距平4种方法对BET的突变检验结果不完全一致，但其中3种方法的分析结果均表明BET在2001年前后有突变，再结合图 2 中年BET变化趋势分析，确定2001 年是排露沟流域BET变化的突变年份。

突变前后BET变化情况见表2所示，BET突变后有较大幅度的增加，由突变前的 32.76 mm增加到突变后的60.25 mm，增加10.12%；变异系数由0.13 变到0.22，表明BET突变后的年际间变异度更大。

2.2 流域蒸散发变化归因分析

2.2.1 敏感性分析

2001年作为研究期1994~2020年的突变点，将其分为基准期1994~2001年和变化期2002~2020 年两段（表3）。由表3可见，BET、P、R、ω的变化期比基准期均有所上升，上升速率分别为8.09%、11.74%、20.95%、8.13%，而PET与基准期相比有所下降，下降速率为16.66%。径流系数α较基准期有所增大，表明降水转化为径流的比例在增大；干燥指数k较基准期有较大幅度下降，表明变化期比基准期空气更潮湿。ε_p为正值，表示降雨量P与BET呈正相关，在P、R、PET、ω不变的条件下，当P每增加1 mm 时，BET平均增加0.72 mm，但ε_p从基准期的0.72降至变化期的0.69，表明正相关呈减弱趋势，这可能与R、PET、ω指标的增强有关。依次类推，当其他影响因子不变的情况下，当R、PET、ω中任何一个影响因子变化1 mm或1时，BET的变化量由基准期的-0.40 mm、-0.26 mm、0.61变到变化期的-0.46 mm、-0.47 mm、0.62（表3），弹性系数ε的绝对值又代表BET受影响的敏感度^［17］，结合表3可知，研究期内BET对P的敏感度最强，ω次之，R、PET的敏感度较弱。

2.2.2 归因分析

由表 4 可知，利用Budyko模型原理计算BET变化量dR′为 29.01 mm，实际BET变化量dR值为 20.79 mm。变化期内P、R、PET、ω对BET的贡献程度差异较大，气候因素（P、R、PET ）对BET的贡献率为 55.69%，其中P、R、PET引起的BET变化的贡献率分别为71.65%、-85.57%、69.617%；流域下垫面因素ω引起的BET变化量为12.85 mm，贡献率为44.31%。

综上可知，BET变化主要归因于气候因素变化的驱动作用，但P、PET的影响是正向的，而R的影响是负向，气候因素对BET的综合影响存在着正负作用相互制衡的作用，ω对BET的影响也是正向促进的作用。

2.3 气候变化对流域蒸散发的影响

气候因子是描述气候变化的重要指标，不仅包括P、R、PET等指标，气温（T_a）也是重要的描述指标。T_a是BET的主要热动力，BET的变化受到P、R、PET、T_a等因子的综合影响^［22］，为了进一步探讨BET变化的驱动作用，将P、R、PET、T_a与BET进行散点相关分析（图4）。

如图4-A所示，在基准期BET与Ta呈现弱负相关，递减率为-5.651 5；在变化期BET与T_a呈弱正相关，增长率为30.77，说明BET从基准期到变化期，T_a对其引起变化的贡献在增强。由图-B可知，BET与P呈较强正相关，BET从基准期到变化期随P变化的增长率从0.369 4增长至0.717 1，说明变化期P对BET的影响程度也逐渐显著。由图4-C可知，BET与R呈较弱正相关，BET从基准期到变化期随R变化的增长率从0.066 6增长至0.447 9，说明变化期R对BET的影响程度也逐渐在增强。由图4-D可知，BET与PET呈较弱负相关，BET从基准期到变化期随PET变化的递减率从-0.050 7增大到-0.159 4，说明变化期PET对BET的影响程度也逐渐在增强。经分析可知，研究区BET从基准期到变化期，T_a、P、R、PET对其变化的影响均在逐渐增强，进一步解释了研究区近27 a以来，BET总体呈现出增长趋势的原因，即BET的增长趋势是受T_a、P、R、PET等因子综合影响的结果，此与相关学者^［21-23］的研究结果基本一致。

2.4 下垫面条件变化对流域蒸散发的影响

2.4.1 下垫面参数对流域蒸散发的影响

由图5可知，ω与BET在基准期呈现较强正相关，从基准期到变化期BET随ω变化的上升率从130.3 增加为265.8，说明变化期ω对BET的影响在增强。

2.4.2 下垫面青海云杉林碳密度对流域蒸散发的影响

有关研究^［28-32］表明，ω与流域地形、土壤、植被覆盖度等因素有关，地形、土壤在较短时间内变化较小，可忽略，因此，在一个人为扰动较小的流域，植被覆盖度是ω的变化主要因子。碳汇功能是反映植被覆盖度变化的重要指标，其变化不仅强烈地依赖BET的变化，而且对BET的变化有较敏感的指示作用；又由于青海云杉是主要建群树种，一定程度上可反映该研究区的植被覆盖度情况。因此，采用研究区2003~2021年固定样地每木检尺数据计算得到的青海云杉碳密度指标，进一步探讨ω对BET的影响作用。

由表5可知，2003~2021 年研究区青海云杉平均碳密度为65.8 mg/hm²，且呈现增长趋势，平均增长速率为1.99 mg/（hm²·a），碳密度呈增长趋势，平均增长速率为2.24 mg/（hm²·a），与BET的变化趋势基本相一致。这充分说明，BET的变化不仅受P、R、PET、T_a等气候因子的影响，而且同时也受到ω变化的影响。

3 结论

1） 1994~2020 年排露沟流域蒸散发总体呈增加趋势，但在2001年有次突变，突变前后变异系数由0.13增大到0.22，属于较强变异；蒸散发对降水、径流、潜在蒸发和下垫面参数的弹性系数分别为0.69、-0.46、-0.47、0.62，蒸散发对降水的敏感度最强，对下垫面条件次之，对径流、潜在蒸发的敏感度较弱；降水、径流、潜在蒸发和下垫面参数对蒸散发的贡献率分别为71.65%、-85.57%、69.61%、44.31%，降水、潜在蒸发、下垫面参数对流域蒸散发变化贡献是正向的，而且依次减弱，径流对蒸散发变化贡献是较强负向。

2） 流域蒸散发从基准期到变化期，包括气温、降水、径流、潜在蒸发等在内的气候因子对其变化的影响均呈逐渐增强趋势。

3） 2003~2021年排露沟流域青海云杉平均碳密度为65.8 mg/hm²，平均增长速率为1.99 mg/（hm²·a），与流域蒸散发的变化趋势基本相一致；在人类活动干扰较小的条件下，森林植被的碳密度变化是探讨流域下垫面参数变化的重要指标，祁连山排露沟流域蒸散发的变化不仅受气候变化的影响，同时在短时期内受到青海云杉生长变化的影响，其变化是气候与植被生长综合影响的结果。

参考文献

[1] Oki T， Kanae S.Global hydrological cycles and world water resources［J］.Science，2006，313（5790）： 1068-1072.

[2] 宋晓猛，张建云，占车生，等.气候变化和人类活动对水文循环影响研究进展［J］.水利学报，2013，44（7）：779-790.

[3] Jia X， Zha T， Gong J，et al.Carbon and water exchange over a temperate semi -arid shrubland during three years of contrasting precipitation and soil moisture patterns［J］.Agricultural and Forest Meteorology，2016，228：120-129.

[4] Chandrasekaran M， Chanratana M， Kim K，et al.Impact of arbuscular mycorrhizal fungi on photosynthesis，water status，and gas exchange of plants under salt stress-a meta-analysis［J］.Frontiers in plant science，2019，10，457.

[5] Chahine M T.The hydrological cycle and its influence on climate［J］.Nature，1992，359（6394）：373-380.

[6] 尹云鹤，吴绍洪，戴尔阜.1971~2008 年我国潜在蒸散时空演变的归因［J］.科学通报，2010，55（22）：2226-2234

[7] 张薇，韦群，吴天傲，等.基于GBDT 算法的参考作物蒸散量模型在江苏省的预测［J］.江苏农业学报，2020，36（5）：1169-1180.

[8] 海涛，严新军，侍克斌，等.吐鲁番地区水面蒸发演化特征及影响因素［J］.甘肃农业大学学报，2022，57（2）：137-144.

[9] Liu C M， Zeng Y.Changes of pan evaporation in the recent 40 years in the Yellow River Basin［J］.Water International，2004，29（4）：510 -516.

[10] Jhajharia D， Dinpashoh Y， Kahya E，et al.Trends in reference evapotranspiration in the humid region of northeast India［J］.Hydrological Processes，2012，26（3）：421-435.

[11] Cong Z T， Zhao J J， Yang D W，et al.Understanding the hydrological trends of river basins in China［J］.Journal of Hydrology，2010，388（3/4）：350-356.

[12] 刘昌明，张丹.中国地表潜在蒸散发敏感性的时空变化特征分析［J］.地理学报，2011，66（5）：579-588

[13] 邹海平，陈汇林，田光辉，等.海南岛参考作物蒸散量时空变化特征及成因分析［J］.中国农业气象，2018，39（1）：18-26.

[14] 曹永强，高璐，袁立婷，等.辽宁省潜在蒸散发量及其敏感性规律分析［J］.地理科学，2017，37（9）：1422-1429.

[15] Goyal R K.Sensitivity of evapotranspiration to global warming： a case study of arid zone of Rajasthan（India）［J］.Agricultural Water Management，2004，69（1）：1-11.

[16] Dinpashoh Y， Jhajharia D， Fakheri-fard A，et al.Trends in reference crop evapotranspiration over Iran［J］.Journal of Hydrology，2011，399（3/4）： 422-433.

[17] 杨林山，李常斌，王帅兵，等.洮河流域潜在蒸散发的气候敏感性分析［J］.农业工程学报，2014，30（11）：102-109.

[18] 王亚琴，杨巍，邢博，等.Budyko框架下黄河流域蒸散发时空变化影响因素研究［J］.水文地质工程地质，2023，50（2）：23-33.

[19] 蒋小满，王国强，薛宝林，等.基于Budyko理论和互补模型的浑河流域实际蒸散发时空变化特征及驱动力分析［J］.生态学杂志，2022，41（7）：1361-1368.

[20] 余加男，李占玲，冯雅茹.基于Budyko理论的黑河流域中游地区实际蒸散发估算及其变化归因分析［J］.节水灌溉，2022（2）：54-58.

[21] 祁添垚，张强，王月，等.1960~2005 年中国蒸发皿蒸发量变化趋势及其影响因素分析［J］.地理科学，2015，35（12）：1599-1606.

[22] 赵宝山，严程明，苏俊波，等.1960~2020年海南岛气温、降水及参考作物蒸散量变化趋势［J］.节水灌溉，2023（10）：83-90.

[23] 文建川，梁平，廖育玲，等.1971~2020年贵州从江县参考作物蒸散量变化特征及成因分析［J］.中国农学通报，2023，39（17）：72-78.

[24] 王秀云，孙玉军.森林生态系统碳储量估测方法及其研究进展［J］.世界林业研究，2008，21（5）：24-29.

[25] 牛赟，毛广雄.林区资源与环境调查教程［M］.兰州：兰州大学出版社，2021.

[26] 彭守璋，赵传燕，郑祥霖，等.祁连山青海云杉林生物量和碳储量空间分布特征［J］.应用生态学报，2011，22（7）：1689-1694.

[27] 刘建泉，李进军，郝虎，等.祁连山青海云杉林生物量与碳储量及其影响因素分析［J］.现代农业科技，2017（12）：140-143+146.

[28] 祁晓凡，李文鹏，崔虎群，等.黑河流域中游盆地地表水与地下水转化机制研究［J］.水文地质工程地质，2022，49（3）：29-43.

[29] 何旭强，张勃，孙力炜，等.气候变化和人类活动对黑河上中游径流量变化的贡献率［J］.生态学杂志，2012，31（11）：2884-2890.

[30] 潘燕辉，张辉，马金珠.气候变化和人类活动对黑河水资源的影响［J］.人民黄河，2012，34（5）：55-60.

[31] 李旭梅，杨钧，李小伟，等.宁夏天然草地土壤有机碳分布格局及驱动因素［J］.草原与草坪，2024，44（4）：1-8.

基金资助

甘肃省自然科学基金重点项目(22JR5RG1029)

甘肃省中央引导地方科技发展资金项目(22ZY2QG001)

国家自然科学基金项目(32060247)

江苏高等学校协同创新计划项目(HSXT30511)

河西学院科研项目(KYQD2023014)