水循环和能量循环之间的相互作用研究是世界气候研究计划(WCPR)4个核心项目之一。径流是水循环的关键环节，深受气候变化的影响，同时其变化直接影响可利用水资源总量。径流变化的归因分析，尤其是气候因素对径流变化的影响一直是研究的重点。其中，降水被广泛认为是影响径流变化的主要气候因素。例如，Liu等指出降水对1966—2013年长江源区径流变化的贡献达99.7%。潜在蒸发对径流变化也有显著的调节作用，其增加促使径流减少，而其减小则使径流增加。已有研究表明，潜在蒸发增加对1966—2013年长江源区径流增加的贡献约为-6.8%，而潜在蒸发减少对1960—2010年阿克苏河径流增加的贡献为27%左右。同时，地表覆被也是影响径流的重要因素。例如，Buendia等的研究表明，1965—2009年森林扩张对诺格拉帕利亚雷萨河径流减少的贡献达40%。此外，在全球变暖的背景下，冰川、积雪融水以及多年冻土对径流的调节作用逐渐增强，其中冰川融水可占1961—2006年塔里木河四源流(叶尔羌河、阿克苏河、和田河和开都河)河流径流增加量的85.7%左右。

总体而言，学术界利用现代观测资料在径流与气候要素变化之间定量关系方面已取得了较为深入的系统性认识。但是，现代观测资料时间较短，现有研究大多是针对年际—年代际尺度径流量变化的归因分析。古气候研究发现气候要素(温度、降水等)呈现多尺度特征，不仅具有年际变化，而且存在鲜明的百年尺度变化。IPCC评估报告指出，在百年尺度上，过去千年存在中世纪气候异常期(MCA)、小冰期(LIA)和现代暖期(MWP)3个暖—冷—暖气候特征时期。然而，目前对百年尺度上径流变化的成因，特别是其与气候要素变化的关系尚不清楚，有待进一步研究。

黄河流域历史悠久，在中国经济社会中具有重要的地位。因其主体属于半干旱—干旱气候区，水资源短缺是该地区突出的环境问题。同时，近年来东亚季风增强导致极端降水增多，洪涝灾害频发。径流作为水量平衡的基础要素，研究其百年尺度变化，对于探究径流变化的多尺度特征及应对气候变化调节水资源配置具有重要意义。因此，本文以黄河中、上游为研究区，重点研究MCA、LIA和MWP三个气候特征期之间径流的差异，厘清降水、潜在蒸发(大气蒸发力)、地表覆被因素对径流变化的贡献，以期增进对黄河流域径流量多尺度变化及其成因的认识。

1．研究区域、资料与方法

1.1研究区概况

黄河流域位于96~119ºE、32~42ºN之间，发源于青藏高原巴颜喀拉山脉，流经青海、四川、甘肃、宁夏、内蒙古、陕西、山西、山东和河南9个省(自治区)的370多个县市，汇入渤海，全长5464km，流域总面积达79.5万km2。流域至西向东地势起伏剧烈，气候差异显著，包括干旱、半干旱、半湿润等气候类型，降水时空分配不均，多年平均降水量由东南向西北递减，流域内大部分地区年降水量介于200~650mm。根据流域结构，将黄河流域分为上游、中游和下游3个研究区域。黄河上游从源头到内蒙古河口镇，全长3472km，流域面积42.8万km2；中游从河口镇到河南郑州桃花峪，全长1206km，流域面积34.4万km2；下游从桃花峪到渤海，全长786km，流域面积2.3万km2(图1)。

1.2研究资料

现代仪器观测的气象与水文资料一般不超过百年，大多数国家的现代观测始于20世纪中叶，因而难以支撑多年代际到百年际尺度的气象与水文变化规律的研究。地球系统模式以大气、海洋、陆面和海冰各分量模型为主要组成部分，同时考虑了圈层之间复杂的相互作用，是理解过去气候与环境演变机理、预估未来全球变化情景的重要工具。

模式模拟的时间尺度可达千年，长时间模拟能够再现多尺度特征，特别是能弥补现代观测资料时间短而导致的不足以揭示多年代际—百年际尺度变化规律的问题。本文收集了第3次古气候模拟比较计划(PMIP3)过去千年气候模拟试验与第5次耦合模式评估比较计划(CMIP5)全强迫历史模拟试验的8个地球系统模式试验结果(模式详细介绍见表1)。

其中，过去千年试验的积分时段为850—1850年，历史模拟试验积分时段为1850—2005年。本文采用了模式输出的月尺度水文气象要素，包括降水、径流、蒸散发、地表气压、温度、风速、叶面积指数、相对湿度、地表潜热通量以及地表感热通量。鉴于模式中大多没有汇流模块，因而缺失水文学意义上的河道径流，本研究采用流域内所有网格平均径流深(即水文学意义上的“产流”)作为代理变量，下文简称“径流”。

现代径流系数，由径流与降水的比值而得，其中降水数据来源于中国气象数据网(http://data.cma.cn)提供的1979—2000年中国地面降水月值0.5°×0.5°格点数据集，径流则来源于国家地球系统科学数据中心(http://www.geodata.cn)提供的1979—2000年中国区域10km分辨率逐日总径流数据集。

1.3研究方法

首先，通过对模式模拟径流系数与观测资料进行对比评估各模式的精度。径流系数是表征气候水文要素和流域下垫面对降水形成径流影响程度的综合指标，具有明显的空间分异，对比模拟与观测的径流系数空间分异特征，可用于指示模式对气候、水文要素的综合模拟能力。因此，本文以根据观测数据计算的1979—2000年径流系数作为参照，将8个模式计算的同时段径流系数与其进行对比，评价模拟结果的误差大小，并绘制8个模式模拟的黄河流域和中国区域的泰勒图，进一步评估模式对径流系数的空间模拟能力，遴选出最优的模式。

进而，基于最优模式模拟的温度变化明确过去千年3个气候特征期(MCA、LIA和MWP)的起止时间，对比3个时期径流的差异。再基于Budyko假说与傅抱璞公式计算径流对各因子的弹性系数，并探究各因子对百年尺度径流变化的贡献。弹性系数是用于度量一个变量的变化幅度对另一个变量变化幅度依存敏感度的指标，被广泛应用于量化径流与降水、蒸散发以及地表特征参数的依存敏感度。计算步骤如下：

根据Budyko假说，多年平均的实际蒸散发量决定于大气对陆地的补给量(降水量)与蒸发能力(潜在蒸散量)之间的平衡，并以其作为实际蒸发的限定边界条件给出了相应的理论模型，即式(1)。近年来，在Budyko理论模型的基础上发展出了大量考虑下垫面因素的经验公式，其中包括傅抱璞提出的经验公式(式(2))；进一步，结合流域水量平衡公式(式(3))，式(2)可简化为式(4)。

式中：P为多年平均降水量(mm)；E为多年平均蒸散发量(mm)；R为多年平均径流量(mm)；E0为多年平均潜在蒸发量(mm)；ω为地表参数，与地形、植被及土壤透水性有关；ΔS为流域蓄水变化量(mm)，多年平均ΔS近似于零。为估算各因子对径流变化的贡献，依据式(4)并结合全微分方程(5)与偏导数计算径流对降水、蒸散发、地表参数的弹性系数，公式如下：

式中：ΔR为径流变化量(mm)；ΔP为降水变化量(mm)；ΔE0为潜在蒸发变化量(mm)；Δω为地表参数变化量(mm)；∂R/∂P表示径流对降水的弹性系数；∂R/∂E0表示径流对蒸散发的弹性系数；∂R/∂ω表示径流对流域特征的弹性系数。弹性系数与各变量变化幅度的乘积则表示该变量变化导致的径流变化量。

在计算弹性系数时，实际蒸散发采用模式模拟结果，潜在蒸散发采用世界粮农组织推荐的彭曼公式进行计算：

式中：E0为潜在蒸发量(mm·d-1)；Rn为地表净辐射量(MJ·m-2·d-1)；G为土壤热通量(MJ·m-2·d-1)；Tmean为日平均气温(℃)；u2为2m高处风(m·s-1)；es为饱和水汽压(kPa)；ea为实际水汽压(kPa)；Δ为饱和水汽曲线斜率(kPa·℃-1)；γ为干湿表常数(kPa·℃-1)。

地表参数ω的变化受多种因素影响，因此其对径流的影响也较为复杂，采用傅抱璞提出的半参数公式进行计算：

式中：p为日降水量(mm)；k、g与h为常数；α为径流系数。

2．结果分析

2.1模式精度评估

利用观测数据与8个气候(地球)系统模式模拟数据计算了1979—2000年黄河流域的径流系数，如图2所示。由观测数据计算结果看出(图2a)，黄河流域的径流系数存在鲜明的地域差异，总体呈现南高北低的地域分异特征。流域西北部干旱区的潜在蒸发大，地表产流微弱，导致径流系数极小，而东南地区降水相对充沛易于产流，径流系数较大。

从8个模式的模拟结果与基于观测数据的计算结果对比来看，BCC-CSM1-1(图2b)和MPI-ESM-P(图2e)的模拟结果基本再现了径流系数南高北低的空间分布特征。不过，BCC-CSM1-1和MPI-ESM-P的模拟结果与实测值之间也存在一定差异。其中，BCC-CSM1-1模拟的黄河上游地区的径流系数偏大，这是因为本文采用的径流系数计算方法是径流与同期降水量的比值，然而BCC-CSM1-1模式中冰川消融形成的径流未从总径流中分离出来，模拟输出的径流中包含“历史时期”积累的冰川消融产生的径流，从而导致模拟的径流系数明显大于实测值。MPI-ESM-P模拟的黄河中游东南部的径流系数约0.28，低于观测值。

另外，BCC-CSM1-1和MPI-ESM-P模拟的径流系数之间也存在一定差异，原因主要是各模式的过去千年模拟试验的设置不尽相同。例如：在是否考虑碳循环与采用的植被覆盖数据方面，BCC-CSM1-1的试验考虑了碳循环过程，植被冠层的各类参数实时动态地响应气象条件的变化，指征植被覆盖的叶面积指数、覆盖度、冠层反照率等参数均是动态的。MPI-ESM-P试验则未考虑碳循环过程，而是采用了基于各类历史资料重建的过去千年农田、牧场、天然林、灌木以及草地变化的数据集，受人类农业开发的影响，同一地区的植被类型存在较大的变化，但对于同一类植被而言，其冠层参数不随气候条件变化而变化。鉴于本文重点分析自然条件下径流的变化特征及其成因，因而下文将采用BCC-CSM1-1的模拟结果。

尽管BCC-CSM1-1和MPI-ESM-P模拟结果存在一定的不确定性，相比而言，其余6个模式模拟的径流系数与实测数据计算结果的差异则极为明显。例如，HadCM3(图2c)、IPSL-CM5A-LR(图2d)与CCSM4(图2f)模拟的黄河中游东南部地区的平均径流系数明显偏小，分别约为0.15、0.21、0.22，实际上该地区的径流系数可达0.42左右。而CSIRO-Mk3L-1-2(图2h)模拟的径流系数呈西高东低的空间特征，与实测南北差异的空间分布特征明显不同，且中游东南部径流系数平均值仅为0.19，明显低于实测值。MRI-CGCM3(图2g)与GISS-E2-R(图2i)模拟的黄河流域大部分地区的径流系数均接近于0，无明显空间分布特征，与实测值差异最大。因而，这6个模式的模拟结果并不能体现整个黄河流域径流系数南高北低的空间分布特征。

为了进一步评估8个气候(地球)系统模式对径流系数的空间模拟能力，对其模拟的黄河流域绘制了泰勒图(图3a)。从空间相关系数来看，所有模式在0.46~0.74之间。其中BCC-CSM1-1和MPI-ESM-P模拟的径流系数与观测场相关最高，分别为0.74、0.68，而MRI-CGCM3模拟得到的相关系数与其他几个模式相比较差。所有模式与观测的标准差之比在0.38~0.88之间，其中BCC-CSM1-1与观测的标准差最为接近，而MRI-CGCM3与观测的标准差差异最大，其他几个模式径流系数与观测的标准差之比处于中等水平。

同时，为了评估气候(地球)系统模式对径流系数大空间尺度的模拟能力，对其模拟的中国区域也绘制了泰勒图(图3b)。从空间相关系数来看，所有模式在0.60~0.74之间。其中BCC-CSM1-1和MPI-ESM-P模拟的径流系数与观测场相关系数最高，均为0.74。所有模式与观测的标准差在0.44~0.64之间，其中BCC-CSM1-1和IPSL-CM5A-LR与观测的标准差最为接近，其他几个模式径流系数与观测的标准差之比处于中等水平。总体来看，BCC-CSM1-1为最优模式，其基本能合理模拟出黄河流域以及整个中国区域的径流系数空间分布。

2.2全球温度变化及黄河中上游径流变化

BCC-CSM1-1模拟的过去千年全球温度变化呈现出鲜明的多尺度特征。如图4所示，过去千年温度变化在百年尺度存在明显的暖—冷—暖的阶段性特征。暖、冷、暖3个阶段在时间上与IPCCAR5指出的中世纪气候异常期(MCA，950—1250年)、小冰期(LIA，1450—1850年)以及现代暖期(MWP)大致对应。不过由于模式的误差及基于代用资料重建结果的不确定性，BCC-CSM1-1模拟与重建温度的起止时间存在一定差异。鉴于本文重点是分析百年尺度冷、暖阶段之间地表径流的差异，为突显冷、暖的差异且确保地表要素之间的自洽，本文依据BCC-CSM1-1的模拟结果选择了950—1250年中最暖的200a和1450—1850年中最冷的200a，即991—1191年和1591—1791年，分别代表MCA和LIA，其平均温度分别为13.8℃和13.6℃，选择1875—2012年代表MWP，其平均温度为14.3℃。下文将就这3个时段之间的地表径流差异及成因开展分析。

随气候变化，过去千年黄河中上游径流也呈现出鲜明的多尺度变化特征(图5)。在百年尺度上，中游地区从MCA到LIA径流呈增加的趋势，到了MWP径流明显减少。LIA径流在3个特征时期最大，达560.99mm；MCA径流最低，为530.69mm，二者较850—2012年平均值分别偏高16mm(+2.94%)、偏低14.31mm(-2.63%)。MWP径流介于LIA与MCA之间，为554.42mm，较850—2012年平均值偏高9.42mm(+1.73%)。

与中游地区相比，上游地区径流变率较小，850—2012年的方差仅为中游地区的23%。就3个特征时期而言，MWP径流最大，为270.58mm；LIA径流最低，为262.36mm，二者较850—2012年平均值分别偏高3.94mm(+1.5%)、偏低4.28mm(-1.6%)。MCA径流介于LIA与MWP之间，为265.33mm，较850—2012年平均值仅偏低1.31mm(-0.5%)。

由此来看，上游地区与中游地区不仅径流变率大小不同，而且3个特征时期径流的相对大小也不一致，最主要的差异为，上游暖期径流较冷期偏多，且径流的变率小。而中游则存在相反的特征，暖期径流少，冷期径流多，且径流变率大。这一特征与基于代用资料重建的干湿变化特征基本一致。从Chen等整理的多条基于湖泊沉积物重建的干湿变化序列来看，黄河上游MCA偏湿，LIA偏干；而黄河中游地区MCA偏干，LIA偏湿。对比MCA与LIA的旱涝频率也表明中国东部暖期易干旱，冷期则易多雨。

黄河上游属于典型的高寒气候，植被类型以草地为主，生态环境脆弱。高寒的气候特征使得上游冰川、积雪和冻土发育，其对径流变化有显著的调节作用，暖期的升温利于三者消融形成地表径流，反之亦然。这可能是导致黄河上游暖期的径流较冷期偏大主要原因。黄河中游则属于大陆性季风气候，且位于干旱半干旱气候区，对气候变化响应敏感。该地的蒸散发主要受水量供应条件限制，在水分充足的情况下，暖期升温将导致蒸散发增加，对径流有显著的消耗作用。因此，黄河中游径流与温度变化相反，暖期径流偏少，冷期径流偏多。

2.3过去千年各因子弹性系数的时空差异

黄河中、上游径流对各因子的弹性系数存在鲜明的差异。如图6所示，中游地区径流对降水与潜在蒸发的弹性系数绝对值高于上游地区，分别偏高约0.07和0.15(相对偏高约8.6%、37%)，而地表变化(surface)的弹性系数绝对值则低于上游地区，偏低约13.36(相对偏低约10%)。这一结果意味着中游径流深对降水与潜在蒸发变化的响应敏感度大于上游，对地表变化的响应敏感度则小于上游。其主要原因是上游与中游径流来源存在差异。中游地区径流主要来源于降水，致使其对降水响应极为敏感；同时中游地处半干旱—干旱气候区，大气蒸发力大，地表蒸发是当地水量平衡的重要分量，致使径流对潜在蒸发响应也较为敏感；而上游地区因其高寒的气候特征，地表径流同时来源于降水与冰雪融水，冰雪融水可以“缓冲”降水对径流的影响，且其高寒的气候特征致使大气蒸发力较弱，径流对潜在蒸散发响应敏感度较低。由此造成上游地表变化的弹性系数大于中游，而降水与潜在蒸发的弹性系数则小于中游。

同时，径流对各因子的弹性系数也受特征期之间气候冷暖转变的影响。上游地区，降水与潜在蒸发的弹性系数受气候冷暖转变影响微弱，二者在MCA、LIA、MWP3个时期大致相当，分别约为0.81、-0.42(图6)。而地表变化的弹性系数则随气候冷暖变化而异，其中LIA到MWP(冷转暖)的弹性系数(绝对值)略小于MCA到LIA(暖转冷)和MCA到MWP(持续偏暖)。中游地区，降水与潜在蒸发的弹性系数受气候冷暖转变影响较上游明显，二者在MCA到LIA(暖转冷)和MCA到MWP(持续增暖)大致相等，分别约为0.88、-0.57，而二者绝对值在LIA到MWP(冷转暖)略有增加，分别增至约0.89、0.58。LIA到MWP(冷转暖)地表变化的弹性系数(绝对值)略小于MCA到LIA(暖转冷)和MCA到MWP(持续偏暖)，此特征与上游地区一致。

由此来看，中游地区降水和潜在蒸发的弹性系数绝对值大于上游，且在冷转暖过程中弹性系数略大于暖转冷过程。这意味着气候增暖对中游地表径流的影响程度大于上游地区，这与王国庆等年际尺度的研究结果一致。同时，中游与上游有一个共同的特征，在持续偏暖过程中地表变化的弹性系数(绝对值)大于暖转冷与冷转暖过程，意味着气候持续增暖背景下，地表对径流的调控能力增强。

2.4特征时期之间径流差异及归因

黄河上游由MCA到LIA径流减少约2.98mm(表2)。其中，由降水主导的径流减少高达3.02mm，地表变化导致的径流减幅大约为0.34mm。同时，由MCA到LIA气候转冷，大气蒸发力下降，潜在蒸发减少，有利于径流增加，其导致的径流增幅约1.48mm，在一定程度上缓解了降水减少主导的径流减少。

由LIA到MWP，径流增加约8.22mm。其中，由降水增加主导的径流增加达10.74mm，是径流增加的主要成因。同时，由LIA到MWP气候转暖，大气蒸发力提升，潜在蒸发增加，由此导致径流减幅约为2.77mm；地表变化导致的径流减幅约为1.10mm。MCA到MWP，径流偏多约5.24mm，降水偏多仍是主要成因，由其导致的径流增幅高达7.72mm。MWP较MCA温度升高，蒸发力提升，潜在蒸发增加，由此导致径流减幅约为1.29mm；同时，地表变化导致的径流减幅约为1.44mm。

由此来看，上游地区由暖转冷，径流减少，由冷转暖和持续偏暖，径流增加。且降水与潜在蒸发的贡献是反向的，但潜在蒸发的贡献大约仅是降水贡献的1/6~1/2；地表变化的贡献持续为负，但其贡献更小，仅是降水贡献的1/5~1/10。因而，上游3个特征期之间径流变化普遍由降水变化主导，降水变化贡献的一部分被潜在蒸散发抵消。

黄河中游由MCA到LIA径流增加约30.3mm(表3)。其中降水增加仍是主要成因，其所致径流增幅高达27.64mm。同时，MCA到LIA气候转冷，大气蒸发力减弱，潜在蒸发减少，促使径流增加约7.52mm，而由地表变化导致的径流减幅为1.82mm。由LIA到MWP，径流减少约6.57mm，其中由降水主导的径流减少达4.75mm，同时气候转暖，潜在蒸发增强，致使径流减少约4.29mm，而由地表变化导致的径流增幅为0.09mm。

MCA到MWP，径流增加约23.73mm，由降水主导的径流增加高达22.89mm。同时，地表变化促使的径流减幅约1.73mm。尽管MWP相对MCA气候更为温暖，但风速等大气条件致使MWP相对MCA的潜在蒸发减少，进而促使径流增加，潜在增幅约3.23mm。

由此来看，中游地区由暖转冷和持续偏暖，径流增加，由冷转暖，径流减少。且降水与潜在蒸发的贡献是同向的，由暖转冷和持续偏暖，潜在蒸发的贡献大约仅是降水贡献的1/7~1/4，由冷转暖，潜在蒸发的贡献与降水的贡献大致相当。地表变化的贡献与降水贡献是反向的，但其数值不足降水贡献的1/10。因而，中游地区3个特征期之间径流变化是降水与潜在蒸发变化共同导致的。

2.5地表变化对径流影响的复杂性

本文采用“隐式”的方式计算了地表参数ω，在归因分析过程尚未明确地表参数指代的具体物理量。已有研究发现植被覆盖与地表参数ω显著相关，并给出了由归一化植被指数(NDVI)与气候季节性指数推算的ω半参数公式。近期，张建云等通过弹性系数法与复合函数链式求导法分析了黄淮海流域径流、NDVI和ω三者的关系，表明NDVI增加会使ω增加，进而导致径流量减少。

如表2~3所示，地表参数ω变化的贡献较小，在MCA到LIA、LIA到MWP之间，地表变化的贡献明显低于分析误差，仅MCA到MWP之间地表变化的贡献明显大于分析误差。为探究百年尺度上植被变化通过改变地表参数ω进而对径流的影响，本文分析了MCA和MWP时期的叶面积指数(LAI)之差。分析发现，中游地区MWP的LAI较MCA偏高9%，上游地区MWP的LAI较MCA偏高10%，这与先前研究指出的现代暖期CO2浓度高所致的施肥效应是一致的。叶面积指数(LAI)与NDVI均是反映植被状况的参数，二者呈显著正相关，较高的LAI会导致较多的冠层截留和气孔蒸腾，从而不利于形成径流。这与表2~3所示的MCA到MWP之间地表参数变化的负贡献是一致的。由此表明MCA到MWP之间地表参数变化的负贡献中的一部分来自于植被动态的调节。这可能意味着未来CO2浓度持续升高的情景下，CO2施肥效应造成的植被生长旺盛对水文过程的调节作用不容忽视。

3．结论与讨论

上述分析结果表明，过去千年黄河中、上游径流变化呈现出鲜明的百年尺度特征，其百年尺度变化与气候冷暖变化在中游与上游存在不同的组合关系，上游地区径流与冷暖变化位相一致，而中游地区二者变化位相则相反。随着百年尺度气候冷暖转换，上游地区径流变化由降水变化所主导，潜在蒸发则为负贡献，部分抵消了降水变化的贡献；中游地区径流变化也主要受控于降水变化，同时潜在蒸发也为正贡献，进一步加强了降水变化导致的径流变化。地表变化的贡献为负，但仅在2个暖期之间对比情形下地表变化的贡献大于分析误差。

已有研究重建了过去黄河流域不同地区的径流量变化，但对其变化的定量归因分析较少。本文弥补了百年尺度径流变化定量分析研究的不足，量化了各因子对过去千年百年尺度径流变化的贡献，明确了不同气候转变期各因子贡献的差异。研究结果可为过去千年不同气候变化背景下径流变化主要驱动因子的认识提供参考，并且丰富了径流变化归因的多尺度研究。更重要的是过去千年气候作为现代和未来气候变化研究的历史背景，厘清其气候冷暖转变与径流变化的关系，对预测未来不同气候变化情景下的径流变化具有重要意义。

但是，当前的研究仍然存在一些不足。首先，BCC-CSM1-1模式中没有汇流模块与水文流域尺度的径流深，本文使用了模式中流域内所有网格径流深的平均值，可能与实际径流深存在一定误差。其次，在计算潜在蒸散发与地表特征参数ω时，有些变量采用的月数据，平滑了一些要素的特征，但并不影响对径流百年尺度变化的研究。最后，由于BCC-CSM1-1陆地模块包含的地表参数较少，本文对ω的分析有限，有关ω对径流变化的影响有待进一步的研究。

（本文原刊《地理科学进展》2022年第7期第1226—1238页，文中原有注释，引用请务必参考原刊。）