伏牛山龍池墁南坡油松徑向生長對氣候變化的響應

李 軒,彭劍峰,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指數間的定量關系。在目前全球氣候變化趨勢下,龍池墁山地低海拔地區將更加不利于幼齡油松的生長,隨著老齡油松不斷死去,最終會危及油松種群在低海拔地區生態群落中的地位,甚至造成龍池墁山地油松帶將向更高海拔遷移,增加該地區生態失衡的風險。