伏牛山龙池墁南坡油松径向生长对气候变化的响应

李 轩,彭剑峰,2,3,*,李静茹,杨 柳,崔佳月,彭 猛,李成蹊

1 河南大学地理与环境学院,开封 475004

2 河南大学环境与规划国家级实验教学示范中心,开封 475004

3 河南省地球系统观测与模拟重点实验室,开封 475004

气候变化对陆地生态系统产生非常大的影响,尤其是对森林生态影响更加显著[1—2],不仅可以通过直接影响树木的代谢过程来影响生物量,还可以影响到森林的结构和功能[3—4], 进而影响着森林的分布格局[5—7],而且也会引发一系列的生态问题,导致树木生长下降、树木死亡事件增多及森林由“碳汇”向“碳源”转变的风险[8—9]。因此,理解气候变化对树木生长的影响机制对于准确预测树木个体的死亡时间及森林动态十分必要,也有助于更好的预测未来气候变化情景下的森林动态及森林演替的方向。树轮以其具有定年准确、分辨率高、连续性强、数据量测精确和样本分布广泛等特点,在全球气候变化研究中已被广泛应用[10],尤其是利用树轮资料(宽度、密度等)重建区域过去气候变化[11—13]。

经历多年的发展,树轮气候学研究在我国也已取得较大的进展,目前主要在高寒的青藏高原[13—16]和干旱的西北地区[17—20],大多研究集中在树木径向生长对气候因子的响应和利用长年表重建过去的气候及水文等变化；东部季风湿润区的研究起步相对较晚,如Cai等[21]在大别山分析了不同海拔的马尾松和黄山松对气候变化的响应；Duan等[22]利用湖南省、江西省5个采样点马尾松树轮资料重建了中国东南地区1849—2008年间1—4月平均气温；甘展峰等[23]研究了福建长汀福建松对气候的响应；陈峰等[24]用福建沙县马尾松树轮宽度分析与夏季-太平洋涛动指数的关系；张芳芳等[25]研究了湖北神农架地区巴山冷杉树轮δ18O序列的气候意义。目前虽然有一些关于全球气候变化对我国南北分界线伏牛山地区森林生态影响的研究[26],但在树轮研究方面仍多集中在该地区的优势种华山松,如Peng分析了白云山不同年龄华山松[27]、不同海拔华山松[28]及不同树种(华山松和日本落叶松)[29]对气候的响应,史江峰等[30]用华山松宽度序列重建了秦岭东缘冬季气温变化。大多数研究集中在树轮指数与气象因子的响应上,从微观尺度上将树木生理观测实验结合树轮的研究[31]以及宏观尺度上树木径向生长对森林管理措施比如间伐造成的小气候变化的响应、利用树轮资料估算区域森林碳汇的固碳能力的研究[32—33]较少。

伏牛山地属于秦岭余脉,是我国东部季风区暖温带的南界,属于生态交错带,具有明显的过渡性、脆弱性、敏感性、尺度性和动态性等特征,被认为是对气候变化反应最敏感的地区[34—35]；油松主要分布在我国的华北、西北和东北地区,伏牛山是油松分布南界[36],主要分布在海拔900—1800m的地区,因此这里也是油松树木年轮研究的理想区域,如Shi等[37]、Liu等[38],Peng等[11]都已在伏牛山及所属的石人山、白云山等地开展了油松树轮与气候关系的相关研究工作。龙池墁属于伏牛山山系,主峰海拔2129m,是油松分布的南界和上界,海拔较高且尚未开发,保存有大面积的原始森林,该地区仅见Shi等[37]、田沁花等[39]做了一些气候重建研究。在伏牛山大多数树轮的研究集中在海拔梯度和不同树种方面[28—29,37],然而树龄是森林结构、成分变化的重要诱导因素[31,40—43],它极大地影响着森林的生长动态乃至分布格局[41,44],尚未发现油松径向生长对气候变化响应年龄差异的研究。

本研究旨在分析龙池墁油松径向生长对伏牛山南、北坡气候的响应差异；探讨油松径向生长—气候响应的年龄差异并找出油松径向生长的显著气候因子；建立不同龄级油松气候—生长模型。这对评估气候变暖背景下龙池墁油松未来分布变化、利用油松树轮重建伏牛山地区气候有普遍的意义,也能为伏牛山地区森林管理和抚育更新提供科学依据。

1 数据获取和研究方法

1.1 研究区概况

龙池墁位于河南省洛阳市嵩县南部,东邻平顶山市鲁山县,西与西峡和栾川县相望(图1),主峰海拔2129m,属于秦岭东支余脉伏牛山部分,山地多由燕山期花岗岩组成,山坡大多悬崖峭壁,地表附有较薄的土壤,土质疏松,上部有较厚的枯枝落叶。气候属于大陆性季风气候区,多年平均气温13.76℃,年降水量844mm,雨热同期,气候属于大陆性季风气候区,位于北亚热带向暖温带过渡地带。植被在低海拔地区以栓皮栎(QuercusvariabilisBl.)、山茱萸(CornusofficinalisSieb. et Zucc.)等居多,高海拔地区以油松(PinustabuliformisCarr.)和华山松(PinusarmandiiFranch.)为主,夹杂高山杜鹃(Rhododendronlapponicum(L.) Wahlenb)等灌木[45]。龙池墁人口密度小,山势陡峭,很难进入,受人为因素干扰较少,因此保有自然原貌和原始森林,因此龙池墁是树轮研究的理想区域。

1.2 树轮数据获取

2020年11月,选择河南省伏牛山自然保护区龙池墁南坡(图1)原始森林区(33.71°N,111.98°E,海拔1600—1750m)的一片油松林区作为采样点,凡胸径约大于20cm的油松皆采,共采集35棵油松共53根样芯。

图1 采样点和附近气象站位置示意图

样本带回实验室后,采用国际树轮库的标准树轮研究方法,对样芯进行固定、打磨、交叉定年、宽度测量,并用COFECHA程序[46]对定年结果进行质量控制,剔除一些年轮异常、腐芯或测量困难的样芯,最终选取29棵树45根样芯。把通过COFECHA程序定年检验无误的45根样芯宽度序列,依据SPSS 26软件做系统聚类分析,明显地将45根样芯分为老龄和幼龄两个组(图2),其中老龄油松共12树20芯,年龄段介于94—183年间,平均年龄131年；幼龄油松共17树25芯,年龄段在44—74年间,平均年龄为60年。

图2 样芯年轮序列的聚类分析谱系图(横坐标为样芯的编号)

根据树轮聚类分析的结果,把整体、老龄和幼龄组的样芯分别用ARSTAN程序[47]形成年表。首先要进行标准化,即去树木生长趋势,线性函数适用于年龄较短、没有明显生长趋势的树轮样本,负指数函数能较好地拟合树轮宽度快速下降和之后缓慢递减的趋势。标准化之后再进行均值化,常见的方法有算数平均和双权重平均法[48—49],最后获得一系列树轮年表,其中3个标准年表(STD)和样本量见图3,各组年表的统计特征见表1。

图3 不同龄级树轮宽度标准年表(STD)及样本量(竖线为SSS>0.85的起始年份)

表1 整体、老龄和幼龄树轮宽度标准年表(STD)统计特征值

1.3 气象数据

选取离采样点较近的西峡站(33.3°N,111.5°E,252m)和栾川站(33.78°N,111.6°N,751.5m) 2个气象站(数据来自中国气象数据共享服务系统(https://data.cma.cn/))1961—2016年的月平均气温、月平均最高温、月平均最低温和月降水量4个气象因子指标作为分析对象。西峡站的多年月平均温15.27℃,年降水量847.92mm,栾川站的多年月均温12.26℃,年降水量841.45mm(图4)所示,它们降水都主要集中在夏季,具有雨热同期的特点。PDSI月份数据来自Climate Explore(http://climexp.knmi.nl/)离采样点最近的格点(33.75°N,111.75°N)数据。

图4 西峡站和栾川站气象站资料(1961—2016)

1.4 研究方法

本研究利用树木年轮学的专业软件Dendroclim2002[50]对树木年轮标准年表(STD)与两个气象站4个气象指标及PDSI指数(1961—2016年)做皮尔逊相关分析,从上年5月到当年11月共19个月来分析伏牛山南、北坡气候对油松径向生长的影响,选择相关较好的气象站；利用相关较好的气象站和不同年龄油松的标准年表做相关分析,对比不同年龄油松径向生长—气候响应差异并找出影响显著的气候因子；用SPSS 26软件对相关显著的气候因子和树轮标准年表做回归分析,建立不同龄级油松的生长模型。

2 结果分析

2.1 树轮宽度年表的统计特征

龙池墁油松树轮宽度年表的平均敏感度(MS)、信噪比(SNR)、总体样本解释量(EPS)较高,说明油松的标准年表保留了较多的气候信息(表1)。为了对比分析不同年龄组的公共统计特征,研制年表时采用统一的公共区间1980—2010年。研究发现老龄油松具有更高的平均敏感度和标准差,且整体组的信噪比SNR(27.959)和样本总体解释量EPS(0.965)最高,幼龄组油松次之,这与Peng 等[27]在白云山的研究相似。由于采样点坡度较大和采样难度较大,获取的样本量相对较少,一些树只取到单芯,油松的样芯间相关系数(R1)、树内的相关系数(R2)和树间的相关系数(R3)相差不大,表明采样点油松的树木年轮宽度变化具有较高的一致性。幼龄油松的相关值稍大,这可能与幼龄油松的样本量较大有关。

2.2 油松径向生长与气候要素的关系

2.2.1整体油松标准年表与气候要素的相关分析

用Dendroclim 2002软件分析了整体组油松标准年表与栾川、西峡两站气象因子的相关性,其结果如图5。

图5 整体油松树轮宽度标准年表与月均温、月均最高温、月均最低温和降水量的相关分析

整体油松与栾川站的当年5月平均温度达到显著相关；与当年5月、7月、以及上年7月、8月的平均最高温显著负相关；与当年5月降水量显著相关,相关系数仅0.29。与西峡站的当年5月平均温度相关显著；与当年5月、7月的月平均最高温显著相关；与当年的5月、6月的月平均最低温相关显著；与当年3、5月降水量显著相关。

与各月份的PDSI指数呈正相关(图6),其中当年4—9月达到显著相关,当年5月份的相关系数最高,达到0.52；与当年5—7月PDSI指数相关系数为0.54。

图6 整体组、老龄组和幼龄组油松树轮宽度标准年表与月PDSI指数(1961—2016)的相关系数

整体油松标准年表与两站春末及夏季的水热因素相关显著,尤其是与当年5月气候因子的相关系数最大,这与彭剑峰等[51]在神农山地区的研究结果基本一致。从单个气象要素来看整体油松与同是位于龙池墁南侧的西峡站的相关系数较高,例如与西峡站当年5气温的相关系数要比栾川站高0.06—0.22之间,月降水量的相关系数要比栾川高出0.11。

2.2.2不同年龄组油松树轮宽度标准年表与气候要素的相关分析

鉴于整体油松树轮宽度标准年表与西峡站气象因子的相关较高,不同年龄油松径向生长仅与西峡站气候因子做相关分析(图7)。

图7 老龄、幼龄油松树轮宽度标准年表与西峡站气候因子的相关分析

从相关分析的结果来看,不同年龄组油松与与西峡站春末夏初的气温降水相关显著。老龄油松仅与当年5月的平均温显著相关；幼龄油松与上年7月和当年5、6、7月平均温达到显著相关。老龄油松与当年4、5月平均最高温显著相关,幼龄油松与上年7月和当年的5、7月平均最高温显著相关。老龄和幼龄油松均与当年的5、6月平均最低温显著相关。老龄油松与当年5月降水量相关显著,幼龄油松与当年3、5月降水量相关显著。

老龄油松与当年4—8月单月的PDSI指数相关显著(图6),其中与当年5、5—7月的相关系数最高都为0.47以上；幼龄油松与当年4—9月单月的PDSI指数显著相关,其中与当年5月份的相关系数最高为0.52,与当年5—7月的PDSI指数相关系数达到0.55。

总之,老龄、幼龄油松受5月温度和降水影响显著,而幼树限制因子较多,不仅受当年5月的温度和降水影响,还会受到6、7月高、低温和3月降水的影响,PDSI指数的限制作用基本一致,5—7月相关系数最高。

3 讨论

3.1 影响龙池墁油松树轮宽度生长的气候因子

春夏初的水热条件是龙池墁油松的径向生长的限制因素,尤其是5月份,与其他相关显著的月份相比,当年5月的相关系数要高出0.04—0.17。春夏之交的高温对油松径向生长有显著的限制作用,水分充足的条件则有利于油松的生长。与Peng等[27—28,51]在河南山白云、伏牛山、神农山；田沁花等[39]在伏牛山的研究及Zhao等[52]在白云山、龙池墁的研究结果基本一致,然而刘禹等[38]在石人山(海拔2010m)的研究表明油松生长与当年1—3月、5月平均最低温正相关,与本研究结果有一定的差异,可能是采样点位置的不同。龙池墁位于暖温带大陆性季风气候区,5月份雨季尚未来临,而此时气温虽然能满足油松光合作用的热量需求,但水分条件欠缺导致油松出现明显的干旱胁迫会使细胞气孔关闭,出现“碳饥饿”现象,导致光合速率下降[53—54]。5月份高温少雨不会使初级生产力(GPP)变大,而夜间温度高,呼吸作用强又加剧能量消耗,导致净初级生产力(NPP)下降[55]。因此5月的水热条件不仅是龙池墁油松径向生长的限制因子,也是整个豫西山地油松径向生长的限制因子。

本研究中油松树轮宽度标准年表(STD)与气候因子、PDSI相关分析表明,单一的温度、降水量因子对油松的影响不如综合性指标PDSI指数显著,如与当年5月份的气温和降水量相关系数最高分别是-0.42和0.4,而与当年5月的PDSI指数则高达0.52。PDSI指数不仅综合考虑了温度和降水,而且还顾及到前期干旱事件对土壤水分的影响[56],因此油松径向生长对PDSI的响应最好。与PDSI指数月份组合的相关系数大于单月,其中与当年5—7月PDSI指数为最高为0.55,可以看出明显的气候因子的累积效应明显,这与Zhao等[52]在伏牛山、Peng等[27]在白云山的研究结果基本一致。

3.2 油松径向生长对伏牛山两侧气候响应的差异

龙池墁南坡油松径向生长受伏牛山南侧气候影响较大。伏牛山大致呈东西走向,对冬夏季风都有明显的阻挡作用,伏牛山两侧的气候差异显著。栾川与西峡的气候特征相似,主要区别在1月气温。栾川气象站位于伏牛山北侧深受寒冷的冬季风影响,1月份月均温在0℃以下,属于暖温带大陆性季风气候；西峡气象站位于伏牛山南侧,冬季冷空气受到地伏牛山阻挡,1月份均温在0℃以上,属于亚热带季风气候。采样点和西峡站都位于龙池墁的南侧因而相关好,与龙池墁北侧的栾川站的相关稍弱,刘禹等[57]在秦岭中部利用南、北坡油松树轮温度重建研究也有类似的结果。在今后的研究中要重视坡向因素对油松径向生长的影响。

3.3 不同年龄组油松对气候因子的响应差异

从不同年龄油松的标准年表与西峡站气候因子的相关系数来看,上年7月平均温、平均最高温对幼龄油松呈现出明显的“滞后效应”。刘禹等[19]在贺兰山油松的研究中也发现气候因子的滞后效应,张芬等[58]在祁连山东部油松径向生长对气候的响应中表明不同年龄油松对气候因子的响应存在差异,幼龄油松也存在“滞后效应”。幼龄树的根系不发达,幼龄树不能从更广阔、更深的区域获取生长必需的水分,因此幼龄油松对外界环境因子的抵抗力弱[33,59—60]。

3.4 不同年龄油松1990年以来生长下降的差异

无论老龄油松还是幼龄油松的树轮宽度在近十几年明显出现降低趋势,这里固然有油松年龄的缘故,但气候因子造成的影响也不容忽视[61]。与老龄油松相比,幼龄油松在20世纪90年代以后的树轮的平均宽度下降迅速、幅度很大(图8),其中1995年和2014年呈现极端窄轮,尤其是2014年年轮宽度值最小。与相应的前5年平均宽度相比,老龄油松在1995、2014年树轮宽度分别下降了16%、50.2%；而幼龄油松在1995、2014年树轮宽度分别下降了26.7%、54.9%,显然不同年龄油松对相同的气候条件树轮宽度下降的幅度存在差异,幼龄油松下降速度远大于老龄油松(图8)。Lambert等[62]在加拿大的树木结合胸径和树高地上生物量方程研究中提到树桩(stump)的生物量可以代表75%的总生物量,因此龙池墁南坡油松树轮宽度出现不同程度的下降现象是否是生长衰落现象呢？值得进一步研究。

图8 5—7月PDSI及西峡站5月平均温和老龄树、幼龄树树轮的平均宽度变化对比

在全球温度变化趋势短期内不变的大背景下[63],一些报道显示北方针叶林的生长发生改变从而导致分布变化的风险上升[8—9]。西峡站气象资料显示,20世纪90年代以后5月平均温度呈上升趋势(图8),夏季初期高温不利于龙池墁油松的生长,尤其是低海拔地区,这与Peng等[11]、Shi等[37]在伏牛山地区的研究相似,Li等[64]用 Maxent模型对我国杉木(Cunninghamialanceolata(Lamb.) Hook.)分布模拟结果表明其分布的重心有北移的趋势,然而刘禹等[38]在石人山的研究中提到海拔2010m处的油松树轮捕捉到了20世纪气候升温的信号且气温上升有利于油松的生长。本研究的采样点海拔较低(1650—1750m),幼龄油松径向生长呈现出与年龄不相符的生长下降,而老龄油松由于年龄因素必然不断衰退,油松的生长量下降甚至死亡数量增加,导致油松的更新能力降低甚至失去更新能力。龙池墁较低海拔地区未来的森林生态因种群改变而降低地区生态稳定性、增加生态失衡的风险。

4 龙池墁油松生长模型

本研究尝试分析(1961—2016年)整体组、老龄及幼龄油松的标准年表指数(Wzt、Wot、Wyt)与Tc5、Tp7、Pc3、Pc5、PDSIC5-c7等显著因子利用多元线性回归分析方法和步进回归方法,分别建立不同年龄油松的生长模型,以对比两种方法结论的异同。

4.1 油松径向生长多元线性回归模型的建立

基于油松树轮标准年表与西峡站气象因子的相关分析,发现影响龙池墁油松径向生长的主要限制因子为当年5月平均温(Tc5)、上年7月平均温(Tp7)、当年3、5月降水量(Pc3、Pc5)。

为了确定所有的显著气象因子与树轮指数间的定量关系,采用多元线性回归的方法,把所有的显著因子Tc5、Tp7、Pc3、Pc5加入生长模型以探讨它们对龙池墁油松径向生长的影响,整体组、老龄及幼龄油松的多元线性回归生长模型为(1)、(2)、(3)式,它们的统计量(表2)。

表2 不同年龄组标准年表与Tc5、Tp7、Pc3、Pc5的多元线性回归模型的统计量

Wzt=30105+0.002×Pc3+0.001×Pc5- 0.049×Tp7-0.046×Tc5

(1)

Wot=2. 541+0.001(Pc3+Pc5)- 0.031×Tp7-0.042×Tc5

(2)

Wyt=3.516+0.002×Pc3+0.001×Pc5- 0.063×Tp7-0.047×Tc5

(3)

4.2 油松径向生长线性步进回归生长模型的建立：

本研究再用步进回归的方法,设定自变量F≤ 0.05进入回归模型,若F的概率≥ 0.1则剔除变量,确定主导气候因子与树轮指数的定量关系。将所有Tc5、Tp7、Pc3、Pc5做步进回归逐步剔除最终保留最重要限制因子,整体组、老龄及幼龄油松的步进回归生长模型分别为(4)、(5)、(6),确定不同油松径向生长与最显著气象因子间的定量关系。

Wzt=2.313-0.064×Tc5

(4)

Wot=2.143-0.057×Tc5

(5)

Wyt=2.409-0.068×Tc5

(6)

不同年龄组标准年表与Tc5、Tp7、Pc3、Pc5步进回归分析的最终都只保留了当年5月平均温,可见5月均温是油松的主要限制因子,其相关统计量见表3。

表3 不同年龄组标准年表与Tc5、Tp7、Pc3、Pc5的步进回归模型中Tc5的相关统计量

鉴于PDSI指数是复合指标,因此将3个标准年表年表单独与当年5—7月PDSI指数建立线性回归生长模型。不同年龄组标准年表树轮指数与PDSIC5-c7指数做线性回归模型的相关统计量见表4。

表4 不同年龄组标准年表与PDSIc5-c7指数的线性回归模型的相关统计量

Wzt=0.968+0.067×PDSIc5-c7

(7)

Wot=0.942+0.057×PDSIc5-c7

(8)

Wyt=0.988+0.072×PDSIc5-c7

(9)

通过以上的生长模型可以看出当年5月平均温和5—7月PDSI指数对各年龄组油松生长模型都能保留较高的解释量,尤其是5—7月的PDSI指数对龙池墁油松径向生长的影响显著,这为今后气候重建打下良好的基础。

不同年龄组油松对气候因子响应分析结果也与生长模型结果一致,都与当年5月平均温度及当年5—7月PDSI指数具有较高的相关性。

5 结论

树轮STD年表的信噪比(SNR)、样本总解释量(EPS)以及敏感度(MS)值都较高,表明龙池墁南坡油松树轮包含较丰富的环境信息,适合进行树轮与气候因子关系研究。相关分析发现龙池墁较低海拔地区油松径向生长的主要限制因子是当年5月水热因子和5—7月PDSI指数,并且与龙池墁南坡气候因子的关系比北坡密切。

不同龄级的油松对气象因子的响应存在差异,幼龄油松的径向生长受环境因素的影响较老龄油松强烈；幼龄油松与气候因子的相关系数较大且达到显著的月份多；幼龄油松的径向生长不仅受当年夏季高温的制约还受上年夏季高温的影响,呈现明显的滞后效应,幼龄油松出现了与年龄不相符的生长下降现象。

通过对与龙池墁油松径向生长的显著因子和树轮标准年表指数做多元和步进回归分析,建立3个年表的生长模型,确定了树轮指数与西峡站的5月平均温和5—7月PDSI指数间的定量关系。在目前全球气候变化趋势下,龙池墁山地低海拔地区将更加不利于幼龄油松的生长,随着老龄油松不断死去,最终会危及油松种群在低海拔地区生态群落中的地位,甚至造成龙池墁山地油松带将向更高海拔迁移,增加该地区生态失衡的风险。