A transiting temperate-subtropical mixed forest: carbon cycle projection and uncertainty

Terrestrial ecosystems respond to climate change in various ways, making it crucial to improve our understanding of these dynamics and uncertainty in projections. Here, we investigate how the species composition in a temperate-subtropical mixed forest on Jeju Island, South Korea, would change by 2099 and analysed the resultant effects on phenological timings and carbon flux using an individual cohort-based model—the ecosystem demography biosphere model version 2. We use the analyses of variance to decompose the contribution of model parameters (four sets) and climate inputs (four global climate models under four representative concentration pathway (RCP) scenarios) to the total uncertainty in the leaf area index (LAI) and net ecosystem productivity (NEP) projections. We find that with increases in temperature, photosynthetically active radiation, and vapour pressure deficit, the dominance of subtropical species will gradually increase by approximately 11%, from 30.2% in 2013 to 41.1% by the end of this century, yet there was a large variation in the projections depending on the model parameter and climate inputs. We also show the increases in the LAI and length of growing season by the end of this century, resulting in an increased NEP at the rate of up to 62.7 gC m−2 yr−1 per decade under the RCP8.5. The uncertainty in the LAI projection was largely due to the model parameter (and its interaction with climate inputs); however, the uncertainty contribution of climate models is as large as the emission scenario in the NEP projection. This study highlights the importance of identifying uncertainty sources for a robust projection of terrestrial ecosystem and carbon cycle.


Introduction
Terrestrial ecosystems have responded to the ongoing climate change in various ways, such as the lengthening of growing seasons or changing vegetation composition (Mottl et al 2021). Accordingly, large efforts have been made to improve the monitoring and modelling of ecosystem responses to climate changes and extreme events over various biomes, ranging from tropical rainforests (Kim et al 2012, Longo et al 2019a to the Arctic tundra (Kim et al 2021, Larson et al 2022. Specifically, a lot of research is aimed at projecting terrestrial carbon cycles under the future climates based on different greenhouse gas emission scenarios, using numerical or empirical ecosystem models. The scales and scopes of studies have varied greatly: Euskirchen et al (2009) applied the Terrestrial Ecosystem Model (TEM; version 7.0) with nine climate scenarios from the International Panel on Climate Change, and suggesting potential changes in dominant plant function type (PFT) in northern Alaska by 2100. These changes were due to differential growth patterns among the PFTs in response to changes in light availability and growing season length, which resulted in increases in net primary productivity of all PFTs. Longo et al (2018) used the ecosystem demography model version 2 (ED2; Medvigy et al 2009, Longo et al 2019b to determine the sensitivity of carbon exchange alterations in terrestrial ecosystems to temperature changes. They projected that 25% of Amazon forests could lose biomass due to increases in the frequency of droughts under the highemission scenario representation concentration path-way8.5 (RCP8.5) by 2100.
Given the increases in model implementations, it has become essential to quantify projection uncertainties arising from various sources, such as climate inputs (e.g. carbon emission scenarios and climate models), ecosystem models (structure and parameters), and their interactions. Multi-model comparisons and multi-climate input comparisons are widely performed methods for uncertainty estimation by quantifying plausible ranges or distribution of outcomes (Scholze et al 2006, Booth et al 2012, Fisher et al 2014, Wu et al 2017, Fisher et al 2018, Park et al 2020. For example, Migliavacca et al (2012) compared phenological projections of temperate deciduous trees and found that different emission scenarios result in the largest differences in the phenological timing, gross primary productivity (GPP), and evapotranspiration (ET), rather than a phenology model structure or model parameters. Zhu et al (2016) showed that atmospheric carbon dioxide (CO 2 ) fertilization was the factor most significantly affecting the growth of global leaf area index (LAI) and carbon fluxes during 1982-2009 rather than climate change. However, most studies have estimated uncertainty as the range of the projected quantities; therefore, the contribution of individual factors (and their interactions) to the total uncertainty was not fully decoupled. Such uncertainty decompositions, however, have been intensively performed in hydrological studies. Bosshard et al (2013) and Aryal et al (2019) showed that the dominant uncertainty source on river streamflow projection varied by season: climate inputs (from different climate models) during summer and autumn (wet season) and hydrological models during winter and spring (dry season). Addor et al (2014) presented that climate models are the main source of uncertainties; however, the contribution of emission scenarios was predicted to increase toward the end of this century, while the importance of hydrological model varies depending on catchment elevation (i.e. more critical at higher elevation). Considering that the composition and distribution of PFT varies under different climate conditions due to competition among the PFTs for resources (i.e. light and water, Ahlström et al 2017), it is crucial to simultaneously investigate the uncertainty contributions from various factors, such as terrestrial model structure and parameterization, climate inputs, and their interactions. Such systematic decompositions of uncertainty sources are necessary for interpreting vegetation growth and the subsequent carbon cycle under the ongoing climate change.
In this study, we investigate how a transiting temperate-subtropical mixed forest on Jeju Island, South Korea would respond to climate change by the end of this century and also decompose the uncertainty contributions of input sources (i.e. emission scenario, climate models, model parameters). We implemented the ED2 with four behavioural parameter sets (PRMs) using meteorological projections from four global climate models (GCMs) under each of the four RCP scenarios (2.6, 4.5, 6.0, and 8.5). Based on these 64 simulation sets, we estimated the relative growth among the different PFTs, and longterm (2021-2099) trends in phenological changes (i.e. greenup, dormancy, and LAI), and carbon flux changes (i.e. net ecosystem productivity (NEP)). We then analysed the relative significance of each factor (i.e. PRMs, GCM, and RCP) to the total uncertainty in the projections.

Study site
The Jeju flux tower (33 • 19 ′ 4.80 ′′ N, 126 • 34 ′ 4.08 ′′ E) was established in May 2014 at the Jeju experimental forest on Jeju Island, South Korea (figure 1). This forest is located on the south-facing side of Mt. Hallasan at an elevation of 635 m above sea level with relatively moderate topography. The study site has a humid subtropical climate with annual mean temperature of 17 • C, total annual precipitation of 1920 mm (of which approximately half falls during the summer monsoon period between June and august), and prevailing winds from the northwest and northeast (Hong et al 2019). According to a field survey in 2008 (data available upon request), the soil type is silt loam with a pH of 5.29. This forest is dominated by several deciduous and evergreen tree species, including Quercus serrata, Quercus stenophylla and Castanopsis sieboldii.

Data overview
The diameter at breast height of approximately 1500 trees was measured in 2013 at nine plots (400 m 2 a plot) set at approximately 50-150 m away from the flux tower (figure 1). Only those trees with diameters greater than 2 cm were recorded. We classified the trees into five groups according to their favourable climate, phenology, and leaf type: temperate deciduous broadleaf (TDB), temperate evergreen needleleaf (TEN), temperate evergreen broadleaf (TEB), subtropical evergreen broadleaf (SEB), and subtropical deciduous broadleaf (SDB) (table S1). Over the nine plots, the dominant PFTs were the TDB and SEB, with average basal areas (BAs) of 24.71 and 11.34 cm 2 m −2 , respectively. The average BAs of the TEN, TEB, and SDB were less than 1.5 cm 2 m −2 (1.44, 0.57, and 0.2 cm 2 m −2 , respectively).
We used daily ET, GPP, and net ecosystem exchange (NEE) fluxes measured at the eddycovariance flux tower in 2015 (Hong et al 2019). The CO 2 concentration, wind velocity, sonic temperature, and humidity were measured using a 3D sonic anemometer (CSAT-3, Campbell Sci., Logan, UT) and closed-path infrared gas analyser (EC-155, Campbell Sci.), which were recorded using a data logger (CR-3000, Campbell Sci.) (Hong et al 2019). The half-hourly averaged data were generated after quality control processes, including double rotation, spike removal, and spectral correction (Hong et al 2019). The daily ET and NEE fluxes were generated by aggregating the half-hourly data after filling gaps using an artificial neural network (more details can be found at Hong et al 2019). The GPP flux was estimated by partitioning photosynthesis and ecosystem respiration (Lee et al 2021).
We also used gridded Moderate Resolution Imaging Spectroradiometer ( , with four radiative forcing cases from a low-emission scenario to a high-emission scenario (i.e. RCP2.6, 4.5, 6.0, and 8.5) at 0.25 • resolution from 2006 to 2099. For bias correction and spatial downscaling of the ISIMIP meteorological datasets, we employed the quantile mapping method using the biasCorrect function of the R package 'hyfo' (Xu 2020) after calibrating the biasCorrect function using the observation data (Automated Synoptic Observing System data between 2009 and 2014 provided by the Korea Meteorological Administration: www.kma.go.kr).

Model implementation and uncertainty analysis
In this study, a process-based model, ED2; (Medvigy et al 2009, Longo et al 2019b, Kim et al 2021, was implemented to investigate the long-term changes in vegetation growth and carbon cycles at the study site. The ED2 model is structured to represent the individual canopy, instead of conventional 'big-leaf ' approximation, in a heterogeneous forest. Therefore, the model is capable of simulating competition among species, including both deciduous and evergreen hardwoods with conifer trees, thereby allowing for the capturing the long-term dynamic processes, including forest succession (Medvigy et al 2009). The sensitivity of ED2 simulation to individual parameters was investigated (Feng et al 2018), but not the sensitivity to other factors (i.e. climate inputs from multiple GCMs under different RCP scenarios). In this study, we used the ED2 with the growing season index (GSI)-based phenology module from Kim et al (2021).
We initialized ED2 key variables related to vegetation structure (e.g. leaf and structural biomass and LAI) using the forest inventory data surveyed in 2013 with allometric equations in ED2. Soil composition was set following a field survey in 2008 at a 20 cm depth (sand:silt:clay = 36.1:62.5:1.4; data available upon request). Soil carbon was initialized based on the UN FAO Global Soil Organic Carbon Map (ITPS 2018), and soil nitrogen was set to the average value of the soil C/N ratios of warm temperate moist and wet forests (25.1 and 15.6, respectively; Post et al 1985). In addition to the existing PFT in ED2 (mid-temperate deciduous hardwood, i.e. TDB), we modified the model to add new PFTs-TEN, TEB, SEB, and SDB-based on the default parameters of TDB with appropriate phenology (i.e. no leaf-fall for evergreen type) and leaf structure (i.e. needleleaf for TEN). For the key variables related to vegetation physiology and growth, such as maximum carboxylation rate (µmol m −2 s −1 ), maximum and minimum photosynthetic temperature ( • C), specific leaf area (m 2 kg C −1 ), growth respiration factor (unitless), fine root allocation (kg kg −1 ), vapour pressure deficit (VPD) control on stomata closure (mol mol −1 ), and allometric parameters, we generated 1000 parameter sets randomly from prior distributions (the socalled Monte Carlo method). Due to the lack of field measurements related to the tree species at the study site, we assumed the prior distributions as the normal distribution having the same mean value as the ED2 default parameter with a standard deviation that is 10% of the mean value. We then implemented the ED2 daily from 2006 to 2014 with the 1000 parameter sets using the downscaled four GCM datasets under each of the four RCP scenarios. After stabilizing soil moisture for the first four years (2006)(2007)(2008)(2009), the simulation results for the following five years (2010-2014) were used to select four best behavioural parameter sets based on statistical measures (r and root-mean-squared-error) as compared to the MODIS LAI, GPP, and ET data (table S2). With the best behavioural parameter sets, we continued to implement the ED2 from 2015 to 2099 and further evaluated the performance of the ED2 using the MODIS LAI, GPP, and ET data (2015-2020), as well as the GPP, NEE, and ET flux data in 2015 (table S2; note that the NEP value of ED2 was compared to the negative NEE flux data). Using the simulation results, we investigated long-term changes in vegetation growth, phenological timings, carbon uptake, and net carbon exchange (i.e. LAI, GPP, NEP) under each RCP scenario from each GCM. In this study, we defined greenup and dormancy timings as when the GSI passes 0.15 during spring and fall, respectively, and peak growing timing as the GSI higher than 0.9.
We applied a three-way ANOVA, (Bosshard et al 2013) to investigate the major sources of uncertainty in vegetation growth (i.e. LAI) and carbon exchange (i.e. NEP) in the future projection. We subsampled all the simulation runs into iterations having two RCPs, two GCMs, and two PRMs. For the long-term annual analysis, annual measures (M f ; peak for LAI and sum for NEP) were calculated, and the climate change signal Y was retrieved against the annual measure (M c ) during the control period, 2010-2020 (i.e. Y = M f −M c ). For variance decomposition in the three factors (i.e. RCP, GCM, and PRMs), we calculated the total sum of squares (SS T , equation (1)), square deviations of the single effect of each factor, X (SS X , equation (2)), and interaction terms (equations (3) and (4)) for each iteration (i) where Y x1x2x3 is the climate change signal Y under the combination of each factor (x 1 , x 2 , and x 3 ) in the iteration i, Y o indicates averaging over the index. N 1 , N 2 , and N 3 are the total number of RCP, GCM, and PRMs in each iteration (i.e. two). For each effect (X), the fraction of variance contribution (η 2 ), ranging from 0 (no contribution to the total uncertainty) to 1 (sole source of the total uncertainty), was calculated (equation (5))

Long-term climate change
Based on the statistical analysis of the RCP scenarios from 2010 to 2099, the study site is expected to experience significant changes in climatic conditions (i.e. temperature, precipitation, photosynthetically active radiation (PAR), and VPD; figures 2 and S1). Annual mean temperature is projected to increase by at least 0.1 • C per decade (p < 0.001, RCP2.6) and up to 0.5 • C per decade (p < 0.001, RCP8.5). It is noticeable that the degree of increase differed considerably among the GCMs; for example, HadGEM mostly had the highest increases in temperature projected, while GFDL had the smallest increase projected under a same RCP. Annual total precipitation is projected to increase as well at rates between 1.7 and 37.4 mm per decade, but these increasing trends were only statistically significant under the RCP4.5 and 8.5 scenarios (p > 0.05 for RCP2.6 and 6.0). The differences in the precipitation projection were larger among the GCM than among the RCPs. The annual total PAR and the annual mean VPD are expected to increase significantly, regardless of the RCPs (p < 0.001), with the biggest increase under RCP8.5 (24.7 MJ m −2 yr −1 per decade and 0.43 hPa per decade). Specifically, these increases are predicted to occur primarily during the growing season (May-October), indicating that carbon uptake at the study site may be amplified during the growing season due to the intensified climatic conditions.

Forest type transition
Presently, study site is a temperate-subtropical mixed forest, with temperate species occupying 69.8% of the total BA in 2013 (figure 3). Our results showed that the projected climate favours the growth of the subtropical species and could, therefore, result in predominant growth rates in BAs by the end of this century. Specifically, the growth rate of TDB (presently the most dominant PFT) is 19.3 ± 11.5%, without noticeable difference regardless of meteorological intensification (i.e. higher RCPs). Meanwhile, the growth rate of TEB appears much higher-up to 67% under RCP6.0, but its proportion is still expected to be very low due to its currently-low presence (1.5% of the total BA). The growth rates of subtropical species (SDB and SEB) appear greater under higher RCP scenarios, ranging from 50.1 ± 27.5% under RCP2.6-88.1 ± 33.8% under RCP8.5, with large variations depending on GCM and PRMs. As a result, the dominance of subtropical species is expected to increase by 10.9%, from 30.2% in 2013 to 41.1 ± 3.1% by the end of this century. Moreover, due to the substantial growth of SEB, the composition of evergreen species is expected to become 42.5 ± 3.7% in 2099, about 11.4% higher than the current composition of 31.1%. The conspicuously high growth rate of SEB under RCP4.5, as compared to the rates under RCP2.6 and 6.0, could be attributed to the significant wetting trend of RCP4.5 (figure 2). The noticeable increase of evergreen species indicates that there will be substantial changes not only in the forest phenology but also in the year-long carbon flux pattern. It should be highlighted that the long-term growth of the BA is not solely attributable to climate but also to other factors such as seed dispersal from surrounding areas. Such effects, however, were not included in this study; thus, further investigations are required.

Phenological change
As a result of the growth of evergreen species, the increasing trend of the dormancy LAI appeared significant for most of projections (0.09 ± 0.08 m 2 m −2 per decade, figure 4(a)). The rate of change appears more sensitive to the PRMs and GCM rather than RCP. A similar tendency was also shown for the longterm trend of the peak LAI, with higher rates of change and larger differences among the projections (0.11 ± 0.08 m 2 m −2 per decade, figure 4(b)). The projected change under RCP2.6 by 2099 (8.8 ± 8.9%) was comparable to that from a multi-model ensemble (about 7.5 ± 2.5% for the 25-50 • N in the northern hemisphere, Piao et al 2020), while the estimation under RCP8.5 (12.3 ± 11.5%) is marginally similar to the multi-model ensemble (about 25.5 ± 8.0%, Piao et al 2020). The differences in the estimations could result from several factors, such as the site scale and model structures used (i.e. ED2 accounts for the resource competition among the PFTs, thereby resulting in less growth than those of the majority of other terrestrial models allowing single PFT in one unit). Meanwhile, the long-term trends of greenup were quite different to those of dormancy and peak LAI. It was clearer that the advance rate of greenup was greater under higher the RCP with less variations among the projections (0.37 ± 1.77 d per decade under RCP2.6 and 0.66 ± 0.29 d per decade under RCP8.5, figure 4(c)). We also found that almost half of the advance trends were not significant (p > 0.05). These results indicate that leaf emergence is largely sensitive to meteorology (i.e. RCP), yet highly uncertain. The start of dormancy appeared to be delayed for most of projections (0.06 ± 0.45 d per decade, figure 4(d)), which was lower than the rates of greenup advance. The dormancy trends were less dependent to the RCP scenarios and were rather driven by the PRMs and GCM. As a result of greenup advance and dormancy delay, the growing season length is expected to increase by 0.43 ± 0.99 d per decade on average (up to 0.65 ± 0.25 d per decade with RCP8.5).

Carbon flux change
Our results show that the study forest is likely to become a stronger carbon sink (i.e. uptaking more carbon) by the end of this century (figure 5). According to most of the projections, the annual NEP of the study site is expected to increase by the end of this century (figure 5(a)), as indicated by the growth of the BA and LAI (figures 3 and 4). The increasing trends in the annual NEP under the RCP2.6 and 4.5 were 18.5 and 23.8 gC m −2 yr −1 per decade, respectively (p > 0.05). Meanwhile, under the RCP6.0 and 8.5, the trends were significant at rates of 49.8 and 62.7 gC m −2 yr −1 per decade, respectively (figure 5(a)). This projection is similar to the estimation from a model ensemble value, where approximately 57 gC m −2 yr −1 per decade under the RCP8.5 was predicted (Cavaleri et al 2015). Our results also show that the long-term trends were more varied during the growing season, and these significant trends were more prevalent during dormancy ( figure 5(b)). This pattern of a wider range of trends during the growing season, as well as more significant trends during dormancy with less variation, indicates that there will be more uncertainty in the carbon cycle at the study site during growing season.

Uncertainty contribution
The characteristics of the uncertainties varied considerably depending on the variable being investigated. For the annual LAI, the contributions of three factors (PRMs, GCM, and RCP) to the total uncertainty at the beginning of the projections (2021-2099) were almost even ( figure 6(a)). The uncertainty contribution of PRMs gradually increased (13.7% in 2021), and the PRMs accounts for 36.1% of the total uncertainty in 2099 (60.2% when including the fractions of two and three interactions). Such dominant uncertainty contributions of the PRMs at the end of the century are expected to be slightly less during the growing season (April-September), as climatic factors (i.e. GCM and RCP) drive carbon uptake, thereby affecting leaf growth. This, in turn, could result in higher uncertainty contributions compared with those during dormancy ( figure 6(b)).
By contrast, 66.2% of the uncertainty in the annual NEP resulted primarily from the climatic factors (GCM, RCP, and their interaction), (79.7% when including interactions with PRMs, figure 6(c)). The contribution of each factor to the total uncertainty in NEP did not vary considerably over time or by season ( figure 6(d)). Those contributions of the climatic factors are expected slightly higher during the growing season (April-October, 80.7%), as compared to the contribution during dormancy (76.5%).
The difference in factors most strongly affecting uncertainty (i.e. PRMs for annual LAI, RCP for annual NEP) indicates that the strategy to improve projection accuracy could largely vary between variables.

Conclusions
The objective of this study was to investigate the response of a temperate-subtropical mixed forest to climate change (higher temperature, PAR, and VPD) by the end of this century. We implemented an individual cohort-based model ED2 with four behavioural PRM sets under the climate inputs from four GCMs based on four RCP scenarios and analysed the simulated long-term (2021-2099) changes in the vegetation species composition (i.e. BAs of PFTs), phenological timings (greenup and dormancy), and carbon flux (i.e. NEP). We further analysed the contribution of each factor (i.e. PRMs, GCM, and RCP) on the total uncertainty in the projected LAI and NEP, seasonally and annually. Our results showed (a) a gradual transition of PFT composition with growth of subtropical species (greater under higher RCP scenarios), (b) increases in LAI during both growing and dormancy seasons owing to the increases of subtropical evergreen species, and (c) that the study site will likely become a stronger carbon sink. The source of uncertainty differed considerably between the LAI and NEP; the main source of LAI uncertainty was the PRMs (and its interaction with climate inputs), whereas the climate inputs (GCM and RCP and their interaction) contributed the most to the uncertainty in the NEP projection. In this study, we focussed on comparing the results and uncertainties from different sets of RCPs, GCMs, and PRMs. Further studies should follow to assess the uncertainty contribution from other factors, such as warm drought (i.e. high temperature and VPD) and CO 2 fertilization, which have been significant for vegetation growth and physiological processes in the past (Zhu et al 2016, Huang et al 2019, Yuan et al 2022. This study highlights that it is crucial to consider the uncertainties resulting from terrestrial models and climate inputs to improve future projections of vegetation growth and the resulting changes in carbon cycles.

Data availability statement
All data that support the findings of this study are included within the article or references.