Rapidly changing high-latitude seasonality: implications for the 21st century carbon cycle in Alaska

Ecosys is an hourly time-step ecosystem model with multiple canopy, snow, and soil layers. The tightly coupled C, energy, nutrient, and water cycles operate across a wide range of temporal, and spatial scales. Key model processes are outlined below, and a comprehensive description is given in the supplementary material of Mekonnen et al (2019).

Each grid cell is forced with meteorological inputs and the model is initialized with prescribed seed densities for five plant functional types (PFTs) (deciduous, evergreen, sedge, moss, lichen). PFT-specific functional traits (i.e. CO₂ fixation kinetics, leaf optical properties, phenology, morphology, and root traits) result in emergent PFT variation in phenology, irradiance, CO₂ fixation rate, and water uptake. Functional traits differ for PFTs seeded in the boreal forest compared with those seeded in the tundra (figure 2). The model represents a multi-layer canopy, a residue layer, and multi-layer soil column. Leaf properties (including azimuth and angle), solar incidence angle, and light availability dictate propagation and absorption of direct and diffuse shortwave radiation in the canopy (Grant and Baldocchi 1992). Canopy energy balance, calculated using first-order closure schemes and by setting the sum of heat fluxes to zero, is tightly coupled to the water cycle through evapotranspiration. Heat and water transfers between the atmosphere, canopy, and soil column determine the temperature, water content, and ice content of each soil layer.

Zoom In Zoom Out Reset image size

Coupling between carboxylation and CO₂ diffusion determines canopy CO₂ fixation (Grant et al 2001). The Farquhar biogeochemical growth model (Farquhar et al 1980) is used to calculate carboxylation as the lesser of the light-limited and CO₂ limited carboxylation rates. Stomatal conductance balances water availability and maintenance of a fixed internal-to-ambient CO₂ ratio (Grant and Flanagan 2007). In the light limited carboxylation rate, potential electron transport rate is a hyperbolic function of radiation. Leaf temperature directly controls fixation rates through a modified Arrhenius temperature dependence of maximum rate constants, Michaelis–Menten coefficients, light-saturated electron transport rates, and an exponential temperature dependence of CO₂ and O₂ aqueous solubility. The modified Arrhenius temperature dependencies include low and high deactivation temperatures and a constant offset to account for thermal adaptation of each PFT. Air temperatures indirectly impact CO₂ fixation by modifying nutrient mineralization and uptake rates, boundary layer resistances, and evapotranspiration rates.

Photosynthesized sugars are allocated to non-structural C pools and oxidized to support first maintenance respiration and then shoot and root growth and active nutrient uptake. Oxidation rates and maintenance respiration requirements depend on a modified Arrhenius function of temperature that includes low and high temperature cutoffs and a thermal acclimation offset.

For winter deciduous PFTs, temperature-driven leaf onset (leaf offset) occurs after accumulated hours above (below) a set canopy temperature cross a set threshold. For evergreen PFTs, leaf dehardening (hardening) occurs after accumulated hours in increasing (decreasing) photoperiods crosses a set threshold. Leaf and root senescence for all PFTs occurs when non-structural C oxidation cannot meet maintenance respiration demands.

Pools of soil organic matter (woody litter, non-woody litter, particulate organic matter, humus, microbial biomass) are partitioned into subsets of varying susceptibility to hydrolysis (Grant 2001). Decomposition of these pools produces dissolved organic C, which drives microbial growth when oxidized. Substrate availability, soil temperature (according to an Arrhenius functional form with a constant offset and low and high temperatures of deactivation), and soil water content control oxidation rates. Oxygen and nutrient availability may impose additional constraints. Plant nitrogen (N) and phosphorus (P) uptake rates are affected by nutrient mineralization and immobilization rates driven by microbial biomass C:N:P ratios (Grant 2014).

The model was run at a 0.25° × 0.25° grid that covers Alaska. Clay and sand fraction, pH, cation exchange capacity, and bulk density were extracted from the Unified North America Soil Map (Liu et al 2013) and values for initial soil organic C content were extracted from the Northern Circumpolar Soil Carbon Database (Hugelius et al 2013). Surface air temperature, precipitation, incoming shortwave radiation, relative humidity and wind speed were taken from the North American Regional Reanalysis (NARR; Wei et al 2014) for the years 1979–2019. The first decade of the NARR record was used to spin-up the model over the years 1800–1978. NARR weather forcing for 2020–2100 was modified using seasonal anomalies from a CCSM4 ensemble member under the Representative Concentration Pathway 8.5 (RCP8.5). Since global carbon emissions are increasing at a rate consistent with RCP8.5, use of this high emissions scenario is common practice (Lee et al 2014, Lawrence et al 2015, Wieder et al 2015, McGuire et al 2018, Parazoo et al 2018). Historic CO₂ concentrations were used for 1800–2019, and CO₂ concentrations from RCP8.5 were used for 2020–2100. These simulations include N deposition taken from global spatially-distributed estimates (Dentener 2006, Wei et al 2014) and stand-replacing fire events, with frequency derived from the Mean Fire Return Interval (MFRI) dataset of the LANDFIRE product (Rollins 2009).

Ecosys representation of ecosystem C, nutrients, energy, and hydrological dynamics has been tested in many high-latitude sites. For example, modeled active layer depth matched long-term measurements at 28 Circumpolar Active Layer Monitoring sites (R² = 0.63; RMSE = 10 cm; Mekonnen et al 2021), modeled North American tundra gross primary production (GPP) matched upscaled EC tower measurements (geographically weighted regression, R² = 0.78; Mekonnen et al 2018), modeled tree composition of the Alaskan boreal forest agreed well with LANDFIRE–FCCS maps (Mekonnen et al 2019), and modeled NEE agreed well with EC tower measurements at 12 North American tundra and boreal sites (0.6 < R² < 0.9; Grant et al 2009, 2011, 2015, 2017). Additionally, ecosys accurately captured thermal and biological dynamics of short-term soil warming experiments at 4 sites across Alaska (Bouskill et al 2020). Fifteen studies of ecosys performance in high-latitude systems are described in the supplementary material (available online at stacks.iop.org/ERL/17/014032/mmedia).

Here we perform further validation of the ecosys model. Simulated soil temperatures at 24 locations were compared to data from the Soil Climate Analysis Network (SCAN; table S1; Schaefer et al 2007) and the Snow Telemetry (SNOTEL) Network (USDA Natural Resources Conservation Service 2020). Simulated NEE at 8 sites in Alaska was compared to data from Ameriflux EC towers (table S2). Where available, weather observations, rather than NARR forcings, were used to force the model for each Ameriflux site.

At the regional scale, we compared ecosys outputs with 7 observation-based estimates of monthly NEE across Alaska: (1) an estimate of NEE by Commane et al (2017) based on observed atmospheric CO₂ concentrations, remotely sensed data, and meteorological inputs; (2) NOAA's CarbonTracker (CT2019) estimates of NEE based on global measurements of atmospheric CO₂ concentration and an atmospheric transport model (Jacobson et al 2020); (3)–(5) three FLUXCOM estimates of NEE (FLUXCOM-RS-METEO-ERA5, FLUXCOM-RS-METEO-CRUNCEP, FLUXCOM-RS) based on machine-learning upscaling of global EC tower measurements using ERA5 weather forcing, CRUNCEP weather forcing, and no weather forcing, respectively (Jung et al 2020); (6) an estimate of winter CO₂ flux by Natali et al (2019) based on machine learning upscaling of site chamber, soda lime, and EC measurements; and (7) an alternate machine-learning upscaling of global EC tower measurements using ERA5 weather forcing produced by Zeng et al (2020).

Temperature and radiation limitations to plant productivity were quantified according to the methodology outlined in Donohue et al (2013), Keenan and Riley (2018), and Ukkola et al (2016). All daily modeled NPP values for the years 2010–2019 and 2090–2099 were combined and grouped by air temperature (incoming shortwave radiation) into bins of 1 °C (0.3 kWh m⁻² d⁻¹). For each bin the 99th percentile NPP was calculated. Using breakpoint regression analysis, air temperatures and radiation levels that limit NPP, and the lowest air temperatures and radiation levels that do not limit NPP, were calculated. Temperature and radiation scalars were linearly interpolated between these values.

In addition to the extensive validation of model performance in high-latitude ecosystems discussed above and listed in table S4, we also compared simulations of soil temperatures and NEE to site observations and regional observation-based estimates. At the site scale, ecosys soil temperatures at 5 cm depth agreed very well with measurements at 8 Soil Climate Analysis Network (SCAN) and 15 Snow Telemetry (SNOTEL) Network locations that are broadly representative of Alaskan climatic zones and land cover types (mean R² = 0.70 ± 0.11, mean bias = 0.09 ± 1.46 °C, and RMSE = 4.22 ± 1.41 °C; figure S4, table S1). We also found excellent agreement with ecosys monthly NEE and EC tower measurements at eight Alaskan Ameriflux sites (mean R² = 0.67 ± 0.11, mean bias = −0.09 ± 0.19 gC m⁻² d⁻¹, and Nash-Sutcliffe coefficient = 0.41 ± 0.56 gC m⁻² d¹; figure 2, table S2).

At the regional scale, we evaluated ecosys NEE against seven observation-based products (Commane et al 2017, Natali et al 2019, Jacobson et al 2020, Jung et al 2020, Zeng et al 2020). These products were generated by others, either through machine learning upscaling of site measurements or through estimation of land surface flux contributions to measured atmospheric CO₂ gradients using atmospheric transport models. We modeled an annual average NEE of −28 ± 25 gC m⁻² yr⁻¹ across Alaska for the years 2000–2019 (i.e. a net CO₂ sink from the atmosphere; uncertainty is expressed as standard deviation across years). The long-term mean modeled NEE is in excellent agreement with four of the six observation-based products that produced annual NEE estimates (figure 2(c)). Previous studies have expressed concern that process-model underestimation of fall and winter high-latitude C fluxes (Commane et al 2017, Natali et al 2019) has led to an incorrect characterization of the region as a C sink. However, ecosys predicted a 27 gC m⁻² (50%) larger combined fall and winter CO₂ flux than the average of the observation-based products, while still predicting that Alaska is currently a net C sink.

Under current climate conditions, modeled Alaska NEE for years 2000–2019 is largest in magnitude during summer (figure 3). Modeled summer dominance of the Alaskan C cycle is corroborated by the site measurements and 4 of the 6 observation-based products discussed above (figures 2(b) and (d)). The NEE seasonality of the observation-based products is broadly consistent with modeled NEE seasonality (table S3). These comparisons, and those described in the supplementary material, give confidence that ecosys is reasonably capturing the C cycle seasonality across our study domain.

Zoom In Zoom Out Reset image size

To assess how C cycle seasonality will change across Alaska with climate warming, we ran ecosys through the year 2100 using a CCSM4 RCP8.5 scenario. In spring (21 March–20 June), modeled NPP increases from 66 to 260 gC m⁻² season⁻¹ by year 2100 (positive NPP signifies positive plant growth; figure 3). This large increase in spring C fixation is driven by increases in air temperature (on average 6.7 °C by year 2100, figure S5) that lead to enhanced C fixation rates and earlier C uptake, particularly since temperature sensitivity of fixation rates is larger at lower temperatures (Kirschbaum 1995). Growing season onset (defined as the first day that modeled NPP is positive) in Alaska advances by 39 d by the year 2100 (5.8 d °C⁻¹; figure S7), consistent with published estimates of leafout advancement sensitivities in the northern hemisphere (Xu et al (2019): −1 to −4.5 d °C⁻¹; Piao et al (2015): −4.3 d °C⁻¹; Linkosalo et al (2009): −2.2 to −7.3 d °C⁻¹). Current observations of the effect of interannual variation in spring temperatures on high-latitude leaf emergence (Arft et al 1999, Pop et al 2000, Bjorkman et al 2020), growing season length (Chapin et al 1996, Keeling et al 1996, Myneni et al 1997), and plant productivity (Hicke et al 2002, Piao et al 2007) confirm that high-latitude plants experience severe temperature limitations during spring.

We quantified temperature and radiation limitations to modeled plant productivity using breakpoint regression analysis of daily air temperatures, incoming SW radiation, and modeled NPP (figure S6, Methods, Keenan and Riley 2018). According to this method, NPP is considered to be limited by, e.g. cold temperatures, if outlying NPP values increase at warmer temperatures. Under current Alaska climate conditions, modeled NPP experiences a 61% limitation due to air temperature in the spring. By year 2100, however, spring temperature limitation to NPP will relax by more than a factor of 3 (to 19%; figure S6). Since incoming shortwave radiation during these months is high, photosynthetic activity can respond positively to warmer spring temperatures. The geographically weighted regression coefficient (R² ) between spring temperature limitation and spring NPP remains high (0.83–0.85) throughout the century for the study domain (figure 4), providing confirmation that temperature remains a primary control of spring NPP throughout the century.

Zoom In Zoom Out Reset image size

Modeled spring microbial respiration (denoted as R_h in figure 3) increases from 48 to 116 gC m⁻² season⁻¹ throughout the century in response to earlier snowmelt and warmer soil temperatures. Soil temperatures increase more slowly than air temperatures in the spring (figure S5), so the increase in microbial respiration is much smaller than the projected increase in NPP. This difference results in a very large increase in the magnitude of spring NEE by year 2100 (−19 to −144 g C m⁻² season⁻¹; figure 3(a)).

In summer (21 June–20 September), modeled NPP increases from 201 to 344 gC m⁻² season⁻¹ throughout the century (figure 3(a)). Similar increases in spring and summer air temperature (∼6 °C–7 °C; figure S5) lead to similar increases in spring and summer Alaska GPP by year 2100 (326 vs. 331 gC m⁻²; figure 3). However, higher baseline air temperatures during summer lead to larger increases in autotrophic respiration (187 gC m⁻² in summer vs. 132 gC m⁻² in spring; denoted as R_a in figure 3), resulting in smaller increases in NPP compared to spring (figures 3 and S5).

Some studies have hypothesized that increases in growing season water stress may lead to a reduction in summer GPP (Buermann et al 2013, Liu et al 2020). Consistent with these studies, the impact of warming on summer water stress has been demonstrated using ecosys in an analysis of boreal forest dynamics (Grant et al 2009). In this study, we find that increasing summer water stress is buffered by increasing precipitation and increasing atmospheric CO₂ concentrations that reduce transpiration. The RCP8.5 climate forcing we applied here (Methods) predicts a 22% increase in summer precipitation across Alaska (figure S5), consistent with studies that suggest rainfall is projected to increase with warming at high latitudes (Bintanja and Andry 2017). As a result, we predicted slightly wetter summer soils, almost no change in summer evapotranspiration (0.4% decrease), and a small increase (∼16%) in summer water stress (calculated as the number of hours that the canopy water potential drops below a threshold value; Methods) across Alaska by year 2100.

Modeled summer microbial respiration increases from 114 to 221 gC m⁻² season⁻¹ throughout the century (figure 3). This increase is larger than that during spring because there are larger plant litter inputs and warmer and drier soils in summer than in spring. Summer net C uptake is projected to increase by only 36 gC m⁻², which is much smaller than the increase in spring net C uptake. By year 2100, summer net C uptake (123 gC m⁻² season⁻¹) is less than spring net C uptake (144 gC m⁻² season⁻¹) across Alaska (figures 3 and 5).

Zoom In Zoom Out Reset image size

In fall (21 September–20 December) and winter (21 December–20 March), there is only a small change in NPP throughout the century (21 gC m⁻² season⁻¹ for fall and 16 gC m⁻² season⁻¹ for winter). In fall, the current complete temperature limitation reduces to a partial limitation by year 2100 (figure S6). However, day lengths during fall in Alaska are very short, and there is not enough sunlight to drive photosynthesis. Since incoming shortwave radiation seasonality is driven primarily by earth-sun geometry (i.e. not climate change), the extreme radiation limitation to fall and winter NPP is not expected to change significantly over the coming century (Holland and Landrum 2015). In fall, temperature and radiation are each currently strongly limiting, but by year 2100 radiation will become the dominant limiting factor to fall C fixation. This prediction is consistent with observed large-scale increases in radiation limitation across northern latitudes (Zhang et al 2020).

By year 2100, average fall soil temperatures increase from −1.3 °C to 0.8 °C, average winter soil temperatures increase from −3.0 °C to −0.9 °C, and the first frost (defined as the date when surface soil temperature first drops below −0.2 °C) is delayed by one month on average. Modeled microbial respiration by year 2100 increases in response to warming soil temperatures from 25 to 61 gC m⁻² season⁻¹ in fall, and from 15 to 31 gC m⁻² season⁻¹ in winter. Net C loss during fall and winter is projected to shift from 76 gC m⁻² season⁻¹ (88% of summer net C uptake) to 163 gC m⁻² season⁻¹ (133% of summer net C uptake) over the course of the century (figures 3 and 5). Adding in the large increase in spring net C uptake, modeled annual Alaska NEE will increase in magnitude from −30 gC m⁻² y⁻¹ in the current decade to −108 gC m⁻² y⁻¹ by year 2100 (figure 5).

Whereas most large-scale land models link plant nutrient acquisition with instantaneous photosynthetic demand, ecosys allows plants to use nonstructural C reserves to uptake and store nutrients whenever they are available (Riley et al 2018, 2021). In the model, rates of mineralization and N fixation (symbiotic and non-symbiotic) control soil N availability and depend on soil temperature and liquid water availability. During fall and winter, modeled plant N uptake varies with the number of days that soil temperatures at 5 cm depth remain above freezing (R² = 0.60). Throughout the 21st century, the number of days above freezing increases from 21 to 60 d, synchronous with an increase in fall and winter plant N uptake from 0.26 to 0.69 gN m⁻² (figure 6(a)).

Zoom In Zoom Out Reset image size

Nutrient acquisition during the fall and winter has been shown to strongly influence year-round vegetation growth and competitive dynamics in northern ecosystems (Chapin and Bloom 1976, Larsen et al 2012, Malyshev and Henry 2012). By the end of the century, 22.5% of modeled plant N acquisition occurs during fall and winter, and it is likely that the large increase in spring productivity discussed above would not be possible without this source of N. Indeed, we find that modeled increases in spring NPP throughout the century are more strongly correlated with increases in fall and winter plant N uptake (R² = 0.69) than with increases in spring plant N uptake (R² = 0.47; figure 6). This result highlights the importance of accounting for fall and winter plant nutrient uptake in predictions of seasonal and annual high-latitude ecosystem response to climate change.

Some processes important to ecosystem C balance and export, such as topography, landscape-scale hydrology, thermokarst, and other geomorphological dynamics, are not represented in these model runs. Additionally, since boreal forest PFT species are not initialized in regions of present-day tundra, boreal treeline advance, which has been observed across the high-latitudes (Harsch et al 2009), does not occur in these simulations. However, trends in the seasonality of terrestrial ecosystem C exchange identified here are expected to be robust since they are attributed primarily to seasonal light availability driven by earth-sun geometry and large-scale shifts in seasonal temperature driven by climate change. As is the case for all model analyses of ecosystem dynamics, there is uncertainty associated with the gridded climate data and soil information we used to force ecosys (Mekonnen et al 2016, Wang and Clow 2021). Biases in temperature forcing data would affect the rates of processes (e.g. maximum fixation rates, electron transport rates) that drive modeled plant productivity and microbial respiration. While this uncertainty has an impact on comparisons between ecosys and observation-based products, it does not affect our conclusions, as we show using a sensitivity analysis (supplementary material).