Evaluation of Carbon Sequestration Capacity of a 57-year-old Korean Pine Plantation in Mt. Taeh wa based on Carbon Flux Measurement Using Eddy-covariance and Automated Soil Chamber System

에디 공분산 및 자동화 토양챔버 시스템을 이용한 탄소 플럭스 관측 기반 태화산 57년생 잣나무조림지의 탄소흡수능력 평가

Abstract

Forests are the largest carbon (C) sinks in terrestrial ecosystems. Recently, as enhancing forest C sequestration capacity has been proposed as a basic direction of the Republic of Korea's "2050 Carbon Neutral Strategy," accurate estimation of forest C sequestration has been emphasized. According to the Intergovernmental Panel on Climate Change guidelines, sequestration quantity is calculated from changes in C stocks in forest C pools, such as biomass, deadwood, litter and soil layer, and harvested wood products. However, in Korea, only the overstory biomass increase is now considered the amount of sequestration quantity, so there can be a significant difference from the actual forest C sequestration. In this study, we quantified forest C exchange through C flux measurement using an eddy covariance system and an automated soil chamber system in a 57-year-old Korean pine plantation located in Mt. Taehwa, Gwangju-si, Gyeonggi-do. Then, the net amount of C sequestration was compared with the amount of the overstory biomass increase. We estimated the annual C stock change in the remaining C pools by comparing the net sequestration amount from the C flux measurement with the overstory biomass increase and C stock change in the litter layer. Therefore, the net C sequestration of the Korean pine plantation estimated from the flux measurement was 5.96 MgC ha^-1, which was about 2.2 times greater than 2.77 MgC ha^-1 of the overstory biomass increase. The annual C stock increase in the litter layer was estimated to be 0.75 MgC ha^-1, resulting in a total annual C stock increase of 2.45 MgC ha^-1 in the remaining C pools. Our results indicate that the domestic forest is a larger C sink than the current methods, implying that more accurate calculations of the C sequestration capacity are necessary to quantify C stock changes in C pools along with the C flux measurement.

산림은 육상생태계에서 가장 큰 탄소흡수원으로 기후변화 대응에 있어 산림의 대기 중 이산화탄소 농도 저감 역할은 중요하다. 최근 '2050 탄소중립 계획'에 산림의 탄소흡수 기능의 강화가 기본 방향으로 제시되면서 정확한 산림의 탄소흡수량 산정이 강조되고 있다. 산림부문의 탄소흡수량은 Intergovernmental Panel on Climate Change 지침을 따라 산림 내 생물량, 고사목, 임상 유기물층, 토양층, 수확된 목재제품 등 여러 탄소 저장고 내 탄소축적 변화량으로부터 산정한다. 그러나 국내 산림의 경우 하층 식생을 제외한 주요 수종의 임목 재적 증가로부터 추정한 생물량 증가량만을 산림의 탄소흡수량으로 산정하고 있어 실제 산림의 탄소흡수량과 큰 차이가 발생할 수 있다. 이에 본 연구에서는 경기도 광주시 태화산에 위치한 57년생 잣나무 조림지에서 에디 공분산 시스템과 자동화 토양챔버 시스템을 이용한 탄소 플럭스 관측을 통해 산림의 탄소 교환량 및 순 탄소흡수량을 정량화하고, 이를 현재 산림의 탄소흡수량 산정 방법에 따라 법정림 임분수확표 내 연평균 생장량과 국가 고유계수를 이용하여 계산한 잣나무 조림지 임목의 생물량 증가량과 비교하였다. 또한 탄소 플럭스 관측기반의 순 탄소흡수량과 잣나무 조림지의 생물량 증가량 및 임상 유기물층의 탄소저장 변화량 등의 차이로부터 나머지 탄소 저장고에서 연간 탄소저장 변화량을 추정하였다. 그 결과 탄소 플럭스로부터 계산한 잣나무 조림지의 연간 순 탄소흡수량은 5.96 MgC ha^-1으로 생물량 증가로부터 계산한 임목의 연간 탄소흡수량 2.77 MgC ha^-1보다 약 2.2배 많았다. 연간 임상 유기물층의 탄소저장 변화량은 0.75 MgC ha^-1로 추정되어, 연간 하층 식생, 고사목, 토양층 등의 탄소 저장고로 유입되는 탄소의 양이 2.45 MgC ha^-1으로 추정되었다. 본 연구의 결과는 국내 산림이 현재 평가 수준보다 더 큰 탄소흡수원임을 보여주며, 탄소 플럭스 관측과 더불어 하층 식생, 고사목, 토양층 등의 탄소 저장고에서 탄소축적 변화량의 정량화를 통해 더욱 정확한 산림부문 탄소흡수량 산정이 필요함을 시사한다.

서론

산림은 연간 50 Pg 이상의 탄소를 흡수하는 육상생태계에서 가장 큰 탄소흡수원이다(Beer et al., 2010). 유엔기후 변화협약에서는 산림을 주요 탄소 저장고로 인정하고 있으며(IPCC, 2014), 기후변화 대응에 있어 산림의 대기 중 이산화탄소(CO₂) 농도 저감 역할이 강조되고 있다. 최근 전 세계적으로 기후변화의 심각성을 인식하고, 2050년까지 탄소 중립을 달성하기 위한 대책 수립에 나서고 있다. 우리 정부도 ‘2050 탄소 중립’ 목표 달성을 위한 정책을 추진하면서 산림의 탄소흡수 기능의 강화를 기본 방향으로 제시하였다(ME, 2020). 세부 내용에 따르면 국내 산림의 노령화에 따라 산림의 탄소흡수량이 2017년 기준 연간 4,570만 톤에서 2050년에 약 30% 수준까지 감소할 것으로 전망되며, 산림 복원 및 조림 등을 통한 탄소흡수원의 확대와 산림경영을 통한 탄소흡수능력의 강화 등을 세부 전략으로 제시하고 있다. 이에 따라 앞으로의 산림 경영활동에 따른 산림의 탄소흡수량 변화를 정확히 추정하고 예측하기 위해서는 그 어느 때보다도 산림의 탄소흡수량 추정 방법에 대한 정확한 검증이 필요한 상황이다.

산림의 순 탄소흡수량은 총 탄소흡수량과 총 탄소 배출량의 차인 순생태계생산량(NEP, Net Ecosystem Production)으로 나타낼 수 있는데, NEP가 0보다 크면 탄소흡수, 0보다 작으면 탄소 배출을 의미한다. 산림의 NEP를 추정하는 방법은 크게 3가지로 나눌 수 있다(Jiang et al., 2020). 첫 번째는 산림 내 각 탄소 저장고에서의 탄소저장 변화량의 합(\(\begin{aligned}\sum_{i=1}^{n} \Delta C_{i}\end{aligned}\))을 이용하여 계산하는 방법(식 1)으로, 탄소 저장고는 크게 상층식생 생물량(C_overstory), 하층식생 생물량(C_understory), 고사목(C_deadwood), 임상 유기물층(C_litter), 토양층(C_soil) 등으로 구분할 수 있다(Jiang et al., 2020). 두 번째는 산림 내 각 구성요소의 순일차생산량(NPP, Net Primary Production)의 합과 종속생물호흡량(R_h, Heterotrophic Respiration)의 차이를 이용하여 계산하는 방법(식 2)이며, 세번째는 산림의 총 탄소흡수량을 나타내는 총일차생산량(GPP, Gross Primary Production)에서 산림생태계 내 각 호흡을 통한 탄소 배출량의 총합, 즉 생태계 호흡량(R_e, Ecosystem Respiration)을 제외하는 방법(식 3)이다.

NEP = ΔC_overstory+ΔC_understory+ΔC_deadwood+ΔC_litter+ΔC_soil+ΔC_others       (1)

NEP = NPP_overstory+NPP_understory+NPP_others-R_h       (2)

NEP = GPP – R_overstory – R_understory – R_h

= GPP – R_above – R_s = GPP – R_e       (3)

기후변화에 관한 정부 간 협의체(IPCC, Intergovernmental Panel on Climate Change)에서는 첫 번째 방법(식 1)과 같이 산림의 탄소 저장고를 지상부 생물량(C_above), 지하부 생물량(C_below), 고사목(C_deadwood), 임상 유기물층(C_litter), 토양층(C_soil), 수확된 목재제품(C_{harvested wood product}) 등 6가지로 구분하여 \(\begin{aligned}\sum_{i=1}^{6} \Delta C_{i}\end{aligned}\)으로부터 산림의 탄소흡수량을 계산하도록 제시하고 있다(Aalde et al., 2006; Ogle et al., 2019).

우리나라도 IPCC의 지침에 따라 임분수확표 내 산림의 주요 수종별 지위지수 및 임령에 따른 정기평균 생장량 또는 연평균 생장량 자료(NIFoS, 2020)와 국가 고유계수(목재 기본밀도, 생물량 확장계수, 뿌리 함량비 등)를 이용하여 C_above 및 C_below의 탄소저장 변화량을 산정하고 있다(NIFoS, 2019). 그러나 C_above 및 C_below의 탄소저장 변화량 산정 시 C_understory는 반영되지 않고 있으며, C_deadwood, C_litter, C_soil 등은 탄소저장 변화량에 관한 연구 활동 자료가 충분히 구축되지 않은 상황이기 때문에 연간 탄소저장량의 변화가 0이라고 가정하고 있다(Aalde et al., 2006; NIFoS, 2015; Ogle et al., 2019; GIR, 2020). 따라서 현행 방법에 따라 산정된 산림의 탄소흡수량과 실제 산림의 탄소흡수량 간에 큰 차이가 발생할 수 있지만, 연간 ΔC_deadwood, ΔC_litter, ΔC_soil 등을 추정하기 위해 낙엽 낙지 생산량 및 고사목과 임상 유기물의 분해량 등을 높은 빈도로 계속해서 측정하는 것은 간단하지 않으며, 환경 요인에도 큰 영향을 받기 때문에 정확한 예측이 쉽지 않다(Saarsalmi et al., 2007; Bradford et al., 2016; Petraglia et al., 2019; Cheng et el., 2020). 또한 토양 내 유기 탄소의 형성과 분해 과정은 세근 생산 및 고사의 생물계절학적 특성과 미생물의 활동, 환경 요인 등의 상호작용에 영향을 받고, 공간적으로도 큰 이질성을 보이기 때문에 산림 내 정확한 ΔC_soil를 추정하는 데 많은 어려움이 있다(Lehmann et al., 2007; Cordon et al., 2010; Luan et al., 2021).

에디 공분산 시스템과 자동화 토양챔버 시스템을 이용한 산림의 탄소 플럭스 측정은 산림-대기 및 산림 내 탄소 이동을 측정하여 식 3과 같이 NEP를 추정할 수 있다. 에디 공분산 시스템은 기체분석기와 3차원 풍향풍속계를 이용하여 대기와 산림 간의 CO₂ 교환량을 측정함으로써 광합성을 통한 총 탄소흡수량인 GPP, 호흡을 통한 총 탄소 배출량인 R_e, 그리고 GPP와 R_e의 차인 NEP를 계산할 수 있다(Kang et al., 2014, 2018; Chang et al., 2021). 또한, 자동화 토양챔버 시스템은 챔버가 닫혀있는 동안 챔버 내부의 CO₂ 농도변화를 이용하여 뿌리 호흡 및 유기물의 분해 등을 포함한 토양표면으로부터의 탄소 배출량을 계산할 수 있다(Sha et al., 2021). 두 가지 방법의 혼용은 산림-대기 및 산림 내 탄소 이동을 비교적 짧은 시간 규모(초단위)에서 연속적으로 측정할 수 있어 직접 \(\begin{aligned}\sum_{i=1}^{n} \Delta C_{i}\end{aligned}\)를 측정하지 않고도 환경 요인의 변화에 따른 산림생태계의 생산 및 분해 반응과 그에 따른 탄소흡수 및 배출량을 보다 정밀하게 분석할 수 있다(Baldocchi, 2003; Chang et al., 2021; Sha et al., 2021).

따라서 본 연구에서는 경기도 광주시에 있는 태화산 6영급 잣나무(Pinus koraiensis) 조림지에서 탄소 플럭스 측정을 통해 산림의 탄소흡수능력을 평가함으로써 산림의 탄소흡수량 추정 방법의 개선 방향을 제시하고자 하였다. 본 연구의 세부적인 목표는 (1) 에디 공분산 및 자동화 토양챔버 시스템을 이용한 탄소 플럭스 측정을 통해 산림의 탄소흡수 및 배출량을 정량화하고, (2) 현행 방법에 따른 ΔC_overstory로부터 계산한 잣나무 조림지 임목의 탄소흡수량과 탄소 플럭스 측정으로부터 추정한 산림생태계 탄소흡수량을 비교하여, (3) 각 탄소순환의 구성요소에 대한 고찰을 바탕으로 산림의 탄소흡수량 추정 시 개선 가능한 부분을 제언하는 것이다.

재료 및 방법

1. 조사지 개황

본 연구는 경기도 광주시 도척면 상림리에 있는 서울대학교 태화산학술림 내 잣나무 조림지(북위 37° 18′17′′, 동경 127° 19′02′′)에서 수행하였다(Figure 1). 기상청 이천 종관기상관측 자료 기준 해당 지역의 최근 30년간(1991~2020년) 연평균기온은 11.4±0.7°C이며, 연평균강수량은 1,316±308 mm이었다. 본 연구가 수행된 잣나무 조림지는 해발고도 188~196 m 내 태화산 동북 사면 산복에 있으며, 사면의 경사는 10.0±0.8°이다. 서울대학교 태화산학술림 산림조사부 자료에 따르면 본 조사지는 1964년에 조림된 57년생 잣나무 조림지로 1988년과 2000년에 두 차례에 걸쳐 간벌 작업이 시행되었으며(Ryu et al., 2014), 현재 평균 임분 밀도는 약 456 tree ha^-1이다. 조사지 내 잣나무는 임분 전체 기저 면적의 90% 이상을 차지하고 있으며, 임분 내 하층 식생으로는 산철쭉과 개옻나무가 주를 이루고 있다(Lee et al., 2019). 잣나무의 평균 수고는 약 21 m이며, 흉고직경은 32.5±6.1 cm이다. 조사지의 토양은 약산성의 사양토 또는 양토이다(Suh et al., 2014; Lee et al., 2021).

Figure 1. Research site at an approximately 57-year old Korean pine plantation on the Mt. Taehwa Seoul National University Forest, Gyeonggi-do, Republic of Korea.

2. 에디 공분산 시스템 기반의 CO₂ 플럭스 관측

본 연구 조사지인 잣나무 조림지의 GPP, R_e, 그리고 NEP를 계산하기 위해 2019년 11월 1일부터 2020년 10월 31일까지 에디 공분산 플럭스 타워 관측자료를 활용하였다. 주된 관측자료는 플럭스 타워의 41 m 높이에서 3차원 풍향 풍속계(CSAT-3, Campbell Sci., Inc., USA)와 적외선 가스 분석기(LI-7500A, Li-Cor, Inc., USA)를 이용하여 10 Hz 빈도로 관측한 풍향, 풍속, CO₂ 농도 자료이다. 관측 장비는 거칠기 저층(roughness sublayer) 높이를 고려하여 난류 혼합에 의한 물질 교환이 활발한 내부 경계층(internal boundary layer) 내 41 m 높이에에 설치되었다(NIER, 2011). 이와 동시에 대기 온도 및 상대습도, 강수량, 토양온도 및 토양수분 등의 30분 평균 산림 기상 자료를 취득하였다. 관측 장비 및 높이에 대한 자세한 내용은 Table 1에 나타내었다. 관측된 자료는 KoFlux 표준화 자료 처리방법에 따라 10단계에 걸쳐 자료 처리 및 품질 관리를 실시하였다(Kang et al., 2018; Lee et al., 2021). 먼저 EddyPro 소프트웨어(Li-Cor, Inc., USA)를 이용하여 3차원 풍향 · 풍속 자료의 좌표 변환(평면 맞추기 회전법, planar fit rotation) (Wilczak et al., 2001)과 공기밀도 변동에 따른 플럭스 교정(Webb et al., 1980) 등을 진행하였다. 이후 MATLAB 프로그램(The MathWorks, Inc., USA)을 이용하여 플럭스 자료의 이상치 제거, 산림 기상 자료의 결측값 보간, 주간 CO₂ 플럭스의 결측값 보간(주변 분포 표본 추출법, marginal distribution sampling)(Reichstein et al., 2005), 야간 CO₂ 플럭스 보정(Lloyd and Taylor, 1994)과 NEE의 R_e 및 GPP로의 배분 작업 등을 수행하였다. 산림 기상 자료의 결측값 보간은 조사지 인근의 기상청 광주 및 용인 방재 기상 관측 자료와 이천 종관기상관측 자료 중 일조시간, 기온, 해면 기압, 증기압, 풍향 · 풍속, 강수량 등의 자료와의 선형회귀식을 이용하여 수행하였다. 야간 CO₂ 플럭스 보정은 Kang et al.(2014)과 같이 광반응곡선법(Light Response Curve)(Lee et al., 1999), 마찰속도(Friction Velocity, u*) 필터링법(Falge et al., 2001; Gu et al., 2005), 이류의 영향을 고려한 van Gorsel법(van Gorsel et al., 2007, 2008, 2009) 3가지를 이용하여 수행하였고, 각 방법에서 추정된 GPP, R_e, NEP의 평균값을 분석에 이용하였다.

Table 1. List of instruments for measurements of CO₂ flux and environmental factors at an approximately 57-year old Korean pine plantation on the Mt. Taehwa Seoul National University Forest, Gyeonggi-do, Republic of Korea.

^a Li-Cor, Inc., USA, ^b Cambell Sci., Inc., USA, ^c Kipp & Zonen, Netherlands, ^d Vaisala, Finland, ^e Texas Electronics, USA.

3. 자동화 토양챔버 시스템 기반의 토양 CO₂ 방출량 측정

잣나무 조림지에서 토양 CO₂ 방출량을 측정하기 위해 30 m×30 m 크기의 조사구 내에 24시간 연속 측정 가능한 자동화 토양챔버(automated open/close chamber, AOCC) 시스템을 구축하였다. 본 연구의 AOCC 시스템은 Suh et al.(2014)과 같이 챔버부, 펌프부, 전력분배부 세 가지 주요 부분으로 구성하였다. 총 12개의 챔버가 10분 간격으로 순차적으로 열리고 닫히면서 챔버 내 CO₂ 농도변화를 기록하며, 모든 챔버가 측정을 완료하는 한 사이클에 총 2시간이 소요된다. 본 연구의 AOCC 시스템은 closed dynamic chamber(CDC) 방법을 기반으로 하며, 밀폐된 챔버에서 단위 시간 동안 CO₂ 농도변화를 가스분석기(LI-850, Li-Cor, Inc., USA)를 이용하여 측정한 후, 다음의 식에 따라 토양 CO₂ 방출속도(R_s; mg CO₂ m^-2 h^-1)를 계산하였다(Bekku et al., 1995; Suh et al., 2006, 2014).

R_s = 60⋅10^-6⋅a⋅ρ⋅V⋅A^-1       (4)

여기서 a는 단위 시간 동안 챔버 내 CO₂ 농도 변화율(ppm min^-1)을 의미하며, ρ는 대기 중 CO₂ 밀도(mg m^-3), V는 챔버의 부피(0.3×0.3×0.2 m³), A는 챔버와 맞닿는 토양 단면적(m²)을 나타낸다. 본 연구에서는 에디 공분산 시스템과 동일한 관측 기간에 계산된 R_s 자료를 분석에 이용하였다. R_s 자료의 결측값 보간은 지온 및 토양수분과의 상관관계로부터 추정한 비선형회귀식을 이용하여 수행하였다(Yuste et al., 2003). 먼저 지온과 R_s 간의 4가지 비선형회귀식을 도출한 후, R²_adj와 RMSE값을 기준으로 최적의 지온-R_s 회귀식을 선정하였다(Figure S1, Table S1). 이후 선정된 지온-R_s 회귀식을 이용하여 지온의 영향을 최소화한 R_{s_normalized}(실제 관측한 R_s/지온으로부터 추정한 R_s)을 계산하였다. 계산된 R_{s_normalized}를 이용하여 토양수분과 R_{s_normalized}간의 2가지 비선형회귀식을 도출하였으며, 동일하게 R²_adj와 RMSE 값을 기준으로 최적의 토양수분-R_{s_normalized} 회귀식을 선정하였다(Figure S1, Table S1). 최종적으로 최적의 지온-R_s 회귀식과 토양수분-R_{s_normalized} 회귀식을 종합하여 지온 및 토양수분 변화에 따른 R_s 예측 모델을 Figure S2와 같이 도출하였으며, 본 모델을 이용하여 R_s 자료의 결측값(366일 중 57일, 15.6%) 보간을 수행하였다.

 

Figure S1. (a) Response of soil respiration(R_s) to soil temperature(T_s) and (b) response of normalized soil respiration(R_{s_normalized}, the ratio of measured R_s to estimated R_s from T_s) to volumetric soil water content(Ɵ) at an approximately 57-year old Korean pine plantation on the Mt. Taehwa Seoul National University Forest, Gyeonggi-do, Republic of Korea.

Table S1. Adjusted R² and RMSE values of four nonlinear models to estimate response of soil respiration(R_s) to soil temperature(T_s) and of two nonlinear models to estimate responses of normalized soil respiration(R_{s_normalized}, the ratio of measured R_s to estimated R_s from T_s) to volumetric soil water content(Ɵ) at an approximately 57-year old Korean pine plantation on the Mt. Taehwa Seoul National University Forest, Gyeonggi-do, Republic of Korea.

 

^a van't Hoff (1884) equation : R = ae^bT

^b Arrhenius (1889) equation : R = ae^{{b/R*[1/(Tref+273)-1/(T+273)]}}

^c Lloyd and Taylor (1994) equation : R = ae^{{b[1/(Tref+46.02)-1/(T+46.02)]}}

^d Ratkowsky et al. (1982) equation : R = (T – T_min)²/(T_ref – T_min)²

^e Doran et al. (1990) equation : R = aƟ² + bƟ + c

^f Mielnick and Dugas (2000) equation : R = a(Ɵ - Ɵ_min)(Ɵ_max - Ɵ)^b

 

Figure S2. Modeled soil respiration(R_{s_Modeled}) from soil temperature(T_s) and volumetric soil water content(Ɵ) vs measured soil respiration(R_{s_Measured}). The equation of R_{s_Modeled} is R_{s_Modeled} = (aƟ² + bƟ + c)((T_s – d)²/(e – d)²). Orange lines represent 95% prediction interval.

4. 잣나무의 생물량 증가량(ΔC_overstory) 기준 산림의 탄소흡수량 추정

현재 국내에서 산림의 임목 탄소흡수량을 계산하는 방법에 따라 잣나무 조림지의 연간 생물량 증가량, 즉 ΔC_overstory를 연간 산림의 탄소흡수량(tC ha^-1 yr^-1)으로 추정하였다(NIFoS, 2019, 2020).

ΔC_overstory = V⋅D⋅BEF⋅(1+R)⋅CF       (5)

여기서 V는 연평균 생장량(m³ ha^-1 yr^-1), D는 목재 기본 밀도(t d.m m^-3), BEF(Biomass expansion factor)는 생물량 확장계수, R은 뿌리 함량비, CF(Carbon fraction)는 탄소전환계수이다. V는 본 연구 조사지의 평균 수고 21 m를 기준으로 NIFoS(2020)의 잣나무 법정림 임분수확표 내 지위지수 16과 임령 55년에 해당하는 연평균 생장량 값을 이용하였으며, D, BEF, R, CF 값은 NIFoS(2019)의 잣나무 국가 고유계수 및 탄소전환계수 값을 이용하였다.

5. 임상 유기물층의 탄소저장 변화량(ΔC_litter) 추정

임상 유기물층으로의 연간 탄소 유입량을 추정하기 위해 조사구(30 m×30 m) 내에 10개의 낙엽수거망을 설치하여 2020년 10월부터 12월까지 매주 낙엽수거망 내 떨어진 낙엽 낙지 등을 수거하였다. 수거한 유기물은 70°C에서 항량이 될 때까지 건조시킨 후, 건중량에 CF를 곱하여 본 연구 조사지에서의 연간 유기물생산량으로 가정하였다. 연간 유기물분해량의 경우 2020년 이전의 유기물생산량 자료가 부재하여 2020년의 유기물생산량과 선행 연구의 잣나무낙엽 분해 상수(k)를 이용하여 아래의 식에 따라 간접적으로 추정하였다(You et al., 2000; Won et al., 2018b).

연간 유기물분해량 = 연간 유기물생산량 – (EXP(– k)⋅연간 유기물생산량) (6)

최종적으로 연간 유기물생산량과 유기물분해량의 차이를 ΔC_litter이라고 가정하고 분석을 수행하였다.

6. 자료 분석

본 연구 조사지에서의 탄소순환을 Figure 2와 같이 간단히 도식화하였다. Figure 2와 식 3에 근거하여 각 탄소 플럭스가 전체 산림의 탄소 교환량에서 기여하는 정도를 확인하기 위해 연간 GPP, R_e, NEP, R_s 값을 정량화하고, 각 탄소 플럭스의 GPP에 대한 기여도를 확인하였다. 또한 R_e와 R_s의 차이를 지상부호흡량(Aboveground Respiration, R_above)이라고 가정하여 R_above가 GPP에서 차지하는 비중도 함께 확인하였다. 실제 산림의 탄소흡수량과 현행 방법으로부터 추정한 산림의 탄소흡수량의 차이를 확인하기 위해 플럭스 관측으로부터 추정한 NEP와 잣나무의 생물량 증가량, ΔC_overstory의 차이를 확인하였다. 식 1과 식 3의 관계로부터 아래와 같이 식 7의 관계를 도출하여 NEP와 ΔC_overstory 및 ΔC_litter의 차이로부터 나머지 탄소 저장고에서의 연간 탄소저장 변화량을 추정하였다.

Figure 2. Schematic carbon cycle of research site at an approximately 57-year old Korean pine plantation on the Mt. Taehwa Seoul National University Forest, Gyeonggi-do, Republic of Korea. Diamond-shaped boxes are gross primary production(GPP), autotrophic respiration(R_a), net primary production(NPP), heterotrophic respiration (R_h), and net ecosystem production(NEP) for the Korean pine plantation. Oval-shaped boxes are carbon stock change in carbon pools including aboveground overstory biomass pool(ΔC_{overstory_above}), belowground overstory biomass pool(ΔC_{overstory_below}), aboveground understory biomass pool(ΔC_{understory_above}), belowground understory biomass pool(ΔC_{understory_below}), deadwood pool(ΔC_deadwood), litter layer(ΔC_litter), and soil layer(ΔC_soil). Carbon pool and flow such as volatile carbon emission, dissolved organic carbon, soil methane uptake, insect are omitted from the cycle.

ΔC_understory+ΔC_soil+ΔC_deadwood+ΔC_others = NEP–ΔC_overstory–ΔC_litter       (7)

본 연구의 모든 통계적 분석은 R(ver. 4.1.1, R Core Team)과 SigmaPlot(ver. 12.5, SPSS Inc. USA) 소프트웨어를 이용하여 수행하였다.

결과

1. 산림 기상 관측 결과

2019년 11월 1일부터 2020년 10월 31일까지 연구 조사지의 평균 기온은 11.9±9.9°C이었으며, 총 강수량은 1,568 mm이었다. 같은 기간 연구 조사지 내 평균 지온은 12.4±7.1°C이었으며, 평균 토양수분함량은 19.2±5.2%이었다. 일평균 기온은 8월 초에 최고 28.6°C를, 12월 말에 최저 –9.3°C를 기록하였으며, 일평균 지온의 경우 8월 중순에 최고 24.8°C를, 2월 초에 최저 2.0°C를 기록하여 기온보다 계절에 따른 변화폭이 작은 것으로 확인되었다[Figure 3(a)]. 일최다 강수량은 133 mm로 장마철인 8월 초에 확인되었으며, 일평균 토양수분함량은 장마철 중 많은 비가 내렸던 8월 초에 최고 49.0%를, 장마가 시작되기 전인 6월말에 최저 9.8%를 기록하였다[Figure 3(b)].

Figure 3. (a) Daily mean air temperature(T_a) and soil temperature(T_s), (b) daily total precipitation(Prec) and daily mean volumetric soil water content(Ɵ), and (c) daily total gross primary production(GPP), net ecosystem production(NEP), ecosystem respiration(R_e), soil respiration(R_s), and modeled soil respiration (R_{s_modeled}) at an approximately 57-year old Korean pine plantation on the Mt. Taehwa Seoul National University Forest, Gyeonggi-do, Republic of Korea.

2. 산림의 탄소 교환량 및 NEP

일 GPP는 겨울철 특히 1월경 일 최저 –0.3 gC m^-2 d^-1를 기록하였으며, 2월부터 꾸준히 증가하여 7월 중순에 일최대 12.7gC m^-2 d^-1를 기록하였다[Figure 3(c)]. 많은 비가 내렸던 8월 초에 일 GPP가 2.2 gC m^-2 d^-1까지 일시적으로 크게 감소하였으나, 장마철 이후 다시 증가하였다가 9월 부터 서서히 감소하는 경향을 보였다. 일 R_e는 GPP와 마찬가지로 2월부터 꾸준히 증가하였으나, GPP와 달리 여름 장마 기간에 계속해서 증가하여 8월 초에 일최대 7.5 gC m^-2 d^-1를 기록한 후, 서서히 감소하는 경향을 보였다. 일 R_s도 R_e와 비슷한 패턴을 보였으나, 지온이 10°C 이상이었던 5월에서 10월 사이 연간 R_s의 약 79%가 대기 중으로 배출되는 것으로 나타났다. 결과적으로 일 NEP는 GPP와 비슷한 경향을 보였으나 GPP가 감소하고 R_e가 크게 증가한 장마 기간에 연중 최저치인 –4.8 gC m^-2 d^-1까지 감소하여 일시적인 탄소 배출원이 되기도 하였다[Figure 3(c)].

월 GPP는 1월에 최저 31.6 gC m^-2, 6월에 최대 268.0 gC m^-2을 기록하였다(Table 2). 그러나 월 NEP는 월 R_e/GPP 비가 46%로 가장 작았던 4월에 최대 102.6 gC m^-2을 기록하였으며, 월 R_e/GPP 비가 102%였던 8월에는 –3.9 gC m^-2로 탄소흡수량이 배출량보다 작았다. 월 R_e는 2월부터 꾸준히 증가하여 8월에는 최대 206.3 gC m^-2를 보였으며, 마찬가지로 월 R_s도 2월부터 꾸준히 증가하여 8월에 최대 101.2 gC m^-2(R_e의 49%)를 배출하였다. 그러나 R_e와 R_s의 차로부터 추정한 월 R_above는 8월보다 이른 6월에 최대 127.7 gC m^-2(R_e의 73%)를 기록하여 R_s와 서로 다른 패턴을 보였다.

Table 2. Montly total gross primary production(GPP), ecosystem respiration(R_e), soil respiration(R_s), aboveground respiration(R_above, the difference between R_e and R_s), and net ecosystem production(NEP) from Nov. 1^st 2019 to Oct. 31^st 2020 at an approximately 57-year old Korean pine plantation on the Mt. Taehwa Seoul National University Forest, Gyeonggi-do, Republic of Korea. Total amount of carbon flux during this study period and contribution of each carbon flux component to GPP are also shown below.

2019년 11월 1일부터 2020년 10월 31일까지 본 연구대상지인 잣나무 조림지의 총 탄소흡수량인 GPP는 1,905.6 gC m^-2로 그중 69%인 1,308.9 gC m^-2이 R_e를 통해 다시 대기 중으로 배출되는 것으로 나타났다(Table 2). 그 중 R_s를 통한 탄소 배출이 GPP의 27%인 512.5 gC m^-2으로, R_e와 R_s의 차로부터 추정한 R_above는 GPP의 42%인 796.4 gC m^-2이었다. 따라서 약 1년간 본 잣나무 조림지의 순탄소흡수량인 NEP는 GPP의 31%인 596.5 gC m^-2이었으며, 단위 변환 시 5.965 MgC ha^-1이었다.

3. 잣나무 조림지의 ΔC_overstory 및 ΔC_litter

법정림 임분수확표 내 지위지수 16의 55년생 잣나무 조림지의 V는 5.98 m³ ha^-1 yr^-1이었으며(NIFoS, 2020), D는 0.408 t d.m. m^-3, BEF는 1.742, R은 0.276, CF는 0.51을 적용하였다(NIFoS, 2019). 그 결과 ΔC_overstory는 2.766 MgC ha^-1이었으며, 단위 변환 시 GPP의 약 15%인 276.6 gC m^-2으로 플럭스 관측으로부터 추정한 산림의 순 탄소흡수량 NEP의 46%에 해당하였다.

한편 2020년 10월부터 12월까지의 유기물생산량은 총 111.7 gC m^-2(GPP의 약 6%)이었으며, k에 0.4(You et al., 2000; Won et al., 2018b)를 적용하여 식 6에 따라 계산한 연간 유기물분해량은 GPP의 약 2%인 36.9 gC m^-2이었다. 이에 따라 유기물생산량과 분해량의 차이로부터 추정한 연간 ΔC_litter은 74.8 gC m^-2(GPP의 약 4%)이었으며, 식 7로 부터 추정한 연간 ΔC_understory+ΔC_soil+ΔC_deadwood+ΔC_others의 총합은 GPP의 약 13%인 245.1 gC m^-2으로 ΔC_overstory의 89% 수준이었다.

고찰

본 연구가 수행된 57년생 잣나무 조림지에서 계절별 탄소 교환량 분석 결과, 대체로 기온이 높은 여름철 생장 기간에 GPP, R_e, NEP 모두 증가하는 경향을 보였다[Figure 3(c)]. 그러나 많은 비가 내렸던 8월에는 6월 대비 GPP가 24%(65.5 gC m^-2) 감소한 반면, R_e는 18%(31.5 gC m^-2) 증가하여 8월 NEP는 6월 대비 97.1 gC m^-2 감소한 –3.9 gC m^-2로 일시적인 탄소 배출원이 되었다(Table 2). 선행 연구에 따르면 동아시아 지역의 몬순 기후로 인한 여름철 집중 강수 및 태풍의 발생은 낙엽량을 증가시키거나 일사량을 감소시켜 연간 산림생태계의 GPP를 최대 30%까지 감소시키는 것으로 보고된 바 있다(Saigusa et al., 2010; Hong and Kim, 2011). 본 연구의 에디 공분산 플럭스 타워에서 관측한 기상 자료에 따르면 8월 강수량은 677 mm로, 7월 강수량 323 mm의 약 2.1배, 6월 강수량 97 mm의 약 7배의 많은 비가 내렸으며, 강수일수는 6월에 13일에서 8월에 20일까지 약 54%가 증가하였다. 이에 따라 월평균 일사량은 6월에 228.2 W m^-2에서 8월에 127.0 W m^-2까지 약 44% 감소하여 광합성을 통한 총 탄소흡수량인 GPP가 감소한 것으로 생각된다. 반면 8월 R_e는 6월 대비 증가하였는데, 세부적으로 8월 R_above는 6월 대비 18%(22.6 gC m^-2) 감소하였지만, R_s는 6월에 47.1 gC m^-2에서 8월에 101.2 gC m^-2까지 약 2.1배 증가하여 R_e의 증가는 R_s의 증가에서 기인한 것으로 나타났다.

2016년부터 2018년까지 본 연구가 수행된 57년생 잣나무 조림지에서의 선행 연구 결과에서도 많은 비가 집중적으로 내리는 장마철에 생태계 호흡량이 크게 증가하였다(Lee et al., 2021). 그 원인으로 장마철 이전 세근 등의 뿌리 생장 및 활동이 낮은 토양수분함량으로 억제되었다가 장마에 의한 많은 강수로 토양수분함량이 회복되어 수분과 양분의 획득을 위해 크게 증가하기 때문인 것으로 고찰하였다(Park et al., 2010; Hagedorn et al., 2016; Wang et al., 2016). 실제로 본 연구에서 토양수분함량이 약 20% 이하로 감소할 때 R_s가 유의하게 감소하는 경향을 보였으며(Figure S1), 선행연구에서는 잣나무의 세근 생산량이 장마철과 그 이후인 7월에서 11월 사이에 크게 증가하는 것으로 확인되었다(Park et al., 2010; Zhou et al., 2014; An et al., 2017). 따라서 본 연구에서 R_above가 6월 이후 완만하게 감소한 반면, R_s 및 R_e가 6월 대비 8월에 증가한 원인 중 하나는 장마가 시작되기 전 6월 말에 9.8%까지 감소하였던 일평균 토양수분함량이 7~8월에 많은 강수로 회복되면서 잣나무가 수분 및 양분을 확보하기 위해 지하부의 세근 생산 및 활동을 크게 증가시켜 지하부 생장호흡(growth respiration)이 증가하고, 증가된 생물량으로 인한 유지호흡(maintenancec respiration)이 함께 증가하였기 때문인 것으로 생각된다[Figure 3(b), (c), Table 2]. 뿐만 아니라, Hagedorn et al.(2016)은 크게 감소했던 토양수분의 회복 이후 세근 주변의 균근 및 미생물 생물량과 토양호흡량이 크게 증가한 것으로 보고하였는데, 본 연구에서도 7~8월 세근 생장과 더불어 주변 균근 및 미생물 군집의 성장으로 인한 생장 및 유지호흡의 증가가 R_s의 증가에 영향을 미쳤을 것으로 생각된다.

본 연구에서 57년생 잣나무 조림지의 연간 총 탄소흡수량 GPP는 1,905.6 gC m^-2이었으며, 그중 69%가 호흡을 통한 총 탄소 배출량 R_e를 통해 대기 중으로 다시 배출되어 연간 순 탄소흡수량은 GPP의 31%인 596.5 gC m^-2이었다. 아시아 지역의 22개 산림생태계에서 연간 탄소 교환량 자료를 이용하여 메타 분석한 결과에 따르면, R_e/GPP 비(%)는 일반적으로 60% 이상이었으며, 잣나무 조림지와 같은 상록수림에서는 약 50~100%의 비를 보였다(Chang et al., 2021). 그중 열대림과 같은 저위도 지역을 제외한 온대 및 아한대 산림의 경우 R_e/GPP 비가 평균 약 80% 정도에 이르는 것으로 나타났다(Kondo et al., 2017). 이에 따라 아시아 지역 상록수림의 NEP/GPP 비는 0~50% 범위를 보였으며, 본 연구 조사지인 잣나무 조림지의 R_e/GPP 및 NEP/GPP 비는 아시아 지역의 분포 범위 내에 부합하였다.

본 연구 조사지인 잣나무 조림지와 비슷한 국내외 소나무속 조림지와 침엽수림 및 혼효림 등의 선행 연구 조사 결과(Table 3), 국내 약 20~40년생 참나무-잣나무 혼효림의 연 NEP는 362 gC m^-2이었다(Lee et al., 2017). 일본의 44년생 편백 조림지의 연 NEP는 477~480 gC m^-2이었으며(Takanashi et al., 2005), 90년생 소나무림의 경우 388 gC m^-2이었다(Mizoguchi et el., 2012). 중국의 200년 이상의 잣나무 우점 혼효림의 경우 연 NEP가 308 gC m^-2이었으며(Zhang et al., 2006), 그 외 유럽과 미주 등 소나무속 조림지의 연 NEP는 임령에 따라 24 gC m^-2부터 765 gC m^-2까지 다양하게 분포하는 것으로 나타났다(Kolari et al., 2004; Vickers et al., 2012; Ziemblińska et al., 2016; Chan et al., 2018). Table 3의 선행 연구 결과를 바탕으로 임령, 연평균강수량, 연평균기온에 따른 연 NEP의 변화 양상을 확인하였다. 먼저 임령과 연 NEP 간의 상관관계에서는 임령이 약 40년이 될 때까지 연 NEP가 가파르게 증가하였으며, 73년 이후 서서히 감소하였다[Figure 4(a), R²_adj = 0.57, p < 0.05]. 그러나 연 NEP와 연평균강수량 및 연평균기온간에 양의 상관관계를 보였으나 통계적으로 유의하지는 않았다[Figure 4(b), p = 0.47; Figure 4(c), p = 0.0516].

Table 3. Site information and annual net ecosystem production (NEP, gC m^-2 yr^-1) for conifer and/or mixed forests(no. 1-8, eddy covariance-based NEP estimation; no. 9-10, biometric-based NEP estimation).

Figure 4. Relationship between (a) forest age and annual net ecosystem production(NEP), (b) mean annual precipitation and annual NEP, and (c) mean annual temperature and annual NEP based on data from Table 3(no. 1-8 with black circles, eddy covariance-based NEP estimation; no. 9-10 with grey squares, biometric-based NEP estimation).

하지만 생태계 천이단계 및 임령에 따른 연 NEP의 변화를 확인한 Gough et al.(2016)의 연구 결과에서는 초기에 변이가 컸으나, 대체로 임령이 증가함에 따라 연 NEP가 서서히 감소하는 결과를 보였다. Yu et al.(2014)의 결과에서도 임령의 증가에 따라 연 NEP가 꾸준히 감소하는 패턴을 보였으며, Gao et al.(2016)에서는 임령과 더불어 연평균기온이 연 NEP 증가에 큰 영향을 미치는 것으로 나타났다. 이로부터 임분의 갱신 초기에는 높은 GPP와 갱신 방법에 따라 많은 차이가 발생하는 R_e에 의해 NEP의 변이가 크지만, 대체로 갱신 40년 이후에는 임령이 증가함에 따라 NEP가 감소하는 경향을 확인할 수 있었으며, 향후 산림의 정확한 탄소흡수량을 추정하기 위해서는 임령과 더불어 위도별 기온의 차이나 기후변화로 인한 기온 증가의 영향을 고려해야 할 것으로 생각된다.

본 연구에서 임분수확표 내 지위지수 16, 임령 55년 잣나무 임분의 연평균 생장량과 국가 고유계수를 이용하여 계산한 연간 ΔC_overstory는 276.6 gC m^-2이었으며, Ryu et al.(2014)의 연구에서 국내 잣나무 조림지를 대상으로 제작된 지상부 생물량 추정 모델로부터 계산한 277.8 gC m^-2과 유사한 수준이었다. 또한 ΔC_overstory는 플럭스 관측으로 부터 추정한 산림의 연간 순 탄소흡수량 NEP 596.5 gC m^-2의 약 46% 정도였다. 두 방법으로부터 계산한 탄소흡수량 차이는 319.9 gC m^-2이었으며, 그 중 ΔC_litter의 74.8 gC m^-2를 제외한 연간 ΔC_understory+ΔC_soil+ΔC_deadwood+ΔC_others의 총합은 245.1 gC m^-2로 NEP의 41% 수준이었다. 선행 연구에 따르면, 소나무속 임분의 연간 유기물생산량은 대체로 10월에서 11월에 집중되었다(Kinerson et al., 1974; NamGung and Mun, 2009; Jeong et al., 2013). 그러나 낙엽 활엽수림의 연간 유기물생산량의 약 85% 이상이 가을철에 집중되는 것과 달리(Kwon et al., 2016), 소나무속과 같은 상록침엽수림의 경우 가을철 이외의 기간에 연간 유기물생산량의 약 40% 이상이 발생하였다(Kinerson et al., 1974). 이에 따라 본 연구에서 ΔC_litter 계산 시 이용한 연간 유기물생산량은 10월부터 12월에 한정적으로 측정하여 임상 유기물층으로 유입되는 탄소량을 저평가할 수 있으며, 또한 연간 유기물분해량의 경우에도 본 연구 기간 이전에 유입된 유기물의 분해량을 고려하지 않았기 때문에 본 연구의 ΔC_litter에는 다소 불확실성이 존재한다.

미국 North Carolina 지역의 약 15년생 Loblloly pine 임분의 연간 유기물생산량은 약 462 g m^-2으로 탄소량으로 계산(본 연구와 동일한 CF=0.51 적용) 시 약 235.6 gC m^-2의 탄소가 임상으로 유입되었다(Kinerson et al., 1974). 이와 동일하게 CF=0.51 값을 적용할 경우, 핀란드 내 약 34개 Scots pine 임분의 연간 유기물생산량은 12.2~167.4 gC m^-2의 범위를 보였으며, 본 연구의 잣나무조림지와 임령이 비슷한 62년생 Scots pine 임분은 연간 최대 155.1 gC m^-2의 탄소가 임상으로 유입되었다(Starr et al., 2005). 국내의 경우, 월악산 국립공원 내 소나무 군락의 연간 유기물생산량은 2005년부터 2007년 사이 평균 157.0 gC m^-2(NamGung and Mun, 2009), 2009년부터 2013년 사이 평균 175.5 gC m^-2이었으며(Won et al., 2018a), 남산 소나무 군락의 2008년과 2009년 사이 연간 평균 유기물생산량은 360.6 gC m^-2이었다(Jeong et al., 2013). 따라서 본 연구의 연간 유기물생산량 111.7 gC m^-2을 국내 월악산 및 남산 소나무 군락의 평균 유기물생산량인 231.0 gC m^-2 yr^-1으로 대체하여 ΔC_litter값을 재계산할 경우 기존의 74.8 gC m^-2에서 154.9 gC m^-2까지 증가하여 연간 ΔC_understory+ΔC_soil+ΔC_deadwood+ΔC_others의 총합은 기존 245.1 gC m^-2(NEP의 약 41%)에서 165.0 gC m^-2(NEP의 약 28%)까지 감소할 수 있을 것으로 추정된다.

본 연구 조사지와 같은 잣나무림 또는 잣나무 혼효림에서의 연간 세근 생산량은 0~30 cm 깊이의 표토에서 약 180 gC m^-2까지 달하였고, 0~90 cm 깊이를 고려할 경우 약 212 gC m^-2까지 이르는 것으로 나타났다(Park et al., 2010; Zhou et al., 2014; Wang et al., 2016). 또한, 표토에서 연간 세근 고사량 및 토양으로의 유입량이 최소 50 gC m^-2에서 최대 253 gC m^-2까지 이르러 본 연구에서 ΔC_litter가 저평가될 수 있는 점을 감안하더라도 연간 ΔC_understory+ΔC_soil+ΔC_deadwood+ΔC_others의 총합인 245.1 gC m^-2만큼 상당한 것으로 나타났다(Park et al., 2010; Kim et al., 2012; Zhou et al., 2014; An et al., 2017). 특히 현재 임목의 탄소 흡수량 추정 시 이용되는 연간 생물량 증가량 계산은 특정 기간 사이의 토양 내 뿌리 함량의 변화량을 바탕으로 추정하는데, 연간 세근 생산 및 고사량의 계절에 따른 큰 변이를 고려하였을 때 실제 세근 생산량 및 고사에 의한 토양으로의 탄소 유입량이 저평가되었을 수 있다(Park et al., 2010; Zhou et al., 2014; Wang et al., 2016; An et al., 2017). 따라서 정확한 산림의 탄소흡수량 추정을 위해서는 시공간적으로 높은 해상도의 지하부 탄소순환에 대한 자료 구축 및 추정 연구가 필요할 것으로 생각된다. 나아가 자동화 토양챔버 시스템 구축과 지하부의 뿌리 및 미생물 호흡의 구분 연구를 통해 지하부 탄소 배출량의 환경 요인에 따른 변화를 정밀 분석하고, 지하부 뿌리 및 미생물 활동에 의한 정확한 탄소 배출량을 추정할 수 있을 것이다.

결론

본 연구에서는 에디 공분산 및 자동화 토양챔버 시스템을 이용한 탄소 플럭스 측정을 통해 산림의 탄소흡수 및 배출량을 정량화하고, 현행 방법에 따른 주요 상층 수종의 생물량 증가량으로부터 계산한 임목의 탄소흡수량과 탄소 플럭스 측정으로부터 추정한 탄소흡수량을 비교하여 산림의 탄소흡수량 추정 시 개선 가능한 부분을 제언하고자 하였다. 경기도 광주시 태화산에 있는 57년생 잣나무 조림지는 연간 약 19.05 MgC ha^-1의 탄소를 광합성 활동을 통해 흡수하였으며, 그 중 13.09 MgC ha^-1을 호흡을 통해 배출하여 연간 순 탄소흡수량은 약 5.96 MgC ha^-1이었다. 반면 현재 국내에서 산림의 탄소흡수량 추정 시 이용하고 있는 주요 상층 수종의 생물량 증가량으로부터 추정한 지위지수 16, 임령 55년의 잣나무 조림지의 임목의 탄소흡수량은 2.77 MgC ha^-1으로 실제 산림생태계 탄소흡수량의 약 46% 정도를 차지하였다. 그 차이는 3.20 MgC ha^-1으로 그 중 ΔC_litter의 0.75 MgC ha^-1를 제외한 연간 C_understory+ΔC_soil+ΔC_deadwood+ΔC_others의 총합은 2.45 MgC ha^-1이었다. 향후 산림 경영활동 및 기후변화에 따른 정확한 산림의 순 탄소흡수량을 추정하기 위해서는 탄소 플럭스 관측 방법을 통해 현행 방법을 보완할 수 있을 것이며, 산림 기상 및 임령에 따른 탄소 교환량 추정 또한 가능할 것으로 생각된다. 나아가 현재 산림의 탄소흡수량 추정 방법에 포함되고 있지 않은 하층식생과 토양층의 탄소순환에 대한 정량화 및 추정 연구에 탄소 플럭스 관측을 적용하여 더 정확한 산림의 탄소흡수량을 추정할 수 있을 것이다.

감사의 글

본 연구는 산림청 산림과학기술 연구개발사업 (교란생태계 적응 및 복원기술 개발 사업, 2020185D10-2122-AA02)의 지원으로 수행되었으며, 본 연구를 도와주신 서울대학교 학술림 직원분들께 진심으로 감사드립니다.