DOI QR코드

DOI QR Code

Future Prospects of Forest Type Change Determined from National Forest Inventory Time-series Data

시계열 국가산림자원조사 자료를 이용한 전국 산림의 임상 변화 특성 분석과 미래 전망

  • Eun-Sook, Kim (Forest Ecology Division, National Institute of Forest Science) ;
  • Byung-Heon, Jung (Forest Policy and Economics Division, National Institute of Forest Science) ;
  • Jae-Soo, Bae (Future Forest Strategy Department, National Institute of Forest Science) ;
  • Jong-Hwan, Lim (Forest Ecology Division, National Institute of Forest Science)
  • 김은숙 (국립산림과학원 산림생태연구과) ;
  • 정병헌 (국립산림과학원 산림정책연구과) ;
  • 배재수 (국립산림과학원 미래산림전략연구부) ;
  • 임종환 (국립산림과학원 산림생태연구과)
  • Received : 2022.01.06
  • Accepted : 2022.10.11
  • Published : 2022.12.31

Abstract

Natural and anthropogenic factors cause forest types to continuously change. Since the ratio of forest area by forest type is important information for identifying the characteristics of national forest resources, an accurate understanding of the prospect of forest type change is required. The study aim was to use National Forest Inventory (NFI) time-series data to understand the characteristics of forest type change and to estimate future prospects of nationwide forest type change. We used forest type change information from the fifth and seventh NFI datasets, climate, topography, forest stand, and disturbance variables related to forest type change to analyze trends and characteristics of forest type change. The results showed that the forests in Korea are changing in the direction of decreasing coniferous forests and increasing mixed and broadleaf forests. The forest sites that were changing from coniferous to mixed forests or from mixed to broadleaf forests were mainly located in wet topographic environments and climatic conditions. The forest type changes occurred more frequently in sites with high disturbance potential (high temperature, young or sparse forest stands, and non-forest areas). We used a climate change scenario (RCP 8.5) to establish a forest type change model (SVM) to predict future changes. During the 40-year period from 2015 to 2055, the SVM predicted that coniferous forests will decrease from 38.1% to 28.5%, broadleaf forests will increase from 34.2% to 38.8%, and mixed forests will increase from 27.7% to 32.7%. These results can be used as basic data for establishing future forest management strategies.

우리나라 산림의 임상은 자연적·인위적 요인에 의해 지속적으로 변화하고 있다. 임상(침엽수림, 활엽수림, 혼효림)면적의 비율은 국가 산림자원 특성 파악에 중요하게 활용되는 정보이기 때문에 임상 변화에 대한 정확한 이해와 전망이 필요하다. 따라서 본 연구에서는 국가산림자원조사 시계열 자료를 이용하여 임상 변화 발생 특성을 이해하고 이를 기반으로 미래 임상 변화 예측치를 도출하는 것을 목표로 하였다. 제5차, 제7차 국가산림자원조사 자료의 10년 기간 임상 변화정보와 임상 변화에 영향을 미칠 수 있는 변수(기후, 지형, 임분, 교란 등)를 이용하여 임상 변화 특성을 분석한 결과, 우리나라 산림은 침엽수림이 감소하고 혼효림과 활엽수림이 증가하는 방향으로 변화하고 있는 것으로 확인되었다. 침엽수림에서 혼효림으로, 혼효림에서 활엽수림으로 변화되는 지역은 주로 지형적으로 습윤하고 강수량이 많아서 수분관련 생육환경이 양호하며 주변에 활엽수림이 많은 지역이었다. 또한 기온이 높은 지역, 임분의 임령과 밀도가 낮은 지역, 주변 지역에 비산림이 많은 지역 등 교란 가능성이 높은 지역에서 변화가 많이 발생했다. 이러한 임상의 변화 특성을 반영하여 기계학습 모형(SVM)을 구축하고 기후변화시나리오(RCP 8.5)를 이용하여 미래의 임상 변화를 전망한 결과, 2015년에서 2055년까지 40년 동안 침엽수림은 38.1%에서 28.5%로 감소, 활엽수림은 34.2%에서 38.8%로 증가, 혼효림은 27.7%에서 32.7%로 증가할 것으로 예측되었다. 본 미래 임상분포 변화 정보는 향후 산림관리 전략 수립의 기초자료로 활용될 수 있다.

Keywords

서론

산림생태계는 기상 · 기후환경, 입지 · 토양환경, 동 · 식물이 상호작용하면서 지속적으로 변화하는 동적시스템이다. 산림 내에서는 다양한 수종들이 생명활동을 하면서 생장하고 경쟁하고 일부는 도태된다. 나무가 타 수종과의 경쟁을 이겨내고 큰 나무로 성장하여 숲의 수관부를 차지하는 우세목이 되면 임분 내 수직 구조와 생태계 프로세스를 결정하는데 가장 큰 영향력을 발휘한다. 따라서 우리나라에서는 광역적인 산림의 특성을 이해하고 구분하기 위한 가장 큰 분류체계로 수관부를 점유하는 수종의 종류 및 비율에 따라 임분의 유형, 즉 임상(Forest Type)을 침엽수림, 활엽수림, 침활혼효림으로 구분하고 있다. 2020년 기준, 우리나라의 전국 산림은 침엽수림 39%, 활엽수림 33%, 침활혼효림 28%로 구성되어 있다(Korea Forest Service, 2021)1).

임상은 고정되어 있지 않고 시간 경과에 따라 지속적으로 변화하고 있다. 과거 우리나라 산림의 침엽수림 면적 비율은 1980년 52%, 1990년 49%, 2000년 43%, 2010년 42%로 지속적으로 감소해왔다. 또한 활엽수림 면적 비율은 1980년 19%, 1990년 22%, 2000년 27%, 2010년 28%로 지속적으로 증가해왔다(Kim et al., 2015). 이렇게 전국 산림의 임상 구성이 변화한다는 것은 산림의 특성을 기반으로 산림을 보호, 관리, 육성해야 하는 산림관리정책의 기본 바탕이 지속적으로 변화하고 있다는 것을 의미한다. 따라서 현재와 미래의 임상분포 현황과 변화 특성을 잘 파악하고 미래 변화상을 전망하는 일은 전국 산림의 관리체계를 세우는 데 매우 중요한 일이라 할 수 있다. 이런 이유로 국가산림정책의 목표와 방향을 설정하는 산림기본계획은 산림자원 및 임산물 수급에 관한 장기전망을 바탕으로 수립하여야 한다(산림기본법 제10-11조). 제6차 산림기본계획 수립 시 반영되었던 산림자원의 범주는 산림면적, 임상면적, 임목축적으로(Korea Forest Service, 2018), 특히 임상면적의 장기전망 결과는 임목축적과 탄소축적에 직접적으로 영향을 주고 목재공급량에는 간접적으로 영향을 미친다. 또한 임상별 면적 정보는 산림 물수지 평가에도 영향을 준다.

임상의 변화는 자연적, 인위적인 여러 가지 요인에 의해서 영향을 받기 때문에, 시간이 경과하면서 기존의 임상이 어떻게 변화할 것인가를 예측하고자 할 때 입지조건과 식생천이, 기후변화, 외부 교란 등 다양한 요인을 함께 고려해야 한다(Lee et al., 1999). 즉, 자연적인 수목의 생장과 경쟁 과정에 의해 해당 입지조건에서 타 임상에게 밀려날 가능성이 있는지, 해당 임상이 유지되기 위한 기본적인 기후조건이 변화할 가능성이 있는지, 외부의 인위적 교란에 의해 임분의 훼손 및 변화 가능성이 있는지 등을 살펴보아야 한다.

이와 관련하여 국내 · 외에서 수종 및 임상분포 변화 예측과 관련된 다양한 연구들이 진행되어 왔다. 주로, 현재의 수종과 임상 분포 특성을 파악하여 모형을 구축하고, 이러한 분포 특징이 미래에도 그대로 유지된다는 가정 하에 기후변화에 따른 미래 분포를 예측하는 접근이 이루어지고 있다. 미국 산림청에서는 FIA(Forest Inventory and Analysis) 표본점에서 조사된 134개 산림수종 자료를 이용하여 수종별 분포적지 평가 및 미래 예측 연구를 수행하고, 그 결과를 국가산림계획 수립에 활용하고 있다(Iverson et al., 2008; USDA Forest Service, 2012). 산림수종별 서식지를 예측하기 위해 기후(기온, 강수량), 지형(고도), 토양(토양형, 물리 · 화학성), 토지이용(산림/비산림비율, 파편화지수 등) 변수 기반의 Random forest 모형을 적용하여 수종별 · 임상별 변화 전망 지도를 제작하여 활용했다. 국내에서는 Choi et al.(2011), Jeoung et el.(2013)과 Yoo et al. (2020)가 기후정보(온량지수, 최저온도지수, 유효강우지수)를 이용하여 수종별 분포범위를 평가하는 HyTAG 모형을 이용하여 기후변화에 따른 임상 및 수종변화를 예측했으며, Chun et al.(2013, 2014)는 수종 분포 정보와 기후, 입지, 임분 특성 정보를 이용해서 GARP 모형에 기반한 미래 잠재수종분포 적지를 예측했다. Shin et al.(2012)은 현재의 임상분포와 환경변수(지형, 최한월온도, 여름평균 강수량, 토양염기포화율, 토양유기물함량)를 이용하여 다항로짓 회귀분석에 의한 임상분포 예측모형을 구축했다. Lee and Kwon(2017)은 전국 시군단위의 임상 비율에 영향을 주는 요인 분석을 위해 기후(강수량, 최고기온), 지형(해발고도, 산림률)와 함께 사회경제적 요인(인구밀도, 임도밀도, 사유림 비율, 시업지 비율, 벌채재적) 등을 추가적인 변수로 활용하였다.

이러한 기존의 종분포모형 연구들은 현재 수종분포와 기후 · 환경 특성과의 관계가 미래에도 동일하게 유지될 것이라는 가정을 기반에 두고 있기 때문에, 시간 경과에 따른 실제 산림의 동적인 변화 메커니즘을 반영하지 못한다는 근본적인 한계가 있다(Vose et al., 2012). 즉 산림의 변화를 설명하는데 있어서 산림생태계는 외부 교란에 대한 적응력이 있고, 기후변화 속도와 수종분포 이동 속도도 동일하지 않으며, 시간 경과에 따라 수목의 생장, 경쟁에 따른 숲의 생태적인 변화 등도 중요하게 고려되어야 한다. 기존 모형연구가 기후변화에 따른 산림의 변화 방향성을 보여주는데 유용하게 활용될 수 있지만, 현실성 있는 산림 변화 전망치 추정을 위해서는 과거 산림이 실제로 어떻게 변화해 왔는지를 파악하고 변화를 유발하는 특성을 도출해 내는 것이 선행되어야 한다.

임상의 실제 변화 특성에 기반한 연구도 진행되고 있으며, 주로 위성영상이나 항공사진 등 원격탐사 자료를 이용한 임상변화 자료에 기반한 연구가 수행되고 왔다. Kim and Lee(2013)는 전국의 구상나무 분포변화, Kim et al.(2017)는 최근 10년간 한라산 구상나무 분포 현황을 조사하였다. Kim et al.(2019)는 과거 20년간 전국 아고산 지역의 침엽수림 변화를 조사하였다. Kim and Jang(2012), Jang and Lee(2013)는 충남, 강원지역을 대상으로 다중시기 위성영상을 이용하여 과거 20년 이상 동안의 임상 변화 지도를 구축하여 변화 추세를 확인하고, 임상변화 지역의 기후 및 지형 특성을 분석했다. 그 결과 전체적으로 침엽수림이 감소하고 활엽수림과 혼효림이 증가하였고 이는 기온상승과 관련성이 높은 것으로 평가되었다. 공간적으로는 침엽수림 가장자리에서 혼효림이 증가하고 활엽수림이 혼효림으로 전환되어 전반적으로 활엽수림의 분포범위가 확장되는 현상도 확인되었다. 산림변화에 대한 실증적 연구는 아직 변화정보 취득을 위한 조사와 기본적인 해석 단계에 와 있으며 이를 기반으로 한 모형개발 연구는 미비하다. 실제 산림변화 특성에 기반한 연구들은 시작단계에 있으며, 특정 지역의 임상변화 정보와 제한된 변수만을 이용함에 따라 다양한 요인들이 복합적으로 연결되어 있는 임상 변화 메커니즘을 설명하는 데는 아직 한계가 있다.

또 다른 접근법으로 국외에서는 수종별 생장과 경쟁을 고려할 수 있는 LANDIS-Ⅱ 모형 등을 이용하여 기후변화에 따른 수종별 생장변화를 시뮬레이션하여 이를 기반으로 분포변화를 예측하는 모형 연구도 수행되고 있다(Duveneck et al., 2017). 수목의 생육, 생장, 경쟁, 확산, 교란 등을 역동적인 생태적 과정으로 모사하여 수종의 변화를 예측할 수 있는 모형을 동적식생모형(Dynamic vegetation model)으로 통칭할 수 있다(Snell et al., 2014). 이러한 모형을 적용하기 위해서는 매우 과학적이고 엄밀한 국내 기반자료들이 구축되고 검증된 이후에 적용될 수 있을 것이다.

따라서, 본 연구는 직접 관찰된 전국적인 임상변화 현장 조사 자료를 기반으로 임상변화를 발생시키는 다양한 요인을 파악하고 변화 특성을 해석하며, 이를 기반으로 임상 변화에 대한 실증적인 모형을 개발하여 미래 임상변화 예측치를 제시하는 것을 목표로 하였다. 이를 위해, 시계열 국가산림자원조사 자료를 이용하여 전국 산림의 10년 동안의 임상변화 정보를 확보하였고, 임상변화에 영향을 미칠 수 있는 다양한 변수(기후, 지형, 임분, 교란 등)를 구축하고, 이를 기반으로 임상의 변화와 유지에 영향을 주는 변수를 평가하였다. 또한 유효변수들을 이용하여 10년 경과에 따른 임상변화모형을 구축하고 기후변화시나리오 자료를 이용하여 미래 임상분포를 예측했다.

재료 및 방법

1. 분석자료

임상변화의 정보는 국가산림자원조사(National Forest Inventory, NFI) 시계열 조사 자료를 통해 취득하였다. 국가산림자원조사에서는 제5차 조사부터 전국 산림을 대상으로 4 km × 4 km의 일정한 간격으로 고정표본점을 배치하고 해당 표본점에 대해 5년 주기로 반복적인 현장조사를 수행하고 있다. 따라서 국가산림자원조사 자료는 산림의 변화를 가장 명확하게 파악할 수 있는 유일한 국가 수준의 현장조사 자료라고 할 수 있다. 본 연구에서는 10년 동안의 임상의 변화 상태를 확인하기 위해 제5차 NFI 자료(2006~2010)와 제7차 NFI(2016~2020) 자료를 비교 대상 자료로 활용했다(Table 1).

Table 1. Status of national forest inventory data used for analysis.

HOMHBJ_2022_v111n4_461_t0001.png 이미지

국가산림자원조사 표본점 자료를 이용한 변화 분석을 위해서는 두 시기의 표본점 정보가 실제로 변화 비교가 가능한 표본점인지를 우선적으로 평가해야 한다. 조사표본점은 여러 가지 요인으로 인해 표식이 소실되기도 하고 표본점 지역의 훼손으로 인해 부득이하게 인근 지역에 신규표본점이 설치되기도 한다. 이로 인해 표본점 번호는 동일하다 하더라도 과거와 동일지역에 대한 조사가 수행되었다고 판단할 수 없을 경우, 해당 표본점 조사정보는 전국 산림통계 산출에는 이용이 가능하지만 변화 모니터링 정보로서는 활용이 어려워진다. 이러한 표본점의 경우 모니터링 불가 판정이 내려지게 된다. 본 연구에서는 동일 임분의 변화에 대한 정보가 기본 전제가 되기 때문에 동일 지점에서 조사되어 모니터링 가능으로 평가된 표본점에 한정하여 분석에 활용했다. 또한 본 연구에서는 자연적인 숲의 변화 특성을 파악하는 것이 목적이기 때문에, 시업 이력이 있는 표본점은 분석에 포함하지 않았다.

NFI6과 NFI7 모두에서 모니터링이 가능하다고 평가되어 최종적으로 NFI5와 NFI7 자료간의 비교가 가능한 표본점은 Table 1에서 제시한 전체 표본점 중에서 12,833개소가 해당되었다. 또한 자연적인 임상변화 특성을 평가하기 위해, 조사수행 기간 내에 산림시업이 수행된 것으로 확인된 표본점을 제외한 총 9,714개 표본점 자료를 이용하여 분석을 수행했다.

2. 임상변화 요인 분석 및 모형개발 변수 구축

본 연구에서 임상(forest type)은 침엽수림(CON), 활엽수림(DEC), 혼효림(MIX)의 구분을 말한다. 임상구분 기준은 수치임상도의 임상구분 기준에 따라 주림목(우세목, 준우세목)을 대상으로 한 흉고단면적합 비율을 기준으로 하였을 때 침엽수의 비율이 75% 이상인 지역은 침엽수림, 25~75%인 지역은 혼효림, 25% 미만인 지역은 활엽수림으로 정의하였고(Cho et al., 2012), 각 표본점 내 개체목별 수종과 흉고직경(DBH) 정보를 이용하여 임상을 구분했다.

10년 간의 임상변화의 방향은 침엽수림 →혼효림, 활엽수림 → 혼효림, 혼효림 → 침엽수림, 혼효림 → 활엽수림으로 규정했다. 침엽수림 → 활엽수림, 활엽수림 → 침엽수림의 급격한 변화를 보이는 표본점의 경우 전체 표본점 중 각각 0.9%, 0.1%로 비율이 매우 낮으므로, 본 연구에서는 별도로 고려하지 않았다.

임상변화 유형별로 별도의 모형을 구축하였으며, 침엽수림과 활엽수림의 경우 모형의 종속변수는 임상의 유지를 0, 임상의 변화를 1로 설정했다. 혼효림의 경우 유지는 0, 침엽수로의 변화는 1, 활엽수로의 변화는 2로 지정했다. 또한 임상변화를 설명할 변수로는 기상 · 기후변수, 지형변수, 임분변수, 교란변수를 구축했다(Table 2). 임상변화 모형 구축을 위한 자료는 기본적으로 표본점 조사정보를 활용하되, 현장조사를 통해 산출되지 않은 자료들은 별도 공간자료로 제작하여 표본점 위치에 해당하는 자료를 추출하여 포함하였다.

Table 2. Variables used to predict forest type change.

HOMHBJ_2022_v111n4_461_t0002.png 이미지

기상 · 기후 부문에서는 임상 변화가 발생한 대상기간(2006~2020년) 내 연평균기온과 연누적강수량 평균 자료를 구축해 활용했다. 자료는 기상청에서 제공하고 있는 1 km 해상도의 과거 격자 기상자료(MK-PRISM ver1.2)를 이용하였고, 해당 자료는 기상청 기후정보포털(http://www.climate.go.kr)에서 다운로드했다. 그러나 과거 격자 기상자료의 경우 2017년까지만 제공하고 있기 때문에, 2018~2020년까지의 격자기상자료는 1 km 남한상세 기후 변화시나리오(RCP 8.5) 자료를 이용하였다. 지형의 경우 고도, 경사, 사면, 지형학적수분지수 자료를 구축했다. 고도, 사면향, 경사 자료는 표본점 조사정보를 그대로 이용하였고, 사면향의 경우 사면에 따른 건조 영향을 정량적으로 분석하기 위해, 산지의 남서사면이 가장 건조도가 높고 북동사면이 낮은 점을 고려하여 남서사면(SW)을 90, 북동사면(NE)을 -90으로 환산하여 전체 사면향 값을 변환하여 활용했다. 즉, 사면향의 값이 90에 가까우면 건조, -90에 가까우면 습윤한 상태를 말한다. 지형학적수분지수(TWI, Topographical Wetness Index)는 사면형태에 따른 수분조건을 반영한 지수(Beven and Kirkby, 1979)로서, 전국 30 m 수치표고모델을 이용하여 산출했다. 임분변수는 표본점 임분조사 자료에서 제시하고 있는 영급과 수관밀도 자료를 이용하였다. 또한 임분의 유지와 변화에 영향을 줄 수 있는 교란요인으로는 가장자리 효과를 고려하였다. 가장자리는 외부와 접하고 있어 교란가능성이 높은 지역으로, 기존 연구에서도 임상의 변화가 해당 임상의 가장자리에서부터 변화하는 현상이 확인된 바 있다(Jang and Lee, 2013). 표본점 주변지역 산림 내 침엽수림 비율, 산림 내 활엽수림 비율, 비산림 비율 및 도로로부터의 거리 정보가 외부 교란 요인으로 검토되었다. 주변 임상비율 산정을 위해 2006~2010년 기간 내의 수치임상도를 이용하여야 했으나 과거 임상도는 1:25,000 축적의 임상도로 현재의 1:5,000 임상도와의 공간해상도 불연속이 있어 활용이 불가했으며, 또한 본 연구에서는 기존 산림통계와의 연계성을 위해 2015년(항공사진 촬영은 그 이전 시기) 제작된 혼효림 보정 임상도(임업진흥원 제공)를 이용해야 하는 제한이 있어서 본 정보 또한 위 임상도를 이용하여 산출되었다. 본 변수는 6차 산림자원조사(2011~2015) 기간에 해당하는 정보를 담고 있으므로 변수가 반영하는 기간은 임상변화 초기 시기인 5차 산림자원조사 기간과의 약간의 불일치가 발생하는 한계를 담고 있다. 도로로부터의 거리 정보는 NFI5 표본점 임분조사 자료를 활용했다(Table 2).

3. 임상변화 영향요인 분석 및 예측모형 선정

임상의 변화를 발생시키는 요인을 평가하기 위해 유형별 경향성을 평가하고 임상변화 여부에 따른 통계적 차이를 분석하였다. 해당 변수들은 독립변수 t-test 및 분산분석을 통해 임상변화 여부에 대한 통계적 유의성을 검증하고, 임상변화 예측모형에 활용할 수 있는 변수를 1차적으로 선정했다(유의확률이 0.1보다 작은 변수). 변수 간 상관성이 높은 변수는 상대적으로 유의성이 높은 변수를 선택했다.

임상변화를 예측하기 위한 후보모형으로 로지스틱 회귀모형, Random Forest(RF), Support Vector Machine(SVM)을 이용했다. 로지스틱 회귀모형은 주어진 설명변수 조건하에서 특정한 선택이 이루어질 확률의 최대값이 1, 최소값이 0인 S자형 곡선 형태를 가정했을 때 로지스틱 함수를 활용하는 모형으로(Cox, 1958), 함수를 로짓으로 변환시키면 선형으로 나타나므로 일반적인 선형회귀모형과 유사하게 활용될 수 있다. 혼효림의 경우에는 다항 로지스틱 회귀모형을 적용했다. 변화와 미변화를 설명하기 위한 가장 간단한 선형모형이라 할 수 있다. RF는 의사결정트리가 여러 개 모여져서 앙상블 기법으로 결과를 산출하는 모형이다. 의사결정트리는 특정 항목에 대한 의사결정 규칙을 나무 형태로 분류해 나가는 분석기법으로, RF는 단일 의사결정트리 모델의 과적합을 방지하기 위해 고안된 방법이다(Breiman, 2001). 샘플링을 통해 모델 구성을 위한 여러 가지 자료 세트를 만들고, 이에 따라 여러 종류의 의사결정트리를 만든 이후에 투표를 통해 결과를 산정하는 방식을 취한다. 또한 변수 중요도 평가를 통해 높은 가중치가 부여된 변수를 평가할 수 있는 장점이 있다. SVM은 데이터들을 가장 큰 거리로 떨어지도록 분리해내는 선, 면, 또는 비선형적 공간구분의 체계를 찾아내는 알고리즘을 말한다(Cortes and Vapnik, 1995). 본 연구에서는 선형 SVM 모형을 적용했다.

최적모형 평가를 위해 임상변화에 대한 동일한 자료 세트를 구축하고 세 가지 모형을 적용했다. 임상이 변화하지 않은 샘플개수는 변화한 샘플 개수보다 훨씬 많은 비중을 차지하기 때문에, 임상 변화/미변화 모형도출 시 임상 미변화 정보 특성의 영향력이 상대적으로 커지는 문제가 발생했다. 이러한 문제를 방지하기 위해, 샘플개수가 적은 경우의 샘플 개수를 늘려 비교대상이 되는 경우의 샘플개수와 동일하게 맞춰주는 Upsampling 기법을 적용했다. 또한 세 모형 간의 비교를 위해 K-fold Cross Validation 방법을 적용하고, 분류정확도를 평가하는 전체 정확도(Overall Accuracy), 카파계수(Kappa Coeff.), 예측된 변화 비율과 실제 변화 비율의 차이값을 이용하여 모형 성능을 평가했다. 통계분석과 모형구동은 R 프로그램으로 수행되었다.

4. 임상변화의 미래 전망

임상변화를 설명하는데 유효한 변수를 이용하여 모형을 구동하고, 가장 성능이 좋은 모형을 선정하였다. 그리고 이 모형을 이용하여 미래 임상변화를 예측했다. 임상변화에 영향을 주는 요인 중 지형 조건은 미래에도 변화하지 않는다고 가정하였으며, 기상기후, 임상, 임분특성은 변화한다고 가정했다.

임상 변화가 발생할 미래 대상기간의 연평균기온과 연누적강수량 평균 자료는 기상청에서 제공하고 있는 1 km 남한상세 기후변화시나리오(RCP 8.5)를 이용하여 산출했다. 초기 모형 구동은 2015년 혼효림 보정 임상도(임업진흥원 제공)를 입력자료로 이용하였으며, 이를 통해 예측된 10년 후 임상분포 자료는 다음차기 모형구동을 위한 임상 분포 입력자료 및 주변지역 침엽수림 및 활엽수림 면적비율 산출 기초자료로 이용되었다.

임상이 변화했을 때 주림목의 영급과 수관밀도가 변화할 수 있다. 그러나 변화 임분에 대한 영급과 수관밀도 변화에 대한 일관적인 특성이 확인되지 않았으며, 이에 따라 본 연구에서는 10년 단위 미래 전망 시 임령이 1씩 증가(10년)하고, 수관밀도는 변화하지 않는다고 가정했다.

결과 및 고찰

1. 임상 변화의 경향

NFI5와 NFI7 간의 변화분석이 가능한 표본점은 총 12,833개였으며, 이 중 인위적 시업활동이 확인된 표본점을 제외한 9,714개 표본점 자료를 활용하였다. 그 결과, 10년 후 임상이 유지되는 비율이 침엽수림의 경우는 84%, 활엽수림은 97%였으며, 임상의 변화는 대부분 혼효림으로 이동하는 경우였다. 혼효림은 84%가 유지되었으며, 임상이 변화하는 경우에는 침엽수림보다 활엽수림으로 변화되는 비율이 4배 이상 높았다(Table 3). 이를 요약하면, 침엽수림은 활엽수림보다 더 많이 혼효림으로 변화하고 혼효림은 침엽수림보다 활엽수림으로 더 많이 변화하여, 전반적으로 침엽수림은 감소하고 활엽수림은 증가하는 방향으로 임상이 변화하였다. 침엽수림 감소 비율은 서울·인천·경기 지역에서 가장 높았고 대구·경북에서 가장 낮았다(Table 4).

Table 3. Forest type change of NFI permanent sample points for 10 years.

HOMHBJ_2022_v111n4_461_t0003.png 이미지

Table 4. Forest type change of coniferous forest for 10 years based on NFI data.

HOMHBJ_2022_v111n4_461_t0004.png 이미지

2. 임상변화 영향 요인과 특징

위와 같은 임상변화를 발생시킨 요인을 파악하기 위해 임상의 변화와 12개 변수들(기상‧기후, 지형, 임분, 교란)과의 관계를 파악하였다. 기상 · 기후요인의 경우, 활엽수림과 혼효림 간의 임상변화 지역은 임상의 유지지역보다 기온이 높은 경향을 보였다. 강수량이 많은 지역은 전체적으로 침엽수림에서 혼효림으로, 혼효림에서 활엽수림으로 변화하는 경향을 보였다(Figure 1).

HOMHBJ_2022_v111n4_461_f0001.png 이미지

Figure 1. Relationship between forest type change and climate factors.

※ C→C: Maintenance as CON, C→M: Change from CON to MIX, D→D: Maintenance as DEC, D→M: Change from DEC to MIX, M→M: Maintenance as MIX, M→C: Change from MIX to CON, M→D: Change from MIX to DEC

활엽수림 변화는 유지지역보다 상대적으로 해발고도가 낮은 곳에서 우세한 것으로 확인되었으나 타 임상에서는 뚜렷한 차이가 없었으며, 경사도의 경우 활엽수림과 혼효림의 변화가 발생하는 지역에서 경사가 다소 낮았으나 뚜렷한 경향성을 보이지는 않았다. 수분 관련 지형요인을 보았을 때, 침엽수림 유지지역은 활엽수림 유지지역보다 사면향 값이 높아 지형조건에 따른 건조도가 높은 지역에 위치하는 특징을 보였으며, 임상변화 역시 건조조건에서는 침엽수림쪽으로, 습윤조건에서는 활엽수림 쪽으로의 변화가 우세했다. 임상의 변화는 유역의 상류쪽(낮은 TWI)보다 하류쪽(높은 TWI)에서 더 높은 경향을 보였다(Figure 2).

HOMHBJ_2022_v111n4_461_f0002.png 이미지

Figure 2. Relationship between forest type change and topographical factors.

임상변화 경향과 임상 상태의 관계를 비교한 결과, 임상이 유지되지 않고 변화하는 지역이 영급이 상대적으로 낮고 수관밀도도 낮은 것으로 나타났다(Figure 3). 즉, 임령이 낮을수록 변화 가능성이 크다는 것을 반영한다.

HOMHBJ_2022_v111n4_461_f0003.png 이미지

Figure 3. Relationship between forest type change and forest stand factors.

교란요인의 경우도 임상변화를 설명하는데 주요 변수로 활용될 수 있다. 산림 내 침엽수림 면적비율이 상대적으로 낮은 지역에서 침엽수림에서 혼효림으로, 혼효림에서 활엽수림으로 변화되는 비율이 높았다. 임상이 변화하는 지역에서 비산림면적 비율이 더 높은 경향을 보였다(Figure 4). 도로에서부터의 거리는 뚜렷한 경향성을 보이지 않았다.

HOMHBJ_2022_v111n4_461_f0004.png 이미지

Figure 4. Relationship between forest type change and disturbance factors.

3. 임상변화 예측모형 설명변수 선정 및 모형 평가

임상변화 예측모형 구축을 위한 후보변수를 선정하기 위해, 우선적으로 침엽수림, 활엽수림, 혼효림 변화/미변화에 영향을 주는 변수(p<0.1)를 1차적으로 선정했다. 전체 변수들 중에서 CON100과 DEC100, TAVG와 ELEVATION의 경우는 변수 간 상관성이 높아(r=-0.75, -0.73) 상대적으로 중요도가 높은 변수를 선별하여 적용하였다(표본점 주변지역 설정의 거리 기준은 100, 200, 300 m 거리를 평가하였으나 본 연구에서는 임상 변화와 관련성이 상대적으로 높았던 100 m 기준으로 설정한 변수를 적용하였다). 결과적으로, 침엽수림의 변화와 관련해서는 교란, 임분, 지형, 기후 관련 변수 6개(CON100, AGE, DENSITY, ASPECT, PRCP, NF100)가 임상변화 구분에 유의한 변수로 추출되었다(Table 5). 특히 침엽수림 면적비율과 임분 특성이 유의성이 높은 것으로 나타났다. 활엽수림 변화와 관련해서는 교란, 기후, 임분, 지형 관련 변수 7개(DEC100, PRCP, DENSITY, TAVG, ASPECT, AGE, NF100)가 임상변화 구분에 유의한 변수로 추출되었다(Table 6). 활엽수림 면적 비율과 기후 특성이 유의성이 높은 것으로 나타났다. 혼효림 변화와 관련해서는 교란, 임분, 기후, 지형관련 변수 5개(CON100, AGE, NF100, TAVG, TWI)가 임상변화 구분에 유의한 변수로 추출되었다(Table 7). DEC100는 CON100과의 상관성이 높아 제외하였다.

Table 5. Variable assessment(t-test) for coniferous forest change model.

HOMHBJ_2022_v111n4_461_t0005.png 이미지

Table 6. Variable assessment(t-test) for broadleaf forest change model.

HOMHBJ_2022_v111n4_461_t0006.png 이미지

Table 7. Variable assessment(anova analysis) for mixed forest change model.

HOMHBJ_2022_v111n4_461_t0007.png 이미지

위에서 1차로 선정된 변수들을 이용하여 LR, RF, SVM 모형을 적용한 결과, K-fold Validation의 전체 정확도, Kappa계수 평균값 모두 SVM이 상대적으로 양호한 것으로 산출되었으며, 예측된 침엽수림 변화 비율과 실제 관측된 변화 비율이 가장 유사하게 산출된 모형 역시 SVM이었다(Table 8). 이에 따라, 최종 모형은 SVM으로 선정하였다.

Table 8. K-fold Validation of forest type change models.

HOMHBJ_2022_v111n4_461_t0008.png 이미지

* Change from MIX to CON, ** Change from MIX to DEC

SVM 침엽수림 모형의 경우 주변의 침엽수림 면적 비율과 임분 특성이 상위 중요변수로 평가되었고, SVM 활엽 수림 모형에서는 주변의 활엽수림 면적 비율과 기온 및 강수량 등 기후정보가 상위 중요변수로 평가되었다.

4. 임상분포 변화 미래 전망

임상변화 예측을 가장 효과적으로 수행하는 것으로 평가된 SVM 모형을 적용하고, 모형구축과 검증을 위한 샘플을 구분했던 것을 모두 포함해 전체 표본점 자료를 이용해 최종 모형을 도출했다. 그리고 공간입력자료와 모형을 이용해서 10년단위로 미래 임상변화를 예측했다. 그 결과, 2015년에서 2055년까지 40년 기간 동안 침엽수림은 38.1%에서 28.5%로 감소, 활엽수림은 34.2%에서 38.8%로 증가, 혼효림은 27.7%에서 32.7%로 증가할 것으로 예측되었다(Table 9, Figure 5). 침엽수림은 첫 10년에 급격히 감소하고 이후에는 점진적으로 감소하는 경향을 보였고, 활엽수림은 2035년까지 증가하다가 그 이후에는 다소 감소하는 것으로 나타났다. 혼효림은 전체적으로 약간 증가하는 경향을 보였다.

HOMHBJ_2022_v111n4_461_f0005.png 이미지

Figure 5. Change of Forest type distribution in future based on RCP 8.5 scenarios.

Table 9. Results of forest type change prediction using SVM model.

HOMHBJ_2022_v111n4_461_t0009.png 이미지

이러한 변화 추세는 지역별로 차이가 많이 발생했다. 전국적으로 모든 지역에서 침엽수림 면적이 감소했고 그 중 강원도와 광주 · 전남 지역의 감소면적이 가장 컸다. 활엽수림은 광주· 전남지역이 가장 많이 증가했고 강원도와 경상남· 북도 등 감소한 지역도 있었다. 혼효림은 전 지역에서 증가했고 특히 강원 지역의 증가가 가장 컸다. 비율로 보면, 침엽수림은 제주에서 가장 많이 감소하였고, 대구‧경북이 가장 적은 비율로 감소했다. 활엽수림은 광주‧전남에서 가장 많이 증가했고, 대구‧경북은 활엽수림이 약간 감소했다(Table 10). 이러한 점을 종합해볼 때, 가장 두드러진 변화가 발행한 지역은 강원과 전남이었다. 강원도는 침엽수림 감소와 혼효림 증가, 전남지역은 침엽수림 감소와 활엽수림 증가가 우세하게 나타났다. 비율적으로는 침엽수림 면적 자체가 적고 파편화되어 있는 제주도나 서울·인천·경기 지역의 감소 비율이 높고 침엽수림 면적이 넓고 군집화 되어 있는 경북 지역의 감소 비율이 낮았다.

Table 10. Forest type change from 2015 to 2055.

HOMHBJ_2022_v111n4_461_t0010.png 이미지

결론

지난 10년 동안의 임상 변화 추세와 특성을 파악한 결과 몇 가지 현상을 확인할 수 있었다. 첫째, 우리나라 산림은 침엽수림이 감소하고 혼효림과 활엽수림이 증가하는 방향으로 변화하고 있다. 침엽수림 감소비율은 서울·인천·경기, 대전·세종·충남, 광주·전남 순으로 높았으며 지역별 침엽수림 면적을 고려하면, 광주·전남의 침엽수림 감소 수준이 가장 컸을 것으로 추정된다. 둘째, 침엽수림에서 혼효림, 혼효림에서 활엽수림 방향으로 변화하는 지역은 주로 지형적으로 습윤하고 강수량이 많아서 수분 관련 생육환경이 양호하며, 주변에 활엽수림이 많은 지역이었다. 셋째, 변화의 방향과 관계없이 임상변화가 상대적으로 많이 발생하는 지역은 기온이 높은 지역, 임분의 임령과 밀도가 낮은 지역, 주변 지역에 비산림이 많은 지역 등 교란 가능성이 높은 지역이었다.

수분 관련 생육환경이 양호한 지역에서 침엽수림이 활엽수림으로 변화하는 경향이 많은 것으로 확인되었으며, 이는 식생천이 현상과 관련이 깊다. 100년 전 남한지역 입목지 면적 중 70% 이상이 소나무였으며 상당 부분이 치수림 상태였다(Bae and Kim, 2019). 1970년대에는 황폐한 산림의 복구를 위해 척박하고 건조한 지역에도 잘 자라는 소나무의 육성이 활발히 이루어졌다. 그리고 이후에는 천이 초기 수종인 소나무를 중심으로 한 산림생태계가 잘 유지되면서 점차 산림 내 토양과 미기후 등의 생육환경이 양호해지고, 이후 이를 기반으로 하여 중성수인 참나무류가 생장이 점차 좋아지면서 결국에는 양수인 소나무와의 경쟁에서 이겨 숲의 상층부를 지배하는 과정이 발생하고 있다(Lee et al., 1999). 그리고 종간 경쟁은 타 수종과 가까이 접하고 있는 임분의 가장자리에서 활발히 이루어진다. 이러한 자연적인 천이과정은 토양비옥도와 수분조건 등 생육환경의 개선과 관련성이 높으나 추후 기후변화의 방향에 따라 숲이 더 건조해지거나 습윤해질 수 있기 때문에, 자연적인 숲의 변화와 인위적인 기후변화를 동시에 고려해서 미래를 전망할 필요가 있다. 또한 임상의 변화 가능성이 높다는 것은 임상이 안정화되어 있지 않다는 의미이며, 숲의 성숙도가 낮아 아직 변화 가능성이 많은 지역, 외부 교란 요인이 많은 지역 등이 이에 해당한다.

이러한 임상의 변화 특성을 기반으로 NFI5와 NFI7의 10년 사이의 임상 변화에 대한 경험적 모형을 구축하였다. 그리고 과거 10년 사이의 임상변화 특성이 미래에 동일하게 지속적으로 이어졌을 경우에 미래 임상은 어떻게 변화할 가능성이 있는지를 평가하였다. 그 결과, 산림 내 생육 환경 개선으로 침엽수림은 40년 기간 동안 지속적으로 감소했으나 임령의 증가로 감소속도는 다소 둔화되었다. 활엽수림은 2035년까지는 증가하였으나, 이후 혼효림으로의 비중 증가로 약간의 감소추세를 보였다. 혼효림은 점진적으로 증가하는 것으로 나타났으며, 이는 산림의 변화압력이 지속적으로 커지고 있는 것이 반영된 것이라 볼 수 있다. 지역별로는 침엽수림의 감소면적과 활엽수림 증가 면적이 모두 큰 광주전남 지역의 임상변화 압력이 가장 클 것으로 예상되었다. 침엽수림 면적 자체가 적고 파편화되어 있는 제주도나 서울· 인천· 경기 지역의 침엽수림 감소 비율이 높고 침엽수림 면적이 넓고 군집화 되어 있는 경북 지역의 감소 비율이 낮았다. 이를 종합하면, 침엽수림 비중이 낮은 지역에서의 변화 압력, 남부지방 기후특성에따른 변화압력이 크게 작용하고 있다는 것을 알 수 있다.

본 연구의 임상변화 분석 및 미래 전망에서 고려하지 못한 부분들이 있다. 첫째, 기온이 높은 남부지역에서 임상변화가 많은 경향을 보였는데 실제로 높은 기온이 수종 간의 경쟁요인을 유발시켜 임상이 변화하는 데 직접적인 영향으로 작용했는지는 분명하게 확인되지 못했다. 또한 수종별로 적정생육기온이 다르고 활엽수와 침엽수 내에서도 수종간 특성 차이가 크므로, 임상이라는 큰 분류군으로 구분된 체계 하에서는 기온의 영향을 명확히 해석하기가 어려웠다. 즉, 침엽수림 내에는 잣나무, 낙엽송 등 중부북부에 주로 서식하는 수종도 있는 반면 편백, 삼나무 등 기후적으로 남부지방에 적합한 수종도 있다. 활엽수림 내에도 낙엽 활엽수림과 남부지방 중심으로 서식하는 상록 활엽수림이 함께 포함되어 있다. 산림의 변화 요인과 메커니즘에 대해서 보다 세밀하게 접근하기 위해서는 수종 단위의 분석이 수행되는 것이 필요하다.

둘째, 숲의 변화는 식생천이라는 큰 변화의 흐름과 병해충, 산불, 이상기후에 따른 수목고사 등 산림교란이라는 이벤트가 함께 영향을 주고 있다. 시간과 공간 규모가 다른 이 두 가지 측면의 변화 요인이 시너지 관계인지 혹은 상충 관계에 있는지는 아직 명백하게 밝혀지지 않았다. 그러나 기후변화와 여러 가지 재해요인으로 인한 상록침엽수의 피해가 숲의 천이와 맞물려 변화의 속도를 가속화시킬 것이라는 예측도 가능하다. 향후 숲에서 발생한 교란이 임상의 변화에 어떤 비중으로 역할을 해왔는지에 대한 정량적 분석을 수행하고 그 결과를 미래 예측에 반영하는 것이 필요하다.

셋째, 본 연구에는 국가산림자원조사(NFI) 표본점 자료 기반으로 분석되었다. NFI는 전국을 격자화하여 촘촘하게 조사되는 가장 신뢰도 높은 전국조사 자료임에도 불구하고, 전체 산림지역 중 대상지를 추출하여 조사하는 방식이기 때문에 관찰하고자 하는 대상이 누락될 수 있는 것은 근본적 한계를 포함하고 있다. 예를 들면, 산불, 병해충 발생 등 교란을 통해 소실되거나 변화하는 지역에 표본점이 배치되어 있지 않으면 이러한 변화는 탐지될 수 없는 것이다. 따라서, 전반적인 변화특성을 세밀하게 파악하기 위한 국가산림자원조사 지점 자료와 항공사진이나 위성 자료 등을 통한 산림교란 모니터링 공간정보가 함께 활용될 필요가 있다. 또한 산림시업 활동에 의한 임상변화 특성 분석, 침엽수림에서 활엽수림으로 큰 변화를 겪은 일부 소수 표본점에 대한 세밀한 요인 분석 연구도 필요하다.

본 연구의 산림의 변화 특성에 기반한 임상분포 미래전망 결과는 공간정보에 기반한 지역규모의 세밀한 산림관리를 위한 구체적 실행 전략을 세우는 데 활용될 수 있다. 단, 임상변화 경험적 모형을 기반으로 한 본 미래 전망 결과는 과거의 변화 특성이 미래에도 그대로 유지된다는 가정하에 분석되는 것이므로 불확실성이 많은 미래를 예측하는 데 있어서는 분명한 한계가 있음을 인지해야 한다. 향후 수종별 변화 특성의 정밀한 해석, 산림 교란을 고려한 미래 전망, 산림생태계 변화 프로세스 모형 연구 등 지속적 연구 수행이 필요하다.

References

  1. Bae, J.S. and Kim, E.S. 2019. Understanding forest status of the Korean peninsula in 1910: A focus on digitization of Joseonimyabunpodo(The Korean peninsula forest distribution map). Journal of Korean Society of Forest Science 108(3): 418-428.
  2. Beven, K.J. and Kirkby, M.J. 1979. A physically based, variable contributing area model of basin hydrology. Hydrological Science Bulletin 24(1): 43-69. https://doi.org/10.1080/02626667909491834
  3. Breiman, L., Friedman, J.H., Olshen, R.A. and Stone, C.J. 1984. Classification and regression trees. Monterey, CA: Wadsworth & Brooks/Cole Advanced Books & Software.
  4. Cho, H.K., Yoo, B.O., Kim, S.H., Ryu, J.H., Kim, J.C., Seo, S.H., Kim, J.S., Lee. J.B., and Kim, H.H. 2012. Production and update of large scale forest type map using the digital airphoto. Korea Forest Research Institute. Seoul. Korea. pp. 135.
  5. Choi, S., Lee, W.K., Kwak, D.A., Lee. S., Son. Y., Lim, J.H. and Saborowski, J. 2011. Predicting forest cover changes in future climate using hydrological and thermal indices in South Korea. Climate Research 40(3): 229-245.
  6. Chun, J.H. and Lee, C.B. 2013. Assessing the effects of climate change on the geographic distribution of Pinus densiflora in Korea using Ecological Niche Model. Korean Journal of Agricultural and Forest Meteorology 15(4): 219-233. https://doi.org/10.5532/KJAFM.2013.15.4.219
  7. Chun, J.H., Lee, C.B and Yoo, S.M. 2015. Shifts of Geographic Distribution of Pinus koraiensis based on climate change scenarios and GARP model. Korean Journal of Agricultural and Forest Meteorology 17(4): 348-357. https://doi.org/10.5532/KJAFM.2015.17.4.348
  8. Chun, J.H., Shin, M.Y., Kwon, T.S., Lim, J.H., Lee, Y.G., Park, G.E., Kim, T.W. and Sung, J.H. 2014. Predicting the changes of productive areas for major tree species under climate change in Korea. Korea Forest Research Institute. pp. 240.
  9. Cortes, C. and Vapnik, V. 1995. Support-vector networks. Machine Learning 20(3): 273.
  10. Cox, D.R. 1958. The regression analysis of binary sequences. Journal of the Royal Statistical Society. Series B (Methodological) 20(2): 215-242.
  11. Duveneck, M.J., Thompson, J.R., Gustafson, E.J., Liang, Y. and de Bruijn, A.M.G. 2017. Recovery dynamics and climate change effects to future New England forests. Landscape Ecology 32(7): 1385-1397. https://doi.org/10.1007/s10980-016-0415-5
  12. Iverson, L.R., Prasad, A.M., Matthews, S.N. and Peters, M. 2008. Estimating potential habitat for 134 eastern US tree species under six climate scenarios. Forest Ecology and Management 254(3): 390-406. https://doi.org/10.1016/j.foreco.2007.07.023
  13. Jang, D.H. and Lee, S. 2013. Analysis of forest type changes for vulnerability assessment of climate changes in the mounyainous resgions of Gangwon-do, Korea. Journal of climate research 8(2): 169-183. https://doi.org/10.14383/cri.2013.8.2.169
  14. Jeong, H.Y., Lee, W.-K., Nam, K. and Kim, M. 2013. Estimating the change of potential forest distribution and carton stock by climate changes - Focused on forest in Yongin-City. Journal of Climate Change Research 4(2): 177-188.
  15. Kim, E.S., Kim, C.M., Kim, K.M., Ryu, J.H., Yim, J.S. and Kim, J.C. 2015. Change of Korea's national forest inventory system (1971~2010). Korea Forest Research Institute. pp. 150.
  16. Kim, J.K., Koh, J.-G., Yim, H.-T. and Kim, D.-S. 2017. Changes of spatial distribution of Korean fir forest in Mt. Hallasan for the past 10 years (2006, 2015). Korean Journal of Environment and Ecology 31(6): 549-556. https://doi.org/10.13047/KJEE.2017.31.6.549
  17. Kim, J.-S. and Jang, D.H. 2012. Analysing of forest types in Chungnam coastal area using multi-temporal satellite imagery and ASTER DEM data. Journal of the Korean Geomorphological Association 19(1): 69-81.
  18. Korea Forest Service, 2018. The Sixth National Forest Plan (2018~2037). Korea Forest Service. pp. 151.
  19. Korea Forest Service. 2021. 2021 Statistical yearbook of forestry. Korea Forest Service. pp. 460.
  20. Lee, H.L. and Kwon, O.S. 2017. Analyzing the impac ts of climate change on forest composition in Korea. Environmental and Resource Economics Review 26(2): 229-255. https://doi.org/10.15266/KEREA.2017.26.2.229
  21. Lee, K.J., Han, S.S., Kim, J.H. and Kim, E.S. 1999. Forest ecology. Hyangmunsa. pp. 395.
  22. Snell, R.S., et al. 2014. Using dynamic vegetation models to simulate plant range shifts. Ecography 37(12): 1184-1197. https://doi.org/10.1111/ecog.00580
  23. USDA Forest Service. 2012. Future of America's forest and rangelands: Forest service 2010 resources planning act assessment. General Technical Report. WO-87. Washington, DC. pp. 198.
  24. Vose, J.M., Peterson, D.L., Patel-Weynand, T. 2012. Effects of climatic variability and change on forest ecosystems: A comprehensive science synthesis for the U.S. forest sector. General Technical Report. PNW-GTR-870. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station. pp. 265.
  25. Yoo, S., Lim, C.-H., Kim, M., Song, C., Kim, S.J. and Lee, W.-K. 2020. Potential distribution of endangered coniferous tree species under climate change. Journal of Climate Change Research 11(4): 215-226. https://doi.org/10.15531/KSCCR.2020.11.4.215
  26. Climate Information Portal of Korea Meteorological Administration. http://www.climate.go.kr. (2021. 7. 8)