森林生态系统对全球气候变化响应强烈, 特别是高海拔及高纬度地区生态系统类型的组成、结构、分布和生产力会随气候变化而发生显著变化^[1]。植被净初级生产力(NPP)作为植被对气候变化响应的主要衡量指标^[2], 可反映研究区的气候变化, 树木年轮的宽度亦可真实的记录气候变化^[3], 通过二者结合来研究森林生态系统对气候变化的响应具有重要的意义。

NPP相关研究模型主要包括统计模型、光能利用率模型、过程模型等类型^[4-5]。统计模型以Maimi模型、Thornthwaite Mermorial模型、Chikugo模型为代表, 利用气温与降水等气候因子来估算NPP, 但估算的误差较大^[6-7]。光能利用率模型主要包括CASA、GLO-PEM等模型^[8-9], 是以卫星遥感数据为基础的模型, 植被参数易于获取, 但是遥感数据可能受到天气等的影响, 且遥感数据时间序列较短。过程模型基于生理生态机制模拟NPP, 以Biome-BGC模型、CENTURY模型、TEM模型为主要代表模型, 估算结果比较准确^[4]。其中, Biome-Bio Geochemical Cycles(Biome-BGC)模型是从FOREST-BGC模型发展来的生物地球化学循环模型, 主要用于模拟生态系统能量、碳、氮、及水分的循环, 能够有效刻画生态系统的生理生态过程及其与外界环境的相互作用关系^[10]。Biome-BGC在国外已经被广泛用于气候变化的研究中^[11-14]。在我国, Biome-BGC模型较多应用于北方森林, 例如：祁连山^[15-16], 大兴安岭^[15], 天山^[17-19], 长白山^[20-22], 华北^[23]等地区。除了在北方的森林生态系统有所应用外, 我国学者在沼泽湿地^[5]、沙漠湿地^[24]、高寒草甸区^[25-26]以及热带(西双版纳)地区^[27]的研究也相继开展, 以上的研究成果均证明了模型具有很强的适应性。

树木年轮可以忠实地记录当地的气候变化, 利用树木轮宽数据来检验和验证Biome-BGC的NPP模拟能力, 分析特定森林生态系统的变化趋势的相关研究已有较多^{[19, 23, 28-31]}。Graumlich研究发现树轮宽度数据作为有力的证据可用来研究火山爆发后NPP的变化情况^[32]；Rathgeber在法国东南部的研究表明树轮宽度数据因可记录长期气候变化信息, 故可用于评估气候变化和CO₂对森林生长力的影响^[28]；在普罗旺斯的进一步研究显示了树轮宽度数据可用于验证区域尺度上BIOME3模拟的NPP数据^[29]。国内特别是华北地区、天山等地区也已有相关研究^{[19, 23, 30]}。苏宏新等使用Biome-BGC模型来研究天山云杉对未来气候变化的响应及其对CO₂浓度增加的响应, 并探讨了模拟NPP值与树轮宽度年表的关系^[19]；彭俊杰等对华北地区油松林生态系统研究发现树木年轮可检验校正Biome-BGC模型在内的各种生态系统过程模型^[23]。

太白山为秦岭的最高峰, 太白红杉主要分布于秦岭太白山一带, 林木生长缓慢, 对区域气候变化反应敏感, 可以较好地反应太白山的气候变化^[33-34]。在已有研究中, 林伟等采用野外调查、遥感影像解译、生物量计算相结合的方法对太白红杉碳储量进行了估算^[35]；李亮等应用LPS-GUESS模型模拟了太白山太白红杉林1958—2008年的净初级生产力, 但是此模型适用于全球或者区域尺度上, 在模拟结果上存在一定的不确定性^[36]。而Biome-BGC模型可以用于模拟不同尺度的生态系统能量、碳、氮、水分的循环^[10]。前人通过对太白红杉树木年轮的研究发现, 太白红杉树木年轮宽度对气候变化有很好的响应^{[34, 37]}。树木年轮宽度作为实测的参数, 可对模型模拟结果进行检验与校正^[23]。因此, 本研究将使用Biome-BGC模型模拟太白山太白红杉的植被净初级生产力(NPP), 并用树木年轮宽度及MOD17A3 NPP数据对模拟结果进行验证, 分析影响太白红杉树木生长的主要气候限制因子, 探讨南北坡及不同海拔高度太白红杉NPP对气候变化的差异响应。

1 研究区域与研究方法 1.1 研究区概况

太白山为秦岭主峰, 位于107°16′46″—107°56′28″E, 33°46′46″—34°12′1″N(图 1)。太白红杉在太白、周至、眉县三县均有分布。太白山北坡地势陡峻, 南坡相对较平缓^[38]。属大陆性季风气候区, 年平均气温为1.8—2.1℃, 降水量为800—900mm^[34], 植物资源丰富, 植被垂直分带明显^[39]。太白红杉分布区土壤类型主要为暗棕壤, 土层厚度在林线附近为30cm左右, 在分布下限地区为80cm左右^[40]。太白红杉作为二级国家保护植物, 是秦岭山区唯一生存的落叶松属渐危种, 主要分布在太白山一带2600—3600m的地区^[41]。因其生境特殊且日益濒危, 因此对太白红杉的研究可以为太白红杉的保护提供一定的参考价值。

1.2 年表建立及与气象数据的相关响应分析

太白山树木年轮于2013年分别在南北坡的不同海拔进行采集(图 1), 根据国际树轮数据库的标准(ITRDB), 每个采样点采集15—20棵树, 使用生长锥在树木的胸径处按照不同的方向取2个样芯。采集好的样芯进行晾干、固定、打磨^[42], 之后使用LINTAB(轮宽分析仪)进行轮宽测量, 并使用COFFCHA程序进行交叉定年及结果检测, 去除腐芯、破损、相关性差的样条序列^[43]。最后通过Arstan程序进行去趋势化和标准化, 采用负指数函数或线性函数进行拟合, 若不符合这两种函数则用样条函数进行拟合(样条函数的步长为序列长度的三分之二), 以保留轮宽序列中更多的低频信息, 采用双权重平均法合成树轮宽度年表^[44-45]。因标准年表(STD)中含有更多的低频信息, 所以本文采用标准年表进行研究。

为了探讨气候因子对太白红杉生长的影响, 利用树轮相关响应分析软件Dendroclim 2002对邻近气象站的气象数据(逐月)与树轮年表进行相关响应分析^[44]。南坡选取佛坪气象站(位于107°58′48″E, 33°31′12″N), 北坡则选取了眉县气象站(位于107°43′48″E, 34°16′12″N), 气象数据来源于陕西眉县和佛坪气象局。

1.3 Biome-BGC模型及数据

本研究使用美国Montana大学Numerical Terradynamic Simulation小组研发的Biome-BGC模型(4.2版本)(http://www.ntsg.umt.edu/project/biome-bgc)来模拟太白山太白红杉的NPP。Biome-BGC是一个基于日步长来模拟陆地生态系统植被及土壤中的碳、氮、水、能量流动和存储的生物地球化学循环模型^{[12, 46]}。Biome-BGC模型包括光合、蒸腾、呼吸、分解等作用以及光合产物的分配^[47]。

模型的输入文件包括初始化文件, 气象数据文件(日最高温、日最低温、日均温、日降水、饱和水汽压差、太阳辐射等)以及生理生态参数文件^{[5, 23]}。

初始化文件包括太白红杉采样点的位置、土壤参数、反照率等基本信息以及文件输入输出相关设定^{[12, 47]}。采样点基本位置及海拔如图 1所示, 土壤深度及土壤参数参考了已有研究文献中太白山南北坡的土壤参数^[48-49]。

气象数据文件采用眉县、佛坪气象站的气象数据(日最高温、日最低温、日降水量), 源于陕西眉县、佛坪气象局。并利用山地气候模型MT-CLIM4.3(http://www.ntsg.umt.edu/project/mtclim), 输入眉县、佛坪气象站的日最高温、日最低温、日降水量数据, 在考虑到海拔和坡度等因子的影响下, 模拟得到研究点的气象数据文件^[50]。近来研究得出, 太白山北坡气温直减率最大值为0.619℃/100m, 最小值为0.390℃/100m, 南坡最大为0.546℃/100m, 最小为0.449℃/100m, 本文在MT-CLIM4.3模型中使用了上文提到的气温直减率, 以提高MT-CLIM4.3模型的模拟精度^[51]。太白红杉的生理生态参数目前鲜有研究, Biome-BGC模型自带的参数是在大量文献及研究上得出的^[16], 故本文采用了模型默认的落叶松的参数(表 1)。

模型运行时首先进行初始化模拟(运行spin-up模式), 此模式采用工业革命前CO₂的浓度(294.842ppm), 氮沉积值为0.0001kg N m^-2 a^-1, 生态系统达到稳定状态后^[39], 运行normal-run模式得到研究点1960—2013年的NPP模拟值。

2 结果与分析 2.1 年表与气象因子的相关分析

太白山6个样点标准年表的统计结果如表 2所示。研究表明, 平均敏感度越大, 树木对外界环境因子的变化越敏感^[53]。太白山各采样点的平均敏感度在0.16—0.27之间, 反映了树木生长对环境响应较为敏感。样本总体代表性在0.89—0.94之间, 且信噪比均较高, 说明采样点树木生长受到外界干扰较少且年表质量较高^[54]。第一主成分所占的方差量均在40%以上, 说明树木生长能够较好的反映气候变化的一致性^[55]。标准年表的统计结果表明太白红杉树木年轮数据对区域气候变化响应敏感, 可作为植被生长的指标来验证Biome-BGC模型模拟的NPP结果。

树轮标准化年表与对应气象数据的相关分析结果表明(图 2), 在总的趋势上太白山大部分采样点与当年2—7月的气温呈正相关, 与当年4—6月的降水呈负相关, 且多个采样点达到了显著相关水平。生长季来临后, 气温的升高有助于促进太白红杉的生长, 使得生长期提前, 形成较宽的树轮, 而4—6月过多的降水则会造成气温下降, 抑制太白红杉的生长, 形成窄轮^[56]。从不同坡向差异来看, 北坡多个采样点年表与3—6月的气温相关系数通过了95%的置信度, 而南坡年表则是与5—7月的气温呈显著正相关。这可能是由于北坡较南坡接收的太阳辐射少, 初春时期树木生长对气温的响应更为敏感, 初春来临时, 北坡气温升高对太白红杉生长的促进作用要更为明显, 之后随着气温不断升高, 气温不再是限制太白红杉生长的主要因子, 故相关性逐渐降低^[57]。

2.2 NPP模拟结果的验证

模型运行时首先运行spin-up模式, 待生态系统达到稳定状态后运行normal-run模式, 最终得到NPP模拟值。Biome-BGC模型模拟的太白红杉NPP值显示, 1960—2013年北坡太白红杉NPP年均值为305.33g C m^-2 a^-1, 南坡为320.71g C m^-2 a^-1, 北坡略低于南坡。为了验证模拟结果的可靠性, 本文与其他研究结果进行了对比(表 3)。李亮等^[36]利用植被动态过程模型(LPJ-GUESS)模拟了太白山南北坡1958—2008年太白红杉的NPP, 结果表明太白红杉南北坡NPP分别为320g C m^-2 a^-1和430g C m^-2 a^-1, 南坡模拟值与本文模拟值结果一致, 但是北坡有所差异, 这可能是由于海拔差异引起的。李亮等^[36]选取了海拔3450m和3100m的研究点, 而本文除了在太白红杉的上限(3346.4m)、下限(3062m、3068m)选取研究点外, 在海拔3207m处也选取了研究点, 对北坡不同海拔高度的NPP值均进行了模拟。

同时, 本文还将模拟结果与分辨率为1km的卫星遥感MOD17A3 2000—2013年的NPP数据(https://ladsweb.modaps.eosdis.nasa.gov/)进行了对比。MOD17A3 NPP数据基于MODIS/TERRA卫星的遥感参数, 本文利用MRT(MODIS Reprojection Tool)软件对MOD17A3 NPP数据进行了拼接、重投影、数据格式转换、重采样等前期处理, 根据前人研究的太白红杉生长范围^[58]利用ArcGIS软件对其进行了裁剪, 得出太白红杉NPP分布图(图 3), 并将采样点的NPP进行提取, 与本文模拟的2000—2013年NPP值进行了比较(图 4)。从图 4看出, 本文模拟NPP值与MOD17A3 NPP数据除SBS5样点外, 整体趋势较为一致, 差值较小, SBS2、SBS3、YWD三个样点差值极小。这再次验证了Biome-BGC模型模拟的可靠性。SBS5样点对比值相差较大可能是由于两种尺度的不同所导致的。

为进一步验证模拟结果的准确性, 将采样点的树轮宽度指数与对应采样点的NPP模拟值进行了对比(图 5), 除SBS5、YWDZ两个采样点外, 其余采样点的树轮宽度指数与NPP模拟值的相关关系均呈显著相关(通过了95%的显著性), 其中, SBS1采样点呈极显著相关(通过了99%的显著性)。同时, 从图 5可以看出, 北坡采样点树轮宽度指数与NPP模拟值年际变化一致性要显著好于南坡, 多数采样点树轮宽度指数与NPP模拟值年际变化趋势趋于一致。这说明树轮宽度指数作为记录树木生长环境的重要指标, 可以较好的验证Biome-BGC模型模拟的NPP值。

2.3 NPP模拟结果的分析

上文验证了太白红杉NPP模拟结果的可靠性, 通过分析模拟结果(图 6), 按照十年尺度计算太白红杉NPP模拟值的变化趋势速率(表 4), 发现1960—2013年太白山太白红杉NPP总体呈波动上升趋势, 北坡的上升速率(0.47g C m^-2 a^-1)要小于南坡(1.29g C m^-2 a^-1)。1960—1979年太白红杉NPP值无明显变化, 南北坡均表现出相对平稳的变化趋势；1980—1989年南北坡NPP值均呈现出上升趋势, 且南坡上升速率高于北坡；1990—1999年为下降趋势且降幅较大, 北坡的下降速率高于南坡, 这可能是由于气候变化导致太白红杉的生态幅在此期间发生改变, 出现濒危现象^[57], 导致NPP值呈下降趋势；2000年以来南北坡NPP值表现出明显的上升趋势。

本文在选取采样点的时候, 分别在南北坡不同的海拔进行采样(表 5及图 6), 以研究太白红杉NPP值的海拔差异。结果显示北坡海拔3062m、3068m处的NPP值为428.01、416.33g C m^-2 a^-1, 海拔3207m处的NPP值为293.00g C m^-2 a^-1, 海拔3346.4m处的NPP值为83.98g C m^-2 a^-1。结果表明北坡太白红杉NPP模拟值随着海拔高度的上升而逐渐下降, 在太白红杉的下限NPP值最高, 随着海拔的升高, 生长环境发生变化, 气温、降水逐渐减少, NPP值随之降低。而南坡太白红杉NPP模拟值则随着海拔高度的上升有所增加。从不同海拔采样点NPP模拟值的年际变化来看(图 6), 南北坡1960—2013年NPP模拟值在不同海拔高度上变化趋势较为一致, 北坡低海拔的变化振幅要大于高海拔地区, 南坡无明显变化。

2.4 NPP模拟值与气象因子的关系及气候变化响应

气候变化会引起植物NPP发生变化, 为了更好地理解太白山太白红杉NPP值与气候的相关关系, 本文将每个采样点的模拟NPP值(月均值)与对应的气温、降水数据(北坡采样点采用眉县气象站数据, 南坡采用佛坪气象站)进行逐月相关分析, 结果如图 7所示。NPP值与气温的相关性显著好于降水, 且主要与太白红杉生长季的气温呈显著正相关, 说明气温的升高有助于NPP值的增加。北坡的NPP模拟值与初春气温呈正相关, 这与树轮宽度系数与NPP值的相关分析结果一致。与北坡不同的是, 南坡NPP模拟值与9月、10月的气温呈极显著正相关, 9月、10月的南坡气温的升高可以在一定程度上延长太白红杉的生长季, 从而使NPP值有所提高。

在全球增暖背景下, 已有研究证明秦岭地区近50年来气温逐渐上升, 降水日益减少^[59]。太白山太白红杉1960—2013年间NPP总体呈波动上升趋势, 说明气温的升高有利于NPP的增加。采样点树轮标准年表、NPP模拟值均显示出了与生长季气温的显著正相关, 全球气候变暖使得太白红杉生长季的气温增高从而导致NPP呈上升趋势；同时也显示出与生长季降水的负相关, 即降水的减少有利于太白红杉的生长, NPP的增加。总体来看, 气候变化作为制约太白红杉生境的重要因素, 影响太白红杉树木的生长, 进而对NPP的变化产生影响。

2.5 模型的适用性与不确定性

Biome-BGC模型可以模拟特定生态系统的生理生态过程及与外界因子的作用^[10]。树木年轮数据作为研究气候变化、树木生长生理参数的理想代用资料, 对Biome-BGC模型的模拟结果有一定的验证作用。且太白红杉作为秦岭特有的高山树种, 对气候变化响应极为敏感^[58], NPP及树木年轮的的变化均受到气候因子、环境因子、植被生理特征等多种因子的影响, 故树木年轮对植被的NPP具有一定的指示作用。本文对太白红杉树轮宽度指数与模拟的NPP进行相关分析表明, 多数样点具有很好的一致性, 达到了显著相关, 与个别样点的差异性可能是由于采集的树轮数量以及树木生长的个体差异性所导致的。Biome-BGC模型模拟的NPP值与MOD17A3 NPP数据的对比一致性较高, 个别差值较大与模型的不稳定性及两种数据的尺度不同有关。基于遥感数据的MOD17A3 NPP数据模拟尺度较大, 而Biome-BGC模型尺度较小, 模拟的是特定地区特定森林生态系统的NPP。总体而言, Biome-BGC模型模拟的太白红杉NPP, 结果较为可靠。

模型也有一定的不确定性, 主要是由于以下几个方面引起的：国内有关太白红杉的生理生态参数研究鲜有, 大部分生理生态参数只能使用模型自带参数, 使模拟结果的准确性有了一定的误差；高山气象站的器测气象数据时间较短甚至空白, 而高山植被对气候变化十分敏感, 本研究模型输入的气象数据虽在考虑到海拔影响下利用MT-CLIM4.3软件进行了处理, 但气象站之间的距离以及海拔差异仍可能影响NPP输出结果的准确性；Biome-BGC模型本身也有一定的局限性, 对森林生态系统的外在干扰活动(如风倒)等因素未考虑^[23]。模型参数以及植被生理生态参数的改进、气象数据的完善将会提高模型模拟的准确性。

由于实验条件的限制, 本项目组未开展太白红杉实测NPP、生理生态参数测定的相关实验, 以及对影响其生长的气候因子(水分、气温、CO₂等)及气候因子协同作用未进行水热条件协同分析、并分析其对生态系统的影响。导致模拟结果有一定的偏差, 这也为未来工作的开展提供了一定的方向。

3 研究结论

(1) 应用Biome-BGC模型模拟了太白山太白红杉林生态系统的NPP, 1960—2013年太白山太白红杉林北坡NPP年均值为305.33g C m^-2 a^-1, 南坡为320.71g C m^-2 a^-1, 北坡略低于南坡。1960—2013年太白山太白红杉NPP总体呈上升趋势, 北坡的上升速率(0.47g C m^-2 a^-1)要小于南坡(1.29g C m^-2 a^-1)。

(2) 太白红杉NPP模拟值对海拔有一定的响应关系。北坡太白红杉NPP模拟值随着海拔的上升而逐渐下降, 而南坡太白红杉NPP模拟值则随着海拔高度的上升有所增加。且南北坡1960—2013年NPP模拟值在不同海拔高度上变化趋势较为一致, 北坡低海拔的变化振幅要大于高海拔地区, 南坡无明显变化。

(3) 太白红杉多数采样点的树轮宽度指数与NPP模拟值呈显著相关, 树轮宽度指数作为记录树木生长环境的重要指标较好的验证了Biome-BGC模型模拟的NPP值。

(4) NPP模拟值与气象数据的相关分析表明：NPP值与气温的相关性显著好于降水, 且主要与太白红杉生长季的气温呈显著正相关即气温的升高有助于NPP值的增加。北坡的NPP模拟值与初春气温呈正相关, 南坡NPP模拟值与9月、10月的气温呈极显著正相关, 说明9月、10月的南坡气温的升高可以在一定程度上延长太白红杉的生长季, 从而使NPP值有所提高。