专题报道：自然保护地森林生态系统研究 || 本刊副主任委员姜姜教授团队开展浙江凤阳山针阔混交林碳动态模拟的研究

关键词：针阔混交林; Biome-BGC模型; 碳动态; 气候变化; 浙江凤阳山

基金项目：百山祖国家公园科学研究项目(2022JBGS03,2021ZDLY01);国家自然科学基金项目(32071612);江苏省碳达峰碳中和科技创新专项（BE2022307）。

引文格式：黄璐瑶,杜珊凤,纪小芳,等.基于 Biame-BGC模型的浙江凤阳山针阔混交林碳动态模拟[J].南京林业大学学报(自然科学版),2024,48(5):11-20.HUANG L Y, DU S F, Jl X F, et al. Garbon dynamic simulation based on Biome-BGC model in mixed coniferous and broadleaved forest of Fengyang Mountain, Zhejiang Province[J]. Journal of Nanjing Forestry University (Natural Sciences Edition),2024,48(5):11-20.DOI:10.12302/j.issn.1000-2006.202211005.

▲ BIOME-BGC模型来模拟浙江省凤阳山碳动态

浙江凤阳山自然保护区位于浙江龙泉(119°06'~119°15'E,27°46'~27°58'N)，南部与福建省接邻。保护区现有面积153 km²，其中林地4 245.2 hm²,森林覆盖率高达90.8%。气候特征为温暖湿润,雨水充沛，雾多。受季风影响，年降水量2 400 mm，年均气温和蒸发量分别为12.3 ℃和1 170 mm；相对湿度超过75%，年日照时长达1 448 h；自然保护区有丰富的光、热和水资源,非常适合动植物的生长繁殖。土壤类型主要为黄壤土,土壤结构主要呈现团粒结构，质地为黏壤土。风阳山自然保护区的主要植被类型为亚热带针阔混交林,包括许多天然分布的木本植物，如:光皮桦(Betula luminifera)、黄山木兰(Yulania cylindrica)、柏木(Cupressus funebris)、多脉青冈(Quercus multinervis)等，平均冠层高度约15 m，平均林龄约为40 a。人工杉木(Cunninghamia lanceolata)林、柳杉 (Cryptomeria japonica) 林多分布于海拔1 400~1 500 m处；在海拔600 ~1 300 m处有黄山松(Pinus taiwanensis)林分布；海拔300、900、1 500 m处均有木荷(Schima superba)分布。

1.1.2 模型选择

1)Biome-BGC模型。Biome-BGC模型是一种经典的生物地球化学模型，能够基于日尺度的气象数据对陆地植被、土壤碳、氮、水以及辐射通量进行模拟。此外，该模型考虑环境因子对全球气候变化和物质循环的影响，有着较高的认可度。本研究使用Biome-BGC 4.2版(Numerical Terradynamic Simulation Group， Missoula，MT， USA)。Biome-BGC模型的运行需要3个输入文件:初始化文件、气象数据文件和植被生理生态参数文件。初始化文件包括纬度、海拔及模型运行方式等基础属性数据；气象文件包括气象数据；植被生理生态参数包含44个参数。Biome-BGC模型运转从初始化模拟开始,利用初始化参数运转模型，使浙江凤阳山生态系统达到稳定状态，即碳输入和输出保持平衡。将初始化文件、植被生理生态参数文件和气候数据文件输入模型，模拟出该稳定状态下的植被生产力等参数，并且输出结果,与实测结果进行对比验证。

2)PEST模型。PEST模型是独立于模型外的非线性参数优化方法，基于Gauss-Marquardt-Levenberg算法求取目标函数(模型模拟值与实测值间的差异函数)的最小值。PEST模型首先识别模型中需要优化的参数后调用模型模拟,通过对比实际观测值与模型模拟值来确定参数的优化方向和优化值,使目标函数达到最优,其目标函数如下:

式中:y为m个(实测数据)元素的向量；x为含n个(参数个数)元素的向量；为模拟结果的向量；为待估算参数向量；H为m行n列的雅克比偏导矩阵；T代表转置符号；Q为具有m行m列的实测值权重矩阵。模型原理和使用详见操作指南(http://www.pesthomepage.org/)。

1.2.2 数据来源

海拔、纬度与土壤数据用于模型初始化模拟，生成稳定状态的初始化文件。纬度和海拔来源于浙江凤阳山通量塔站点(119°10'15″E,27°54'22″N)；土壤质地和土层深度数据来源于中国科学院资源环境科学数据中心平台(http://www.resdc.cn)。气温与降水数据来自CMFD，数据集中1979—2018年每3 h的气温与降水数据。运用MATLAB根据浙江凤阳山自然保护区坐标范围提取相应栅格数据,并运算出日累计降水量、日最高气温和日最低气温。随后将降水量、最高最低气温输入到MTCLM模型计算得到昼长、辐射以及VPD。大气CO₂浓度数据来源于CO₂数据库(https://www.co2.earth/historical-co2-datasets)。

植被生理生态参数有44个，其中有些用来描述站点植被类型以及特定植被类型下参数的缺省值(如描述植被是否为木本植物，是否为C₄或C₃植物等)，除去此类缺省参数，对其余参数采用PEST模型进行参数优化。PEST模型通过构建站点实测数据与模型模拟值间的优化函数,进而对 Biome-BGC 模型的生理参数进行优化。本研究样地针阔混交林中，阔叶树种占比超过80%，混合少量的杉木和黄山松。因此，采用Biome-BGC通用的阔叶树参数作为缺省值，在此基础上用PEST模型进行参数优化。参数优化结果参照表 1。

1.2.3 情景设计

根据IPCC第5次评估报告，未来气候变化将朝着增温、极端降水时间频发和大气CO₂浓度增加的趋势发展。在所有的RCP情景中，全球平均气温相对于1986—2005年期间均有不同程度的升高。在RCP2.6情景下，21世纪末全球气温将上升0.3~1.7 ℃；在RCP4.5 情境下，可能接近2 ℃；RCP6.0和RCP8.5情景下可能超过2 ℃。而降水对气候变化的响应是不均匀的，在干旱和潮湿地区之间的对比愈加强烈。

依据《第三次气候变化国家评估报告》，中国区域在20世纪末，相对于1986—2005年，RCP2.6情景的增温强度为每10 a增加0.08 ℃，降水强度为每10 a增加0.6%，且依据地球系统模式预测CO₂浓度到2100年约为421 μmol/mol；RCP4.5情景的增温强度和降水强度分别为每10 a增加0.26 ℃、1.1%，到2100年CO₂约为538 μmol/mol；而RCP6.0情境下的增温强度和降水强度分别为每10 a增加0.61 ℃、1.6%,到2100年的CO₂约为670 μmol/mol。本研究依据上述预测，以2011—2020年的气温、降水量平均值作为2021年的气温和降水量值，以线性增加的方式从2021年起每10 a增加0.06%，直到2100年，此情景为RCP2.6。RCP4.5、RCP6与RCP2.6增加方式一致，得到3种情景组合(表2)。

本研究利用一元线性回归方程对浙江省凤阳山1979—2018年的年总降水量、年平均辐射、年平均气温、净初级生产力、总初级生产力和净生态系统生产力开展时间趋势分析。为探究气象因子与净初级生产力(NPP)间的相互关系，运用MATLAB在日尺度上对4个气象因子(空气气温、降水量、饱和水汽压差、辐射)和NPP之间做二次函数拟合分析。NPP与各个月份平均气温和总降水量间的关系采用相关性分析进行，应用SPSS 25.0计算相关性是否显著(0.05水平)。

本研究应用纪小芳等提供的2017年通量塔实测GPP数据，对模型运行结果进行对比验证(图1)。

结果表明，模型对浙江凤阳山亚热带针阔混交林的GPP模拟值在通量塔实测值的范围内，并且变化趋势较为一致。本研究模拟的2017年GPP总量为1 471.26 g/(m²·a)，与纪小芳等测算凤阳山通量塔2017年实测总量(除去9月份)1 377.75 g/(m²·a)结果基本一致。同时，模拟结果与相近纬度其他亚热带森林生态系统如鼎湖山GPP[1 460.10 g/(m²·a)]和千烟洲针阔混交林生态系统GPP[1 669.7 g/(m²·a)]结果也较为相似。结合样本所在纬度和海拔等要素综合考虑，认为Biome-BGC模型的模拟结果具有一定的参考价值。

浙江凤阳山在过去40年的气候变化较为复杂，1979—2018年降水量、气温等变化见图2。线性拟合结果表明，浙江凤阳山多年降水量均值为169.54 cm，降水量最高值出现在2010年(236.10 cm)，最低值出现在1991年(107.20 cm)，降水量年际变化较大，总体呈现增加趋势。凤阳山多年平均气温为15.40 ℃,最高气温年份在2007年(16.45 ℃)，最低气温年份出现在1984年(14.20 ℃)，变异系数为3.78%，气温总体呈显著增长趋势(P<0.01)，线性拟合结果表明,其增温速率(0.32 ℃每10 a)较整个中国大陆的年均增温幅度(0.22 ℃每10 a)略高。凤阳山多年平均辐射量为298.96 W/m²，辐射最高的年份为2003年(346.49 W/m²)，最低为1990年(273.34 W/m²)。整体变化较为显著，呈现增加的趋势(P<0.01)。

2.3.1 凤阳山森林生态系统生产力年际变化

1979—2018年浙江凤阳山针阔混交林的生产力的动态变化见图3。GPP、NPP、NEP年际波动明显，变化特征基本一致。三者均存在3个波动下降期和3个波动上升期。三者的波动下降期为1979—1984、1986—1995和1998—2001年，上升期为1984—1986、1995—1998和2001—2018年。模拟时间段内，凤阳山森林GPP、NPP、NEP的极小值均在1984、1995和2001年，可能与这3个时期降水量增加、云量遮盖度增加、辐射量减少有关。而三者的极大值均在1998年，可能由于该年份的平均气温和平均辐射量都处于较高水平。

由图3可以看出，GPP和NPP的年际变化保持一致，先分别在1 392.94和451.25 g/(m²·a)处上下波动，直到2001年开始呈现逐渐上升趋势。1979—2018年的NPP波动范围为420.84~485.55 g/(m²·a)，变异系数为3.53%。最大值和最小值相差64.74 g/(m²·a)，平均值为451.25 g/(m²·a)。而1979—2018年的GPP平均值为1 392.94 g/(m²·a)，最高值出现在2018年[1 531.55 g/(m²·a)]，最低为2001年[1 306.39 g/(m²·a)]。由于凤阳山亚热带针阔混交林NEP波动显著，碳源/碳汇功能并不稳定，模拟年份在1984、2002、2005、2008及2010年出现5次负值，为碳源期，其他年份均表现为碳汇期(图3)。

2.3.2 凤阳山生产力日尺度变化特征

凤阳山模拟年份内的GPP、NPP日变化范围分别为1.55~6.31 g/(m²·d)和0.31~1.98 g/(m²·d)，NEP的日变化范围为-0.52~0.42 g/(m²·d)(图4)。

GPP、NPP和NEP三者40 d日平均值呈现极强的规律性。GPP、NPP前80 d平稳或上升缓慢、80~190 d大幅上升、190~230 d较为平稳、230~310 d大幅下降、310 d之后呈现平稳趋势,这与纪小芳等的通量塔研究结果较为一致。可以推断出:第80天和第310天为凤阳山针阔混交林生长开始和停止的时间节点。生长期从春季中期持续到秋季中期，这期间都处于碳生产状态。生长期内NPP整体变化趋势与GPP保持一致，但其变化幅度小于GPP。模拟年份内NEP均值日变化呈现多峰现象，日变化在数值上表现为“负—正—负—正—负”，说明凤阳山亚热带针阔混交林呈现出“碳源—碳汇—碳源—碳汇—碳源”的特征。NEP在夏季出现负值，可能由于夏季气温升高，促进土壤呼吸，异养呼吸部分排放的碳大于植被固定的碳，使凤阳山生态系统在此期间表现为碳源特征。

2.4.1 日尺度上NPP与环境因子的关系

为分析影响浙江凤阳山亚热带针阔混交林碳动态的影响因素，将1979—2018年的日NPP累计值与VPD、空气气温、辐射、和降水量进行二次函数拟合与显著性检验见图5。

由图5可知，4个气象因子均与NPP显著相关(P<0.01)。NPP与空气气温的相关性最高(R²=0.817 4),其次是饱和水汽压差(VPD)，而与降水量的相关性最微弱(R²≤0.001)。在-5~35 ℃的范围内,NPP随气温的增加而增加(图5b)；VPD与NPP呈现二次函数关系，VPD的值低于1 693 Pa时，NPP随VPD的增加而增加(图5a);辐射与NPP同样呈现二次函数关系,当日平均短波辐射小于523 W/m²时，NPP随辐射增加而增加，而后随辐射增加而处于降低趋势(图5c)。

2.4.2 月尺度上NPP与环境因子的关系

为进一步分析NPP对各个月份气象因子的响应程度，对NPP的年累计值与各个月份的平均气温和总降水量进行相关性分析。总体来看，浙江凤阳山针阔混交林的NPP与气温极显著正相关(R=0.902，P<0.001)、与降水量极显著负相关(R=-0.025，P<0.01)。在夏、秋两季具有较强的相关性，尤其是在4月和8月，而在春冬两季部分月份与气温呈现出负相关(图6a)。与气温相比，凤阳山NPP与降水的相关性较弱，但是与9月降水量呈现较明显的正相关(P<0.01)。4月和8月都呈现一定程度的负相关(图6b)。可能原因是降水超过一定额度便会引起气温和太阳辐射的下降，此时的气象因子不利于凤阳山针阔混交林的生长，尤其不利于光合作用的进行。所以NPP在4月和8月与降水呈现负相关。

本研究模拟RCP2.6、RCP4.5、RCP6.0情境下的碳动态变化，结果见图7。RCP情景到21世纪末大气都处于增温增湿、二氧化碳浓度增加的状况，只是增加的程度不同，这3个RCP情境下，浙江凤阳山的植被总生产力都处于不断增加的趋势。RCP6.0增加较为缓和，到后期处于平稳波动的状态，而RCP4.5增加幅度较RCP2.6增加幅度大，RCP6.0即增温强度和降水强度分别为每10 a增加0.61 ℃、1.6%，到2100年的CO₂约为670 μmol/mol，浙江凤阳山在该情境下GPP增加最为强烈。截至21世纪末，RCP2.6、RCP4.5和RCP6.0情景下凤阳山GPP分别达到1 552.73、1 660.30及1 960.41 g/(m²·a)，相比2018年的GPP分别增加1.38%、8.41%和28.00%。

凤阳山森林生态系统GPP、NPP在1979—2018年内总体呈现增加趋势，除2002、2005、2008及2010年，其余年份为碳汇，且呈现“碳源—碳汇”季节交替的特征。该区碳通量特征受水热条件的限制，其中受气温的限制作用最强。且在不同时间尺度上，受限制程度也各不相同。此外，不同气候限制条件下，凤阳山森林总初级生产力在21世纪均呈现上升趋势，RCP6.0的增加幅度大于RCP2.6和RCP4.5。

Biome-BGC模拟结果显示，凤阳山森林年平均NEP在正常情况下大于零，生态系统表现为碳汇。这与纪小芳等的研究结果相同，即凤阳山森林生态系统在2017年表现为碳汇。而模拟年份内NEP均值日变化呈现双峰现象，日变化在数值上表现为“负—正—负—正—负”，NEP在夏季出现负值，这与张越等的研究略有不同，秦岭地区太白红杉NEP在年内呈现单峰，究其原因是凤阳山地区夏季气温过高，且降水量充沛，会促进土壤呼吸，异养呼吸部分排放的碳大于植被固定的碳，使凤阳山生态系统在此期间表现为碳源特征。

气象因子是影响森林生长的重要因素之一，对森林碳通量特征的年际变化起着重要作用。气象因子与NPP间的相关性分析结果显示，气温是影响森林NPP积累的主要限制因子，认为气温升高会在一定程度上促进NPP的积累，而降水量与NPP的相关性最为微弱，这与前人研究结果相印证。张凤英等利用线性回归分析法分析认为，长江流域森林NPP与气温之间呈现正相关，且相关性比降水与森林 NPP 的相关性强；Liu 等研究发现喀斯特地区年NPP与年平均气温呈显著正相关，而与年降水量无显著相关；李亮等通过模型模拟发现气温升高有利于太白红杉生长。

但该结果仅体现在日尺度和年尺度上，然而具体到月尺度气象因子与NPP的相关关系，认识到该地区森林生长对气候的响应具有季节差异性。夏季和冬季凤阳山森林生态系统对气候的响应不同。夏季气温与NPP呈现正相关，夏季降水量与NPP之间呈现负相关，这与前人的研究结果一致，江南及江淮地区夏季降水和NPP之间呈现较强负相关，夏季气温和NPP的相关关系与降水相反，可能因为浙江省夏季降水较为充沛，减弱了浙江省地区太阳到达地表的有效辐射，从而降低了光合作用的强度；夏季是我国植被生长的主要季节，如果气温适当偏高，光合作用强度增大，则会加快NPP的积累速度。此外，凤阳山冬季降水量与NPP之间呈现正相关，冬季气温与NPP间却呈现负相关。究其原因:NPP的大小等于GPP减去植被呼吸,纪小芳等研究表明植被呼吸随气温的升高而加强。所以，冬季气温上升一方面会促进植被呼吸加强，从而减少NPP的积累；另一方面，气温上升光合作用增强，但同时植被的蒸腾作用随之增强，为减少植物水分的过分损失，植物气孔关闭，这又会导致植物吸收固定二氧化碳的能力减弱,从而减少NPP的积累。这可能是冬季气温与NPP呈现负相关的主要原因。而江淮地区在冬季气候相对干燥，降水是冬季植被生长的主要限制因子，一定程度的降水增加有利于NPP的积累。

同时，应认识到未来气温、降水量和二氧化碳浓度同时增加的不同情境下，凤阳山地区森林总初级生产力均呈现上升趋势，但RCP6.0的增加幅度大于RCP2.6和RCP4.5，说明森林碳动态和3个气象因子耦合作用间的关系为正反馈。这与张越等、李传华等研究结论一致，气温、降水量和二氧化碳浓度同时增加的组合下，GPP呈现上升趋势。值得注意的是，本研究虽运用Biome-BGC模型很好地揭示了凤阳山过去、现在及未来的碳动态特征及其对气候变化的响应，结果具有一定可信性，但也存在较大的局限性，首先模型的参数较多，则模拟结果的不确定性就会增大；其次，假设未来气候条件的不同情景可能会与实际气候有所出入，从而对预测结果造成误差，并且由于目前对浙江凤阳山针阔混交林植被生理生态参数的研究极少，利用PEST优化算法优化的参数范围仍需进一步考究。

综上所述，本研究采用Biome-BGC模型对凤阳山森林生态系统40 a的碳交换进行模拟，分别在日尺度和月尺度上对NPP与气象因子的关系进行探究，并预测未来气候变化下凤阳山的碳动态变化趋势和范围，得出主要结论如下:

1)凤阳山亚热带针阔混交林1979—2018年的平均GPP、NPP和NEP分别为1 392.94、451.25、16.21 g/(m²·a)。1984、2002、2005、2008及2010年表现为“碳源”期，其他年份均表现为碳汇期。年内呈现“碳源—碳汇—碳源—碳汇—碳源”的特征。

2)凤阳山生态系统的NPP对气温的变化最为敏感，并且在一定程度内，NPP随着气温、辐射、VPD的增大而增大；而与降水量的相关性最微弱，且在一定程度内，NPP随降水量的增加而减小。但是与空气气温、辐射、VPD和降水量均显著相关。

3)夏季气温的上升对NPP的增加起积极作用，而冬季气温的升高却对NPP起到反作用。一定程度内，冬季降水量对NPP有促进作用，而夏季降水对NPP为负作用。

4)RCP2.6、RCP4.5、RCP6.0未来气候情景变化下，截至2100年，GPP均呈现增加的趋势，增加的幅度有所差异，其中RCP6.0的增加幅度最大。且至2100年，RCP2.6、RCP4.5和RCP6.0情景下凤阳山GPP分别达到1 552.73、1 660.30及1 960.41 g/(m²·a)，相对于2018年GPP分别增加1.38%、8.41%和28.00%。