The Coweeta Experimental Forest at the Coweeta Basin was established in 1934 by USDA Forest Service near Otto, North Carolina and was then renamed as Coweeta Hydrologic Laboratory in 1948 (35°03′ N, 83°25′W, Figure 1). The climate in this region is classified as marine, humid temperate, and features cool summers and mild winters, with less than 5% of annual precipitation occurring as snow (Swift et al., 1988). Two reference watersheds with detailed meteorology, vegetation, soil, and stream chemistry data were selected for this study: lower-elevation watershed 18 (WS18) and higher-elevation watershed 27 (WS27) (Figure 1). WS18 is a 12.5 ha watershed drained by Grady Branch, with an elevation range from 726 m to 993 m and a northwest facing slope of 52%. From 1986 and 2015, over half of the precipitation was exported in streamflow (average of 96 cm yr^–1 streamflow vs. 182 cm yr^–1 total precipitation). The average annual air temperature was 14.0^oC with an average minimum of 8.9°C and an average maximum of 19.6°C. WS27 is a 39-ha watershed drained by Hard Luck Creek, and its elevation ranges from 1,061 to 1,454 m with a northeast facing slope of 55%. From 1986 to 2015, over two third of precipitation was lost as streamflow (an average of 163 cm yr^–1 streamflow vs. 229 cm yr^–1 total precipitation). The average annual air temperature was 10.2^oC with an average minimum of 5.9^oC and an average maximum of 14.5^oC (derived from USDA Forest Service Southern Research Station² accessed on 08/01/2020). Soils are deep sandy loams underlain by folded schist and gneiss. Under the uppermost true and biologically active soils consisting of Ultisols and Inceptisols, there is a porous, friable, and unconsolidated saprolite layer above bedrock, which is believed to be a primary source of base flow and stream geochemistry (Velbel, 1988). Soil depth decreases with elevation with whole basin averaged soil depth around 3 m (Swank and Crossley, 1988). The dominant vegetation at both watersheds is southern mixed deciduous forests with overstory codominance by oaks (Quercus), maples (Acer), hickories (Carya) and tulip poplar (Liriodendron tulipifera L.) and an evergreen understory of rosebay rhododendron (Rhododendron maximum L.) and mountain laurel (Kalmia latifolia L.) (Day et al., 1988; Elliott and Swank, 2008).

These two watersheds have not been subject to human disturbances since the establishment of the experimental forest in 1934; however, they have experienced a variety of natural disturbances, including chestnut blight in the 1930s (Day and Monk, 1974; Elliott and Swank, 2008), fall cankerworm infestation in WS27 in the 1970s (Swank et al., 1981), extended drought in the 1980s (Elliott and Swank, 1994), and Hurricane Opal in 1995 (Elliott et al., 2002). In addition, the region has experienced increasing air temperature and extreme precipitation events over the years (Laseter et al., 2012).

PnET-BGC model has been applied, primarily in the northeastern region of the United States, to study forest primary production, water and nitrogen cycling, hydrology, as well as soil and lake/stream acidification in response to acidic deposition, climate change, and disturbances (Gbondo–Tugbawa and Driscoll, 2002; Chen and Driscoll, 2004, 2005a,2005b; Chen et al., 2004; Zhai et al., 2008; Wu and Driscoll, 2009, 2010; Zhou et al., 2015a; Pourmokhtarian et al., 2016, 2017; Robison and Scanlon, 2018; Valipour et al., 2018). Recently the model has been applied to watersheds in the southeastern and mid and western U.S. (National Smoky Mountain National Park, Zhou et al., 2015b, Niwot Ridge and Loch Vale Watershed, Colorado, and The H. J. Andrews Experimental Forest, Oregon, Dong et al., 2019). With the work presented herein, Coweeta Basin represents the southernmost U.S. application of PnET-BGC. This model includes two major modules: (1) PnET, a submodel that simulates water, carbon, and nitrogen cycling through vegetation, soil and drainage waters, and (2) BGC, a module that includes vegetation and organic matter interactions of abiotic soil processes, solution speciation, and surface water process involving other major elements (Gbondo–Tugbawa et al., 2001). The main state variables that are simulated include primary productivity, evapotranspiration, streamflow soil exchangeable cations, carbon and nitrogen mineralization, litter biomass, concentrations of sulfate, nitrate, potassium, calcium, magnesium, sodium, and acid neutralizing capacity (ANC) in stream at a monthly time step. We applied the model version that accounts for the impact of CO₂ on primary productivity based on the results from the Duke University’s Free-Air Carbon Dioxide Enrichment (FACE) project (Schlesinger et al., 2006) where gross primary productivity (GPP) increases with atmospheric CO₂ concentration until 600 ppm. Increasing atmospheric CO₂ has two confounding effects that are depicted in this model: an increase in maximum photosynthetic rate (A_max) and a reduction in stomatal conductance (g_s). Stomatal conductance and photosynthesis are coupled (Jarvis and Davies, 1998; Ollinger et al., 2002, 2009). Stomatal conductance changes in proportion to the difference between atmospheric CO₂ (C_a) and internal CO₂ (C_i) across the boundary of stomata. C_i is estimated from C_i/C_a (Ollinger et al., 2009, C_i/C_a ∼ 0.7 in Lavergne et al., 2019). The stomatal conductance and maximum photosynthetic rate are used to calculate CO₂ assimilation. Additionally, water use efficiency (WUE) is a function of CO₂ assimilation and vapor pressure deficit (VPD), adjusted by C_i/C_a and increases in A_max (ΔA_max). When atmospheric CO₂ increases, A_max will increase and g_s will decrease, which results in an increase of WUE. Plant transpiration is constrained by both ΔA_max and adjusted conductance (Delgs) when C_a is high and is calculated from GPP and WUE. We further modified the PnET-BGC model by adding a base flow component in order to obtain more accurate simulations of streamflow. The previous PnET-BGC model would simulate zero streamflow if no precipitation occurred in that month, which is not consistent with the observed streamflow pattern (such as December 2006 at WS18, and October 2000 at WS27). Based on observed non-zero streamflow when monthly precipitation was zero between 1936 and present, we added base flows of 0.6 cm/mo at WS18 and 0.9 cm/mo at WS27 to the streamflow originally simulated by the PnET-BGC.

In order to run the PnET-BGC, meteorology and atmospheric deposition data of major elements are required as inputs (Supplementary Tables 1, 2).

To reduce the data dimension of the many coupled vegetation-soil-stream processes studied, we applied principal component analyses based on 17 vegetation-hydrology-soil-stream variables simulated from 1931 to 2100 under both RCP4.5 and 8.5 climate scenarios with the PCA analysis in base R (https://www.r-project.org/, last accessed on 08/01/2020). The 17 variables considered are GPP, NPP, total litter mass, streamflow, evapotranspiration, soil base saturation, soil Al:Ca, net nitrogen mineralization rate, gross nitrogen mineralization rate, gross nitrogen immobilization rate, nitrogen uptake rate, concentrations of NO₃^–, SO₄^2–, Ca²⁺, Mg²⁺, and K⁺, and the acid neutralizing capacity (ANC) of streamwater (μeq/L). Note our application of PCA is for data dimension reduction and identification of important watershed process variables. Due to the exploratory nature of this analysis, normality of the data is not a strict requirement (Vaughan and Ormerod, 2005; Jolliffe and Cadima, 2016).

We evaluated the changes of the ecological processes/states under changing climate for spring (March to May), summer (June to August), fall (September to November), and winter (December to February). We tested whether the 17 variables simulated from the PnET-BGC model are significantly different under changing future climate (2071–2100) compared to the present (1986–2015) in each season using paired t-tests. We first checked the normality and temporal correlation (acf function in R) for annual time series (one seasonal average for each season of each year) between 1986 and 2015 and between 2071 and 2100, and found that data followed normal distribution and there did not exist temporal correlation, therefore justifying use of paired t-tests.

We focused on the variables that are important in Principal Component 1 and 2 and these include transpiration, GPP, NPP, soil base saturation, N mineralization, base cation (particularly potassium and calcium) concentrations, and ANC in streamwater.

Identifying a change point in the time series is the first step in identifying a potential driver of change, and, therefore a mechanism for a potential regime shift (Anderson et al., 2009), even though the existence of an abrupt change point does not necessarily lead to instability, hysteresis or a regime shift. To further investigate the impact of climate change on watershed processes, we conducted change-point analyses in the time series (1931–2100) of PCA score means of the first principal components (PC1s) that explained the largest variance of the variables (Anderson et al., 2009), to identify the years when PC1s exhibited or are projected to exhibit significant shifts, using the R-package of “changepoint” (Killick and Eckley, 2014). We have selected the binary segmentation (“BinSeg” in R) method. We checked for the normality between each change point (i.e., years) identified and they all met the normal distribution requirement (on annual basis). The same change point detection methods were used for the climate data between 1931 and 2100 with climate data after 2016 being the predictions for each RCP scenario (RCP4.5 and 8.5). For both PC1s and climate data (1931–2100) change-point detections, we implemented the program to allow the number of change points to vary from three to ten. However, our analysis consistently revealed only five change points.

For the two studied watersheds in the Coweeta Basin, observed air temperature increased from the past (1936–1965) to the current condition (1986–2015) and is projected to continue to increase between 2.5 and 4.8^oC from current conditions to the future (2071–2100) (Figure 2 and Table 1A). WS27 shows a slightly greater temperature increase than WS18 from current to the future conditions, particularly under RCP8.5 (+4.6^oC at WS18 vs. +4.8^oC at WS27). Under RCP4.5, temperature increase in spring is slightly larger than in summer followed by fall and winter (increase of 2–3°C except winter with an increase of < 2°C). Under RCP 8.5 summer has the largest temperature increase followed by spring, fall, and winter (increase of 4.5–5.5°C except winter with an increase of 4°C) (Table 1A and Figure 2A).

Observed precipitation decreased from past to current conditions during all the seasons except for fall. In contrast, future precipitation is projected to increase ranging from 2.2 to 4.0 cm/mo in all seasons in the future compared to current conditions with the fall showing the smallest increment (Table 1). Both air temperature and precipitation projections for Coweeta follow a similar pattern of regional climate projections for the southeast U.S. (Intergovernmental Panel on Climate Change [IPCC], 2014). Noticeably, the temperature and precipitation projections for each season under changing climate scenario (2071–2100) are statistically significantly different from those under current climate scenario (1986–2015) (Table 1B).

The PnET-BGC mostly effectively simulated monthly and annual vegetation-soil-stream processes at both WS18 and WS27 based on the stream chemistry data from 1970s to 2014 and streamflow data from 1930s/1940 to 2014, with NME ranging from –0.06 to 0.09 (Table 2). However, NO₃^– was underestimated at both WS18 and WS27 and SO₄^2– was overestimated at WS18. The model also underestimated the variance of NO₃^– and SO₄^2–. Measured concentrations of SO₄^2– at WS18 showed an increase between 1980 and 1990s that are seemingly inconsistent with long-term decreases in atmospheric deposition of S during the same period after implementation of the Clean Air Act (1970) and Clean Air Act Amendments (1990) (USEPA⁷ accessed on 08/01/2020). However, this pattern may reflect desorption of soil SO₄^2– to streamwater that was absorbed during prior decades of elevated atmospheric S deposition (Rice et al., 2014). The nitrogen cycle is more complex than other element cycles and it strongly linked to processes of vegetation and microbial activities, and therefore, PnET-BGC simulation of stream NO₃^– involves larger uncertainties than other solutes (Fakhraei et al., 2016; Shao et al., 2020). In addition, contrary to a decrease in both observed nitrogen deposition and simulated nitrate leaching since the 1970s, both watersheds show an unexplainable increasing trend of observed nitrate concentration in streams, which the model does not reproduce.

In addition to the long-term observed streamflow and stream chemistry concentrations, our model simulations were consistent with other published observed data (Supplementary Table 3), including NPP, soil base saturation, and net nitrogen mineralization rate. Model simulated average aboveground NPPs at WS18 and WS27 between 1971 and 1980 are close to published data from Day and Monk (1977) (8,056 vs. 7,965 kg/ha-yr). For net nitrogen mineralization rate, the 10-year averages (1971–1980) are 2.5 and 3.1 mg N/kg soil-mo at WS18 and WS27, which are within the range reported by Knoepp and Swank (1998) based on the samples collected along an elevation gradient (<1.2 – 3.8 mg N/kg soil-mo), and consistent with the pattern that the highest values occur at elevations over 1,000 m. WS18 at the lower elevation has higher soil base saturation (BS) of 24.9% than higher elevation WS27 with a lower BS of 4.0% with simulated data from 2001 to 2010 (Supplementary Table 3). This simulated range is close to the data range obtained from USDA Forest Services in 2008 soil survey from 9.2 to 19.7%.

The PCA analysis indicated that the first three principal components explained ∼ 85% of variance of the 17 variables analyzed (Supplementary Tables 4, 5 and Figures 3, 4). Particularly, PC1 explained 63.2 and 68.8% of total variances under RCP4.5 and RCP8.5 respectively at WS 18. PC1s explained 50.6% and 57.1% of the variances under RCP4.5 and RCP8.5, respectively, at WS 27. Based on the eigenvectors, the first principal components (PC1s) represent the contrast between primary production and transpiration vs. soil and stream ANC, and the second principal components (PC2s) depict nitrogen mineralization rate. The third principal components (PC3s) represent the concentrations of sulfate in stream at WS18 under RCP4.5 and at WS27 under both RCP4.5 and RCP8.5, while the PC3 at WS18 under RCP8.5 represents streamflow.

Based on the PCA analysis, the evaluation of the impacts of climate change on watersheds of Coweeta focused on the processes that are important in PC1s and PC2s during different seasons. Particularly we focus on seasonal changes of transpiration, NPP, and stream ANC in the main text. The figures for other process variables are described in Appendix Figures (Supplementary Figures 1–7). As transpiration, GPP, and NPP are negligible in winter, the probability density distributions are very concentrated, unlike in other seasons.

For both watersheds, transpiration (Figure 5) is projected to increase under changing climate. However, the increase in transpiration at WS18 is smaller than at WS27, probably due to higher water stress. Even though transpiration increases, annual streamflow (Supplementary Figure 1) under future climate is projected to increase with increasing precipitation (26.8 and 32.7% at WS18 under RCP4.5 and 8.5; 11.9 and 14.7% for WS27 under RCP4.5 and 8.5, respectively). Seasonally, winter is the season that shows pronounced increase at both watersheds.

The simulated impact of climate change on GPP is not as dramatic as for the water cycle (Supplementary Figure 2). However, it shows increases in spring and fall under both climate scenarios, increases in summer under RCP 4.5, but slight decrease under RCP 8.5. The increase of GPP at WS27 is more dramatic, especially in spring, summer and fall. Overall, annual GPP is projected to increase due to the extended growing season (2071–2100 compared to 1986–2015: increase by 13.8 and 17.9% for WS18, and 26.9 and 43.0% for WS27 under RCP4.5 and 8.5 scenarios, respectively). The simulated impact of climate change on NPP is moderate (Figure 6). On an annual basis, NPP is projected to increase although decreases are evident in the summer due to increases in respiration. Under RCP8.5, NPP shows decreases in summer and a larger variance in winter, with some high values. During summer, NPP is lowered in WS18 as much as 43.2 g biomass/m²-mo or 23.1% under RCP8.5. For WS27, simulated NPP increases during all seasons under RCP4.5. With RCP8.5, NPP shows decreases during summer. NPP is projected to increase to a greater degree in WS27 than WS18 due to less water stress associated with temperature increase.

At WS18, stream ANC is projected to decrease during all seasons under future climate with the larger decrease under RCP 8.5 at WS18 (Figure 7). The impact of climate change on stream ANC is minimal at WS27, with small decreases under RCP8.5.

Under climate change, the most conservative element, Cl^–, as well as basic cations (Na⁺, K⁺, Ca²⁺, and Mg²⁺) and SO₄^2– exhibited decreases in concentration but increases in total flux associated with streamflow increases. The flux of the base cations in soil is mainly associated with release from wood litter. Calcium concentrations (Supplementary Figure 3) in streamwater decrease dramatically in all seasons at WS18. A similar decrease is also evident in WS27, but not as marked. Potassium concentrations decrease (Supplementary Figure 4) in stream water during all seasons at WS18 with the most distinct change occurring in winter and spring. Decreases in potassium are also evident at WS27 in all seasons with the most distinct change occurring in winter. Sulfate concentration (Supplementary Figure 5) in stream water decreases in both watersheds, especially in winter under RCP8.5 and with lower values than projected under RCP4.5 at WS18. Compared to the conditions under current climate, future sulfate concentration will be less varied at WS27 but showing more distinct change than at WS18. WS18 has the highest median sulfate concentration in summer while at WS27, the highest concentrations occur in fall and winter.

Soil processes and properties are predicted to change as well, including soil base saturation, Al:Ca concentration in soil, and nitrogen mineralization rate. Soil base saturation is indicative of the soil acid-base status with lower values symptomatic of more acidic soil. High concentrations of Al:Ca in soil can be toxic to biota.

Soil base saturation (Supplementary Figure 6) shows decreases during all the seasons at both watersheds under RCP4.5 and 8.5. Decreases are pronounced under both climate scenarios at WS18, but only marginally under RCP8.5 at WS27. Al:Ca in soil shows a dramatic increase in all the seasons, but increases are most pronounced during summer at WS18 under RCP4.5 and 8.5. Al:Ca in soil increases particularly during summer at WS27. The simulated impact of climate change on N mineralization is small (Supplementary Figure 7). In WS18 under RCP8.5, the N mineralization rate shows some decrease. At WS27, the mineralization shows some decrease in winter but increase in summer.

Overall, the response of individual forest processes to climate change varies by season. Furthermore, vegetation seems to be more sensitive to climate change at WS27 than WS18, while soil processes and stream chemistry seem to be more impacted by climate change at WS18 than at WS27.

To more comprehensively present the impact of climate change, we summarized the changes of the 17 model-derived key watershed process variables in each season between current (1986–2015) and future (2071–2100) for both RCP4.5 and 8.5 scenarios, particularly, the directions of the changes and the significance levels of the statistical tests (Table 3).

In the vegetation processes, GPP and NPP are predicted to increase significantly except summer under RCP8.5. This can be explained by higher water stress during this season.

In the hydrology process, transpiration will significantly increase during almost all seasons as both precipitation and temperature are projected to increase from present day to the end of this century. Streamflow at the Coweeta basin is trending to increase in winter and spring seasons and decrease in summer and fall seasons, but the predictions in the future are not significantly different from current conditions except winter under RCP8.5.

The soil processes captured in base saturation across both watersheds and seasons are predicted to significantly decrease, which indicates a greater loss of base cations from the soil under changing climate. Along with the reduction of base saturation, the Al:Ca ratio is predicted to increase at all seasons under both climate scenarios, and the increase is statistically significant in all seasons. Additionally, some nitrogen processes are negatively impacted, such as gross nitrogen mineralization and immobilization. For the net nitrogen mineralization, results are complicated due to counteracting processes being involved at the same time.

In streamwater chemistry, the concentrations of major cations and sulfate are predicted to significantly decrease, even though soil losses of cations to nearby water flows (streams) are predicted to increase, probably due to the higher streamflow. ANC shows significant decrease in four seasons at WS 18 under both climate change scenarios and spring, summer, and fall at WS 27 under RCP8.5. The change of nitrate concentration mostly shows decrease.

We also found broad-scale climatic variability like El Niño-Southern Oscillation (ENSO) affected seasonal biogeochemical processes. Oceanic Niño Index (ONI that indicates ENSO) had a negative relation with the integrated biogeochemical processes (PC1) in the winter and early spring (January to April) at WS18 and WS27 (Supplementary Figure 8). Strong El Niño years generally had lower PC1 (Supplementary Figures 8, 9), indicating lower productivity and higher ANC in soil and streams, probably due to lower temperature in winter-early spring.

From the time series of the scores of the PC1s that explained more than 50% of the total variances, we derived five similar decadal change points for WS 18 and 27 (Table 4 and Figure 8). The change points occurred earlier at WS18 than at WS27 during 1980s and 2010s but are projected to occur at similar time points or even earlier at WS 27 relative to WS18 under future changing climate. This shift may imply that future climate might accelerate changes of biogeochemical processes at the higher-elevation WS27 site. The last change points are predicted to occur 20 years earlier under RCP8.5 (2059 at WS18 and WS27) compared to those under RCP4.5 (2079 at WS18 and 2074 at WS27) at the lower- and higher-elevation forests, respectively. The groups based on the time of change points of the PC1s confirmed their validity due to minimum overlaps between these different groups (Figures 3, 4).

When we compared the change points of temperature and precipitation to those of PC1, we found the first few change points of PC1 matched those of temperature at WS18, implying temperature change is likely an important driver of watershed response (Table 4 and Supplementary Figure 10). However, precipitation change appears to be the main driver for the change points of PC1 in 2079 at WS18. At WS27, precipitation is likely the main driver for the change points in 2010s and the 2070s, while temperature is the main driver for the change points in 1989, 2030s and around 2060. The climatic driver is not clear for the change point in 1969. In general, the last change points of PC1 under RCP4.5 (2070s) at both watersheds are likely driven by precipitation, while they seem to be driven by temperature under RCP8.5 (2050s–2060s).

In this analysis, we found vegetation is affected by climate change more at the higher elevation than at lower elevation, consistent with Lamprecht et al. (2018). The impacts of climate change on the hydrochemical processes have shown large spatial variability. Evapotranspiration and streamflow are predicted to increase under both moderate (RCP4.5) and high-end (RCP8.5) climate change scenarios at both WS18 and WS27. The predicted increase of evapotranspiration with climate change is consistent with patterns found in seven forest watersheds across four northeastern states in the United States (Pourmokhtarian et al., 2017). Annual streamflow, however, is highly variable among these watersheds, with some watersheds predicted to increase while others predicted to decrease, depending on whether the increase of evapotranspiration is offset by the increase of precipitation (Campbell et al., 2011). The hydrological predictions are different at the Andrews Experimental Forest in the northwestern U.S., where it is predicted that there will be a 20–71% decrease in annual transpiration under future climate change corresponding to 49–86% decreases in foliar biomass due to drought stress (Dong et al., 2019). The decrease in evapotranspiration could contribute to higher streamflow at this forest. In the western U.S. where future climate is projected to become warmer and dryer, gross and net nitrogen mineralization are predicted to increase, and NPP is also predicted to increase, due to feedbacks between warmer temperatures and enhanced nitrogen mineralization in forests. In the Coweeta Basin, Knoepp and Swank (1998, 2002) found that temperature and temperature-moisture interactions significantly affected net soil nitrogen mineralization based on a series of long-term studied plots at different elevations. When temperature was warmer and soil moisture content was higher, nitrogen mineralization was reported to be greater, consistent with our findings (Table 3). From the relatively low elevation oak-pine dominated watershed to the higher elevation watershed dominated by northern hardwoods, net nitrogen mineralization rate showed a strong spatial pattern with the highest rate observed at the highest elevation site, even though other factors (e.g., C:N ratio in soil) besides temperature and soil moisture influenced soil nitrogen processes both in the laboratory and in situ.

Climate change is projected to have variable effects on ecosystem responses across the seasons. Even though annual streamflow is projected to increase, these increases are less pronounced in the growing season. Spring may be impacted substantially because of loss of snowpack (Reinmann et al., 2018) and early leaf-out timing that can extend the growing season and exacerbate drought stress. A hotter summer together with a warmer spring, coupled with temporally unevenly distributed precipitation, could increase drought stress and risks of wildfire (Table 3). The drought stress could extend to other seasons. The fall season presents the lowest percent increase of precipitation among all the seasons. With higher evapotranspiration driven by higher temperature, streamflow is projected to decrease thereby causing more droughts. During winter, relative temperature (three times greater at WS27 in winter under RCP8.5) increases are particularly pronounced, which decreases precipitation inputs as snow and the magnitude and duration of snow coverage (Table 1A, Campbell et al., 2011; Pourmokhtarian et al., 2017).

Effects of climate change are difficult to evaluate during the “shoulder seasons” due to the large year to year variability. As a result, impacts during these periods are generally ignored, however, more and more research shows transitional seasons are affected by climate change to a greater degree than previously thought (Xie et al., 2015; Goss et al., 2020). Traditional definition of seasons and transition between seasons are predicted to become less clear as we may experience shorter spring and fall, while winter and summer will become longer, and most of the year will fall into more extreme hot and cold temperatures (Thomson, 2009). The impact of climate change on biogeochemical process during the shoulder transitions needs to be studied further.

Elevated CO₂ associated with climate change tends to increase plant WUE, i.e., increases primary productivity and reduces evapotranspiration (Reyer et al., 2015), which benefits forests by increasing their drought tolerance although there are counteracting effects to consider. Elevated CO₂ could also increase plant leaf area, therefore, making trees less drought tolerant (Ghannoum and Way, 2011). Reduced evapotranspiration could also lead to increases in leaf temperature, thereby increasing temperature stress. Increased productivity under elevated CO₂ may increase litter (Hyvonen et al., 2007), leading to increased litter accumulation (if decomposition cannot catch up) and higher vulnerability to fire. Many studies on the response of forest habitats to increased CO₂ are made at the small scale of individual trees or small plots. Scaling up the effect to the ecosystem scales involves considerable uncertainty. Furthermore, the fertilization effect of CO₂ is commonly found only in young trees (Korner et al., 2007). In our study, we assume elevated CO₂ increases primary productivity at the entire watershed without considering other counteracting effects (e.g., leaf temperature increase or age of trees), therefore our prediction of primary productivity under climate change is likely an optimistic estimate.

It is important to predict climate as accurate as possible in order to make valid predictions of ecosystem function change. This has been facilitated using downscaling method in our study. However, broad-scale climate variability like ENSO may enhance changes in temperature or precipitation (El Niño and La Niña Years and intensities at⁸ and accessed on 08/01/2020) at local spatial scales, and therefore has impact on biogeochemical processes as what we found out in this study. Generally, during El Niño years, the southeastern U.S. experiences cooler temperatures and wetter conditions in winter-early spring⁹ accessed on 08/01/2020. Together with the North Atlantic Oscillation (NAO), it was found that the southern Appalachians have higher snowfall and cooler temperature in winter during El Niño and negative NAO phases (Eck et al., 2019). In contrast, La Niña and positive NAO years result in warmer and drier winter weather. As the prediction of El Niño’s intensity using current climate models involves large uncertainty (Wang et al., 2019), predictions of temperature, precipitation and biogeochemical processes are particularly challenging under climate change. In order to reduce uncertainties of predicted biogeochemical processes, broad-scale climate variability needs to be accounted for in the climate change models.

PnET-BGC is a watershed-scale model that is operated on a monthly time step, so it is not well suited to address finer spatial and temporal resolutions (Merganicova et al., 2019). This model treats the vegetation in the watershed as a big leaf, therefore, it does not consider the difference among tree types and fine structures (Elliott et al., 2015), adaption (Jandi et al., 2019), or explicitly accounts for shift of vegetation types. Previous studies show that the critical transitions of vegetation type in mountain forests (Reich et al., 1995; Albrich et al., 2020) may change the fundamental relationships between foliar N concentration and maximum photosynthesis rate, which is one of the two fundamental relations that PnET-BGC is built on (Gbondo–Tugbawa et al., 2001). Furthermore, the function of understory shrub vegetation in cycling nutrients, carbon, and water (Rothstein, 2000) is another component not considered in this model. Different stages of vegetation not only have different capabilities of nutrient uptake, photosynthesis, cycling of carbon and water, but larger and older trees may be also more susceptible to stress such as drought and storms (Clinton et al., 1993, 2003; Clinton and Baker, 2000), which is not depicted in the model.

This PnET-BGC model requires a number of empirical data to calibrate, some of which are not available for the study watersheds (such as Al and DOC concentrations), making calibration for these parameters less robust. In addition, the model under-predicts nitrate concentrations in streamwater, suggesting the need to critically review this submodule for applications at Coweeta watersheds. Other simulations using PnET-BGC have been challenged to accurately simulate nitrogen dynamics (e.g., Shao et al., 2020). This model is designed to predict long-term hydrologic and biogeochemical function at the watershed scale (Aber and Federer, 1992; Gbondo–Tugbawa et al., 2001), combined with the lack of feedback between vegetation and climate, it is not an ideal tool to derive regime shifts. Other components lacking in the PnET-BGC further contribute to uncertainties in the simulations. These include the lack of multiple soil depths (two-layer model available but harder to calibrate so we used one-layer model), and lack of consideration of phosphorus dynamics.

Despite these aspects that can be improved, PnET-BGC provides a useful tool to evaluate the impact of climate change on forested ecosystems over long time scales (multiple decades to centuries). When depicting complex ecosystems, as many other ecological models have done, it is necessary to make assumptions and acknowledge inherent uncertainties (Aber and Federer, 1992). In order to capture the uncertainty, we have applied two climate change scenarios. The consistent response of biogeochemical cycles provides higher confidence in our predictions on the climate change impact.

This research will provide useful information for resource managers in anticipation of the potential risks forests may face under climate change. Forest managers can consider possible measures such as tree planting and forest restoration, wildfire control, establishment of mixed stands, planting better adapted species or varieties to mitigate adverse impacts in the future. For monitoring, the forest managers may consider watershed specific and seasonal specific evaluations instead of focusing on growing seasons or one process in soil or vegetation only.

We are facing an uncertain future based on (1) the greater uncertainty in prediction for precipitation than for temperature (Intergovernmental Panel on Climate Change [IPCC], 2014), and (2) precipitation will be the main driver for future change points under RCP4.5. Considering the change that may occur in high-elevation forests directly influences ecosystems downslope and downstream and the highly uncertain future projections, measures should be taken to protect, conserve and restore these high-elevation forests. The improved understanding of forests’ non-linear responses to a changing environment gained from this research can potentially guide decisions in forest management even with large uncertainties, for example, detection of early warning signs of potential regime shifts and selection of species planted in restoration.