基于红皮云杉（Picea koraiensis）重建大兴安岭南麓历史径流量

2023-09-28魏英楠孙柏林

干旱区地理(汉文版) 2023年8期

魏英楠，马龙，孙柏林，张晶

（内蒙古农业大学水利与土木建筑工程学院，内蒙古呼和浩特 010018）

气候变暖导致冰川持续退缩，高纬度地区、高海拔山区多年冻土层融化［1-2］，水资源供需矛盾加剧，对全球和区域尺度上降水量、径流量等水文要素产生深刻影响。大兴安岭地处高纬寒旱地区，是我国半湿润区与半干旱区、季风区与非季风区的分界线［3］。独特的气候条件和地理位置导致其成为全球气候变化敏感地带之一，揭示该地区长历史时期水文变化特征的研究亟待开展。然而器测数据往往少于百年，极大地阻碍了对历史水文变化的深入了解。

树木年轮定年准确、易于获取及分辨率高等优势［4］使其成为延长水文资料的重要手段之一。自20世纪30 年代美国内华达州首次利用树轮延长历史年径流量资料后［5］，树木年轮在全球范围被广泛应用于重建长时间尺度水文资料。北美［6-8］和欧洲地区［9-11］研究成果较多，国内主要集中在新疆［12-14］、青藏高原［15］及大兴安岭［3,16］等水文气候要素限制作用显著的高纬度寒旱区域。大兴安岭地区树轮年代学的研究起步较晚，多以兴安落叶松（Larix gmelinii）和樟子松（Pinus sylvestris）为代用材料，例如分析兴安落叶松和樟子松径向生长与气候变化的响应及差异，发现气候变化改变了树木生长与气候间的响应关系，樟子松受气候变暖的影响更加显著［17-18］；基于樟子松宽度年表重建伊敏河1868—2002 年年径流总量［19］，重现了伊敏河过去135 a 经历的6 个干旱时期和4 个湿润时期；利用樟子松首次建立海拉尔河整个流域及11个支流过去202～216 a年径流量序列［20］，发现海拉尔径流量的变化与厄尔尼诺-南方涛动（El Niño-Southern Oscillation，ENSO）密切相关，并预测随着ENSO 活动的增强，海拉尔河流域将面临极端水文事件。

白音敖包自然保护区地处大兴安岭山地南麓，位于东北、华北和蒙古植物区系交汇地带［21］，原始林木资源丰富，且处于半干旱地区，发源于大兴安岭最高峰的贡格尔河是保护区内主要河流，水分条件对该区域树木径向生长限制作用显著［12］，在此开展树木年代学研究具有得天独厚的优势。红皮云杉（Picea koraiensis）是松科云杉属常绿乔木，耐荫耐寒耐干旱，在不同立地条件下均能生长［22］，被列入《世界自然保护联盟》（IUCN）濒危物种［23］，是保护区内树龄较长、树轮边界明显且对水文气候要素响应敏感的珍贵树种，具有较高的树木年轮学研究价值。前人利用红皮云杉重建了白音敖包过去160余年的降水量序列［24］，但关于径流量历史变化的研究较少。基于此，本文利用红皮云杉树轮宽度年表，重建大兴安岭南麓近2 个世纪的径流量历史序列，研究结果重现了大兴安岭南麓贡格尔河流域无数据记载时期的径流量变化，为其今后水文相关研究提供参考。

1 资料与方法

1.1 研究区概况

采样点白音敖包国家级自然保护区（图1）地理位置介于大兴安岭南端和锡林郭勒高平原之间［21］，地形为大兴安岭山地向蒙古高原的过渡地带。气候类型属于大陆性寒温带半干旱森林草原气候，寒暖分明，水热同期［24-25］。贡格尔河横卧于保护区北侧，境内流经长度14 km，敖包河是其主要支流。

图1 研究区采样点和水文气象站分布Fig.1 Distributions of sampling point and hydrometeorological stations in the study area

1.2 水文与气候数据

本文使用的气温、降水量资料为林西气象站1951—2016 年逐月平均、平均最高、平均最低气温和降水量（图2a），数据来源于中国国家气象信息中心。径流量资料为距离研究区最近的锡林浩特水文站1972—2016年逐月径流量（图2b），数据来源于当地水文勘测局水文年鉴。

图2 研究区多年月平均气温、降水量与径流量变化Fig.2 Changes in monthly average temperature,precipitation and runoff for many years in the study area

研究区具有明显的大陆性气候特征，冬夏气温较差大，降水量年内分配极不均匀。多年平均气温4.73 ℃，7 月气温值最高，1 月气温为全年最低值。多年平均降水量约为373.17 mm，呈单峰分布，7 月降水量为全年最高，5—9 月降水量占全年的90%。多年平均径流量约为0.17×108m3，4 月径流量最大，12月至次年2月为枯水期。

1.3 树木年轮采样及年表建立

年轮样本采集时间为2017 年8 月，采样时严格遵循国际树木年轮数据库（ITRDB）标准［26］，共采集37 株树的80 根样芯。待树轮样本在阴凉通风处自然晾干后，使用由粗至细目数不同的砂纸逐级打磨，直至年轮边界清晰可见。利用WinDENDROTM年轮分析系统［27］测定每一生长轮的日历年。使用COFECHA 程序［28］检验交叉定年和宽度测量结果，对缺轮、假轮、定年错误等问题进行检查，最终78根树轮样本用于年表的构建。

树木的径向生长不仅受自身遗传因子的控制，还与其生长过程中所处环境密切相关。为了解树木承载的水文气候变化信息，需去除样本自身生长趋势［4］。使用ARSTAN 程序［29］负指数函数［30］去除树木本身遗传因子及竞争干扰产生的生长趋势，最后利用双权重平均法合成标准化年表（图3a）。

图3 红皮云杉标准化年表特征Fig.3 Standardized chronological characteristics of Picea koraiensis

基于树木年轮学研究的复本原理，利用子样本信号强度（Subsample signal strength，SSS）确定年表的最低复本量（图3b）［31］，为保证重建序列可靠性并保留年表最大有效长度，限定SSS＞0.85 时为年表的可靠时段。

通过分析年表的统计特征值［32］，可以判断其是否适用于树木年轮学研究。年表的主要统计特征值见表1。平均敏感度（MS）是判断一个年表质量优劣的重要指标之一，其数值越大，树木径向生长对周围环境变化越敏感。一般认为MS＞0.2 的序列适合进行树木年轮的研究［33］；标准差的数值越大，树木群体受到气候因素的控制越剧烈；一阶自相关系数反映气候条件对树木生长的持续性影响，其数值越小年表质量就越高；信噪比是用来衡量样本中包含环境信息量大小的统计量，其数值较高时，年表中含有较多气候信息，适用于进行树轮气候学研究；样本总体代表性数值越大，样本的代表性就越好，一般样本总体代表性＞0.85 时所建立的年表比较可靠。红皮云杉年表统计特征值均满足树木年轮学标准，说明年表质量较好，包含较多气候信息，可用于树木年轮气候学的研究之中。

表1 红皮云杉年表主要统计特征Tab.1 Main statistical indicators of Picea koraiensis standardized chronology

1.4 研究方法

（1）使用Pearson 法对年表与气温、降水量、径流量等要素进行相关分析［34］。

（2）采用逐一剔除法检验重建模型的稳定性和可靠性。

（3）为研究径流量变化趋势，利用Z-score 法对重建径流量历史序列进行标准化处理，定义标准化后的径流量连续2 a 以上位于0m3之上的时段为变丰期，连续2 a 以上位于0m3之下的时段为变枯期。计算公式如下：

式中：R为径流量重建值（108m3）；mean 为重建径流量序列的平均值（108m3）；SD 为重建径流量序列标准差（108m3）。

（4）使用变异系数（Cv）反映径流量历史重建序列变化的剧烈程度，计算公式如下：

（5）定义重建径流量高于mean+σ的年份定义为丰水年，低于mean-σ的年份定义为枯水年，连续经历2 a及以上丰水年的时期为丰水期，连续经历2 a以上枯水年的时期为枯水期。其中：

式中：σ为重建径流量序列的标准差（108m3）；Xi为第i年对应的径流量（108m3）；μ为重建径流量序列的平均值（108m3）；N为样本数量（个）。

（6）径流量重建序列的周期变化分析，使用Morlet 小波分析法［35］进行分析，在Matlab 中，其计算公式为：

式中：ψ(t)为小波函数；t为时间（年）。

Morlet小波伸缩尺度（a）与周期（T）有如下关系：

式中：ω0为无量纲频率，通常取6.2附近的经验值。

（7）采用Mann-Kendall 检验法对径流量突变年份进行确定。

2 结果与分析

2.1 树木径向生长与水文气候要素间响应关系的分析

考虑到树木的生长在10月底之前已经结束，冬季水文气象因子对树木年轮生长具有滞后效应［36］，因此利用皮尔逊相关法分析树轮宽度年表与研究区前一年11—12 月和当年1—10 月水文气候资料的响应关系（图4）。

图4 年表与水文气候要素的响应关系Fig.4 Response relationship between chronology and hydroclimatic elements

除降水量外，树轮宽度年表与当年的径流量、气温关系更加密切。年表与各月平均气温、平均最高气温和平均最低气温均呈现显著负相关关系，即过高的气温不利于红皮云杉的生长。年表与6月气温负相关关系极显著，夏季研究区气温处于全年最高值，树木正值生长季旺期，此时气温升高，加强了树木的呼吸、蒸腾作用，导致土壤中本就匮乏的水分不能满足树木的生长所需，从而起到负面效果［37］；年表与降水量表现为正向响应关系，但除前一年11月以外，其余月份的降水量均未能通过0.05 的显著性检验，可能是降水年内分配极不均匀，对树木径向生长限制作用微弱。

红皮云杉宽度年表与径流量呈显著正相关关系，3—10 月径流量通过了0.05 显著性检验，所以年表对径流量的响应更加敏感。径流并非直接作用于树木生长，而是通过气温、降水等因子与其连接起来［38］。研究区属于半干旱地区，水分对红皮云杉生长过程起促进作用，但12 月研究区气温较低，树木处于休眠状态，此时水分过多将导致土壤冻结，树根发生冻害，不利于树木生长［39］；春季气温回升，气候适宜树木生长，树木开始萌发幼芽，但根系尚处于休眠状态，吸收功能较弱，径流可促进根系吸收水分；夏季树木进入生长旺盛时期，形成层细胞活跃，分裂活动旺盛，充足的水分条件有利于植物进行光合作用，形成宽轮；秋天气温下降，树木蒸腾作用减弱，准备进入休眠状态，此时充足的水分有助于植物为次年的生长储存能量［40］。

2.2 径流量重建

基于年表与水文气候要素的响应分析结果，结合红皮云杉的生长特性，本次选择以3—10 月组合月径流量为因变量，树轮宽度年表为自变量，利用幂指数回归方程重建了研究区1845—2016 年径流量序列，重建方程如下：

式中：yi为第i年组合月径流量（108m3）；xi为树轮宽度标准化年表序列；r为组合月径流量与年表的相关系数；N为月径流量实测资料长度（年）；R2adj为调整后r2；F为重建方程方差检验；P为重建方程显著性值。

利用国际树木年轮研究中的逐一剔除法对此重建模型进行检验［30］。逐一剔除法主要采用的统计参数有相关系数（r）、误差缩减值（RE）、有效系数（CE）、符号检验、一阶差符号检验、乘积平均值检验（t）等。其中RE大于0.3代表拟合程度较好，越接近于1，方程可靠性越高［41］；符号检验利用每对样本的数据差值符号进行检验，可检测出样本间是否存在显著性差异；t不但包含了对符号的检验并增加了序列数值的检验，能够直接反映出重建值是否包含有水文气候信息，其值越大表明2 个序列越相近。重建序列和实测序列的r=0.815，一阶差r=0.769，达到了95%的置信水平。以1972—1994 年为校准期、1994—2016 年为检验期对重建方程进行检验，校准期r=0.677，检验期r=0.881，均通过了0.01 显著性检验；重建方程的RE=0.411，CE=0.640，二者远远大于0，说明重建序列可靠；重建方程的原序列符号检验（32+/13-）通过了0.05 显著性检验，但一阶差符号检验（23+/21-）未能通过，说明方程对低频变化的重建更加准确［42］；乘积平均值（t=3.756）通过了0.05 显著性检验，表明径流量重建序列包含较多水文信息。

通过对比锡林浩特水文站观测时段（1972—2016年）径流量实测值与重建值（图5），二者的拟合效果较好，在低频的变化上能较好对应，变化的趋势总体一致，在个别时段，实测值与重建值基本吻合，进一步证明重建方程的可靠性。

图5 径流量重建值与实测值的对比Fig.5 Comparison of the reconstructed values and measured values of runoff

综上所述，重建方程的各项检验参数均通过了检验标准，可用于重建大兴安岭南麓贡格尔河流域历史径流量。

2.3 重建径流量历史序列的变化特征

贡格尔河1845 年以来径流量历史序列存在明显的枯水和丰水变化（图6a～b），径流量序列整体呈下降趋势，变化频率较高，M-K 检验显示在1849、1854年和1993年发生过径流量突变。近172 a里贡格尔河流域经历了7 个变枯期和6 个变丰期，出现了1853—1855 年、2000—2010 年2 个枯水期和1869—1873 年、1911—1918 年、1952—1958 年3 个丰水期，存在30 个丰水年及18 个枯水年，占比分别为17.4%和10.5%。

图6 径流量历史重建序列变化特征Fig.6 Change characteristics of historical reconstruction runoff series

可以看出，丰水年、枯水年出现年份较集中，1845—1852 年径流量表现为下降趋势，变化波动较小（Cv=19.80%），年代际均值较低；1853 年出现首个连续了3 a 枯水年，1854 年径流量由枯水期快速向平水期转变，1856—1866 年丰、枯水年零星出现；1867—1873 年共经历了6 个丰水年，1874 年径流量由丰水期向平水期转变，变化波动较剧烈（Cv=34.42%）；20 世纪10 年代丰水年大量出现，研究区经历持续时间最长的丰水期，年代际均值达到最高，序列波动较小（Cv=23.12%）；1919 年径流量表现出变枯趋势，直至1929 年缓慢回升，1934 年研究区迎来丰水年，年代际均值较上一时期有所提升，20世纪30 年代后期径流量向平水期转变；20 世纪50年代研究区经历最后一个连续7 a的丰水期后，径流量表现出显著的变枯趋势，变化波动（Cv=16.40%）较上一阶段（Cv=26.89%）轻缓，年代际均值明显下降；在经历了约30 a 的平水期后，径流量在1993 年发生突变，变化剧烈（Cv=39.10%）且下降速率较快，最终20世纪末期枯水年集中出现，导致贡格尔河流域历经了持续时间长达11 a 的枯水期，年代际均值达到近2 个世纪的最低值；21 世纪末期径流量缓慢回升，2016年径流量达到丰水年标准。

小波分析发现贡格尔河径流量重建序列在95%的置信水平上存在着显著的3 a、7～12 a、15～22 a及30 a周期（图7a），小波方差图（图7b）显示15～22 a周期贯穿整个历史时段，信号最强，表现最显著。15～22 a 周期可能与太平洋年代际涛动（Pacific decadal oscillation，PDO）有关，亚洲地区气候变化与PDO 关系密切［43］，研究区气候变化受其影响较大；7～12 a 周期与ENSO 变化周期基本一致，表明ENSO活动对研究区存在影响。径流量可以表征大区域的水文气候效应，是检测太阳活动对气候影响的指标之一［44］，研究区11 a 周期与太阳黑子活动周期对应；准3 a 周期是我国半干旱地区重要周期［45］，研究区水文变化与其密切相关。

图7 径流量重建序列周期变化特征Fig.7 Periodic variation characteristics of reconstruction runoff series

2.4 与历史事件及其他重建结果的比较

贡格尔河发源于3000 a 前，是研究区第一大内流河，其丰富的水源为沿岸居民与动植物提供良好的生存条件，是北方极为重要的生态屏障。为验证径流量重建序列的准确性，将重建结果与当地历史典籍［24,46］等资料进行对比，发现重建丰枯时期与历史气象灾害存在对应。

19 世纪50 年代至70 年代，研究区降水量较少导致贡格尔河径流量呈变枯趋势，连续经历3 a枯水年；民国元年至九年，赤峰多地接连暴雨，河水暴涨冲毁堤坝，导致山洪暴发，贡格尔河流域迎来持续时间最长的丰水期。20 世纪50 年代夏、秋雨涝，数千公顷农田被毁，径流量重建序列最后一个丰水期在此时出现。20 世纪60 年代内蒙古降水量时空分布极不均匀，1965 年研究区遭受建国以来最严重的旱灾，降水量比常年少50%～90%，连年干旱和上游牧民筑坝引水导致贡格尔河的水量大幅度减少，径流量序列呈变枯趋势，1968 年研究区处于枯水年。1970 年修建五道石门水电站，下游牧场灌溉供水减少，径流量显著回升，水电站发电记录显示在1992年发电量达到历史最高峰，1992—1993 年贡格尔河径流量连续2 a处于丰水年标准之上。然而20 世纪末期全球气候变暖，大兴安岭冰川融化和连年干旱导致贡格尔河水量骤减，同时河流上游黄岗梁地区开发铁矿，清洗矿砂需要大量河水，工业废水对河水造成污染，导致生态严重退化，河水水位持续下降，贡格尔河历经长达10 a的枯水期；当地政府为解决水面缩减、水位下降等问题，实施了保护治理、生态补水等一系列修复措施，研究区环境得到改善，21 世纪初径流量开始缓慢回升，在径流量重建序列终点已达到丰水年标准。

为增强此次重建结果的可信度，将大兴安岭西部利用樟子松重建的伊敏河流域1868—2002 年年径流量［19］和利用落叶松重建的海拉尔河流域近2个世纪径流量［20］与本次重建结果的重叠时段进行比较，贡格尔河与伊敏河（r=0.278，N=136，P＜0.01）及海拉尔河径流量重建序列（r=0.216，N=162，P＜0.01）均呈显著正相关关系，通过了0.01 显著性水平检验，说明三者关系密切，具有相似的变化趋势，验证了本次重建结果的可靠性。