The changing carbon balance of tundra ecosystems: results from a vertically-resolved peatland biosphere model

Observations from three eddy covariance towers located within the continuous permafrost tundra of the Imnavait Creek watershed, Alaska (68 37ʹN, 149 18ʹW) were used to evaluate vertically-resolved model formulation. Imnavait Creek is an ideal site for analyzing the patterns and drivers of Arctic carbon dynamics due to the abundance of meteorological, vegetation, and soil measurements available at the site (Kade et al 2012, Euskirchen et al 2017). The three towers, spaced ∼500 m apart, span a shallow topographic and associated edaphic gradient from 850 to 900 m with a heath tundra site on the ridge, a tussock tundra site on the slope, and a wet sedge tundra site near Imnavait Creek. The ridgetop has thin, well-drained soil with a shallow (2.3 cm ± 0.3) soil organic matter (SOM) layer, and the organic layer deepens down-slope with 11.7 ± 1.2, and 34 ± 2.6 cm of SOM, respectively, at the tussock, and wet sedge sites (Kade et al 2012). Further details of the study site and measurements can be found in (Euskirchen et al 2012, 2017, Kade et al 2012). During the study period, 2008–2016, the mean annual air temperature measured by the towers was −7.6 °C, and the mean annual precipitation provided by North American Regional Reanalysis across this period was 506 mm (Mesinger et al 2006). The site experienced substantial warming over the study period with the annual temperatures in 2014–2016 being 2 °C–3 °C warmer than 2008–2012 at all three tower sites (figure S6 available online at stacks.iop.org/ERL/17/014019/mmedia). This warming occurred predominantly during the winter months. There was also a significant increase in precipitation in 2013 and 2014 with approximately 50% more precipitation during the growing season (figure S6).

The simulations in this study were performed using the ecosystem demography model v2.2 (ED2) (Longo et al 2019). The ED2 model is a cohort-based terrestrial biosphere model that calculates the fluxes of carbon, energy, and water, all conserved quantities, between the vegetation, soil, and atmosphere (Moorcroft et al 2001, Medvigy et al 2009). ED2 simulates dynamic vegetation communities and accounts for competition between plant functional types (PFTs) and size, and has been developed and tested across a range of different biomes including boreal, temperate, and tropical forests (Medvigy et al 2010, Trugman et al 2016, Fisher et al 2018). In this study, we introduce a new vertically-resolved treatment of soil carbon processes that is critical for simulating the Arctic ecosystem. We also created new tundra PFTs; deciduous shrubs, evergreen shrubs, and graminoids, based on the prevailing vegetation types found at Imnaviat Creek to parameterize the vegetation formulation for an Arctic setting (see supplementary section S1).

The new soil biogeochemistry implementation developed in this study consists of a vertically-resolved representation of soil organic carbon (SOC) decomposition and accompanying rates of HR. SOC accumulates in the model from plant litter and is removed through HR based on a simplified version of CENTURY decomposition model in which the SOC within each layer is partitioned into three types; fast, structural, and slow, that respire at different rates (Bolker et al 1998, Longo et al 2019). Existing versions of the model use temperature and moisture values from the top 20 cm soil to determine the effects of soil temperature and moisture on the decomposition rates of carbon pools. In the new implementation, the three SOC pools are resolved vertically, and the HR in each discrete soil layer is then determined by the abundances and type of carbon in each layer, and the temperature and moisture conditions in each soil layer following the approach of Ise et al (2008). Each type of carbon is assumed to be input at the surface, and each day the depth of the carbon is recalculated based on the total mass of each type of carbon. To calculate the depth of the SOC, we assume a density for each carbon pool; we use average carbon densities of 29.1, 45.9, and 58.7 kgC m⁻³ for fast, structural, and slow carbon pools, respectively based on organic soil densities observed in the Alaskan tundra (Michaelson et al 1996, Yi et al 2009). A schematic of this model is shown in supplementary figure S3. The temperature and soil moisture sensitivities are plotted in figure S4. The decay rates and parameters relating the temperature and moisture dependencies of decomposition were adjusted based on an offline optimization to minimize the deviations between the simulated and tower NEP (see supplement S2). The optimal temperature in this respiration model is 35° C. The optimal relative moisture content, 0.75, is derived from the optimization.

Soil carbon in the model is dynamic, with carbon added at the top of the soil column and removed at each layer through HR. SOC is redistributed daily with SOC types preferentially moved down based on their density since deeper carbon tends to be denser and more recalcitrant (Yi et al 2009). In addition, soil properties, such as thermal conductivity that depend on the soil type, also dynamically change with depth based on soil carbon content. Following the method of Trugman et al (2016), the soil type transitions between mineral and organic peat depending on the soil carbon fraction. Ten transitional soil types are represented, whose properties are linearly interpolated between peat and mineral soil types defined in Longo et al (2019). The dynamic soil properties affect the profiles of soil temperature and moisture which can feeds back onto the vegetation.

The simulations at each of the three eddy covariance sites were driven by half-hourly meteorological data, including air temperature, pressure, precipitation, specific humidity, wind speed, and up- and downward short- and long-wave radiation. The simulations were initialized using information on vegetation composition from Kade et al (2012) and PFT-level measurements of above ground biomass (AGB) at nearby Toolik Lake for each vegetation community (Shaver and Chapin 1991, Chapin and Shaver 1996). The accompanying size-structure of each PFT was specified from the output from an equilibrium simulation in which the model was initialized with near bare-ground vegetation state and run for 200 years. To account for the discrepancies between the observed vegetation and the equilibrium simulation, the model was run for nine years using the observed vegetation mass and the simulated equilibrium size structure to allow the energy, carbon, and water fluxes to equilibrate. Total soil carbon at each site was initialized from Euskirchen et al (2017) with 17, 16, and 11 kg m⁻² for the wet sedge, tussock, and heath respectively. The soil carbon was partitioned into the three pools in the following ratios:1, 3, 6 for the fast, structural, and slow carbon based on the equilibrium ratios for peat found using an earlier version of the ecosystem demography model (Ise et al 2008). The three carbon pools are defined by their bulk density, decomposition rate, and sources and sinks. Soil temperature was initialized using observations from nearby borehole measurements at 0.34, 0.5, 0.9, and 3 m and soil moisture was initialized with field capacity (Euskirchen et al 2017).

In addition to simulations at the three tower sites, we performed a series of sensitivity tests to investigate how the model predictions vary in response to uncertainties regarding the initial soil carbon abundance, initial vegetation abundance, and surface water runoff times. Surface water runoff time, determines timescale at which surface water is reduced by 1/e, and is a proxy for slope in the model (Longo et al 2019). Runoff affects the water budget by removing surface water not immediately absorbed by the soil after precipitation events. In addition, we quantify the sensitivity of the model predictions to vertically resolving the respiration and allowing the soil texture to be dynamic based on the organic layer thickness (see table S3 and supplementary section S4).

The rate at which SOC is transformed into CO₂ via HR is determined by sub-surface soil moisture and temperature conditions. As figure 1 illustrates, there is a clear warming trend in soil temperatures over the analysis period. During the first five years (2008–2012), a strong seasonal cycle of soil temperature is apparent: surface temperatures are above freezing for a brief three-to-four month period during the summer, followed by strong seasonal cooling over the fall and winter in which cold temperatures (down to −10 °C) permeate to the bottom of the soil column (figure 1). The zero-curtain period is evident in the zero-degree contour line, which indicates thawed soil that persists well into the winter. However, from 2013 to 2016, the duration of the zero-curtain period increases and the magnitude of the seasonal cooling during the fall and winter weakens considerably (figure 1). For example, at 50 cm depth, the simulated thawed or mixed phase soil increases from about 180 days in 2013 to year-round in 2016 forming a talik. Furthermore, the coldest temperatures at 50 cm depth increase from about −6 °C in 2012 to 0 °C in 2016. This trend is driven both by warming air temperatures (figures S6–S9) and increasing snowpack which insulates the soil from the cold wintertime temperatures (figures S7 and S8). Consequently, the thaw depth (black contour line in figure 1) decreases from approximately 30 cm in 2008–9 to over 110 cm in 2015. This warming trend is corroborated by borehole measurements of soil temperature at Imnavait, which also show a substantial warming trend over this period, particularly during the winter months, consistent with our simulations: the 50 cm borehole minimum winter temperatures increase from −9 °C in 2012 to −3 °C in 2015 (Euskirchen et al 2017).

Zoom In Zoom Out Reset image size

Vertically resolving the SOC increases wintertime respiration at depth and decreases wintertime respiration at the surface compared to the original decomposition formulation in which decomposition rates vary in relation to the temperature and moisture conditions in the top 20 cm soil (figure 2). During summer, the vertically-resolved model has increased respiration near the surface and decreased respiration at depth compared to the original formulation (figure 2, compare panels (b) and (d)). Underlying this are large differences in the decomposition factor at depth in the original and vertically-resolved version of the model (figure 2, compare panels (a) and (c)). These translate into modest differences in HR rate (figure 2, panels (b) and (d)) however, because: (a) the deeper slow carbon is more recalcitrant, and thus has a much lower respiration rate; (b), below the organic layer, HR is zero due to the absence of SOC at these depths at the three flux tower sites. The difference in the magnitude and seasonal cycle of HR between the two decomposition formulations increases for soils with deeper organic layers because of the greater contribution from the deep thawed layers (figures S15 and 3(a)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Model sensitivity tests indicate that the amount, and thus depth, of soil carbon greatly affect the seasonal cycle and magnitude of HR and NEP (figure 3). Simulations initialized with soil carbon abundances from 5 to 135 kgC m⁻²—reflecting the typical range observed in the North American Arctic (Hugelius et al 2014)—indicate that abundant soil carbon leads to more HR and a smaller seasonal amplitude, i.e. relatively more HR during the fall and winter months (figure 3(a)). This altered seasonal amplitude arises because significant amounts of carbon now experience a different pattern of soil temperature and moisture. The pattern of increased fall and wintertime respiration and accompanying reduced seasonal amplitude seen in the vertically-resolved model in essence reflects the seasonal pattern of temperature at depth (figure 1). Furthermore, soil carbon abundance does not greatly affect plant growth, thus the effects on HR are also apparent in NEP (figure 3(b)). Soil carbon abundance has a strong effect on the net carbon flux and can determine the sign of the cumulative NEP, or whether the ecosystem is a net source or sink of carbon.

The decomposition factor in the top 30 cm responds strongly to the seasonal cycle of soil temperatures with notable departures due to changes in soil moisture (figures 4(a)–(c)). At depth, the decomposition factor is relatively high throughout the winter as deeper soil temperatures linger near the freezing point due to the release of latent heat as freezing fronts advance from the surface and permafrost (figure 1). Similarly, HR has a strong seasonal cycle in the top 20 cm; however, despite moderate decomposition factors, there is little HR deeper in the column due to the absence of organic material at these depths (figures 4(d)–(f)).

Zoom In Zoom Out Reset image size

Soil moisture accounts for more variation than soil temperature in HR across seasons and years, because HR is nonlinearly related to soil moisture—with decomposition rates being highest at a relative soil moisture of 0.75 (figures S4 and S5). Consequently, deviations both below and above this optimum, due to either drying or saturating of the soil, decrease the rate of decomposition and resulting levels of HR. For example, during the early summer, drying of the soils due to evapotranspiration reduces HR (see days 180–200 in figure 4, panels (d)–(g); see also figures S5 and S11). Consequently, timeseries of total HR at each site show substantial interannual variability (figure 4(g)). For example, during the winter of 2014–2015, HR decreased 85% in the wet sedge, 50% in the tussock, and 38% in the heath compared to the prior winter despite relatively warm winter temperatures due to saturated soils during this period (see figures S5 and S11). The 2014–2015 winter had the highest soil moistures of any year analyzed.

The temporal patterns of HR, discussed above, are an important contributor to the temporal trends in NEP at the three sites (figure 5, panels (a), (c), and (e)). As figure 5 illustrates, there is a strong seasonal cycle in the observed NEP with growing season peaks between 0.5 and 1 kgC m⁻² yr⁻¹ and overwinter rates of −0.1 to −0.4 kgC m⁻² yr⁻¹ with the heath having noticeably lower rates of wintertime respiration and the wet sedge having the highest wintertime rates (black lines in panels (a), (c), and (e)). At the wet sedge site (panel (a)), the model underpredicts the summertime peak NEP by approximately 0.4 kgC m⁻² yr⁻¹ likely due to low simulated primary productivity however, the model predictions are generally within 0.1 kgC m⁻² yr⁻¹ of the observed wintertime fluxes except for the winters of 2012–13 and 2014–15 in which the model predicted too much and too little HR respectively. As discussed above, these wintertime discrepancies were predominantly driven by changes in soil moisture.

Zoom In Zoom Out Reset image size

At the tussock site (panel (c)), predicted peak growing season NEP was within 0.1 kgC m⁻² yr⁻¹ of the observations except during 2012–13 and 2013–14 when summer maximum NEP was 0.5 kgC m⁻² yr⁻¹ below the observed NEP and wintertime NEP was 0.3 kgC m⁻² yr⁻¹ lower than observations. These discrepancies are due to anomalously high predicted rates of HR during these years as microbial activity responded to warming summer and winter soil temperatures and near-optimal soil moisture (0.75 relative saturation) conditions for decomposition (figures S5 and S11). At the heath site, the model overpredicted both the maximum NEP and the width of the growing season, suggesting an overly-long growing season phenology, potentially related to early offset and onset of snow cover in the model. Predicted wintertime NEP values are within 0.1 kgC m⁻² yr⁻¹ of the observations with the notable exception of 2012–14, when the model underpredicts wintertime NEP; however, the magnitude of the discrepancy during these years is lower than that seen at the tussock site.

As figure 5 illustrates, the tower measurements (solid black lines in panels (a)–(c)) exhibit a significant gradient in cumulative NEP across the sites: the heath site is a small net source, the tussock site is relatively neutral, and the wet sedge is a strong net source of carbon to the atmosphere. The model's NEP predictions (orange shading panels (a)–(c)) also vary across the three sites; however in the model predictions the upland heath site is a net carbon sink, the tussock site is a small net source, and the wet sedge site is a large source of carbon to the atmosphere over the observation period.

There is an inflection point in the predicted cumulative NEP at all sites after 2013. At the tussock and wet sedge sites, warmer and wetter conditions (see panels (a) and (b) in figures S10 and S11) increase HR substantially, resulting in increased carbon losses at these two sites as indicated by the declines in the orange curves in panels (b) and (d) of figure 5. In contrast, the heath site transitions from being a net carbon sink to a weak carbon source in 2013 and 2014 as indicated by the slight dip in orange curve in panel (f) of figure 5 during this period, before returning to a net carbon sink in 2015 and 2016. Examination of the effect of underlying temperature and moisture conditions on decomposition rates (figure S5) indicate that is due to relatively warm soil temperatures and soil moisture values reaching near optimal values for decomposition during this period.

Despite the close proximity of the three sites, there is distinct variability in their NEP and net carbon budget over this time period (figure 5). The variability in net ecosystem carbon change spans zero, with the heath site a carbon sink, the tussock site near neutral, and the wet sedge site a net carbon source (figure 5). A series of sensitivity tests exploring the key biological and physical differences between these sites (vegetation abundance, surface water runoff time, and soil carbon abundance) indicate that the abundance of soil carbon is the largest driver of the response of the NEP to climate due to increased HR (figures 3(a) and (b), figures S14, S17(a)–(d)).

The second largest driver of NEP variability between the study sites is the initialized vegetation abundance. However, the sensitivity of NEP to initial vegetation was weak compared to initial soil carbon (see figures S17 (i)–(l)). GPP varies significantly between the sites, with the differences predominantly driven by differences in initial soil carbon and vegetation communities (figure S12). Metrics of vegetation productivity and abundance (GPP, AGB and leaf are index) were predicted to increase across all sites over this time period as temperatures warm (see supplement S5, figure S19). Most of the increases in plant productivity come from enhanced shrub growth (figure S19). This pattern is consistent with observed 'Arctic greening' trends (Forbes et al 2010, Myers-Smith et al 2015, Martin et al 2017); however, other factors such as changes in snow cover, standing water, cloud cover can also cause changes in vegetation indices (Myers-Smith et al 2020).

Complex biophysical interactions between soil temperature, soil moisture and HR are a key control on the net carbon flux of tundra ecosystems. Our simulations indicate a thawed layer persists in the soil below a frozen surface for months into the fall and winter (figure 1), consistent with nearby borehole measurements (Euskirchen et al 2017). The active layer depth increased and the duration of the zero-curtain period extended between 2008 and 2017 due to warming air and soil temperatures over this time period. Given current climate trends in the Arctic region (Cohen et al 2014), we expect continued increases in the duration of the zero-curtain period as the climate warms and snow cover is reduced in the future (Collins et al 2013). A warming soil column could lead to more CO₂ emissions further enhancing rates of Arctic warming.

While soil temperatures tend to show relatively stable seasonal patterns (figure 1), soil moisture is typically more variable on short and long timescales (figure S11). For example, snowmelt and associated wetting of the soil can lead to rapid changes in the rates of HR in our simulations, which can increase the rate of HR in dry conditions or greatly reduce the rate of HR if the soils become saturated (figure S5). Deviations from the optimum moisture, both wetting and drying, can decrease the rate of HR. This complex relationship between HR and soil moisture can confuse attempts to attribute changes in the NEP, which are often attributed to temperature alone (e.g. Nottingham et al 2020).

The seasonality and type of precipitation are also important drivers of the Arctic carbon budget. Wintertime snow increases HR as deep snowpack acts as an insulator to cold wintertime temperatures (Liston and Elder 2006, Morgner et al 2010) and may reduce productivity by reducing the growing season (Galvagno et al 2013, Bieniek et al 2015, Scholz et al 2018), such as in the spring of 2014 and 2015 in our simulations (Figures S6–S9). These processes shift the carbon balance towards a net source to the atmosphere. This finding is consistent with a study of heath tundra in Greenland, which found that winter precipitation is correlated with HR and negatively correlated with GPP (Zhang et al 2018). More growing season precipitation may increase vegetation growth, however, this depends on the rate at which evapotranspiration and temperature increase (Wania et al 2009, Collins et al 2013). Increased precipitation could also increase soil moisture potentially inhibiting HR and CO₂ release (Chivers et al 2009, Zona et al 2009). Both of these processes would shift the carbon balance of the land and biosphere towards a net sink. Future increases in precipitation associated with global warming have an uncertain effect on the carbon budget due to these complex processes.

Translating soil temperature and moisture into realistic rates of HR requires vertically resolving soil carbon abundance and decomposition. The vertical structure of soil temperature and moisture determine the decomposition rate that acts on the vertical profile of soil carbon to determine the HR (figure 2). During the zero-curtain period, e.g. January (figure 2(b)), the rate of HR is higher using vertically-resolved soil dynamics. Capturing zero-curtain emissions arising from a subsurface thawed layer is important for correctly simulating the seasonal cycle and amplitude of NEP (Zona et al 2016, Commane et al 2017, Euskirchen et al 2017). The wintertime decomposition rate is much higher at depth, which translates into higher wintertime HR where soil carbon was present (figure 4). Due to the shallow organic layer depth, this effect was relatively modest in the tundra sites we simulated; however, sensitivity tests indicate that soils with deeper organic layers, such as can be found in boreal peatlands, had much higher wintertime respiration (see figures 3(a) and S15, panels (a)–(d)). Increasing rates of wintertime HR are consistent with recent observational studies in the Arctic (Natali et al 2019, Hashemi et al 2021).

Another important process affecting the dynamics of the Arctic carbon cycle is cryoturbation: the physical mixing of material in the soil column as a result of freeze/thaw processes. This process tends to vertically transport carbon deeper into the soil, omitting it may concentrate soil carbon near the surface altering rates of HR and the resulting carbon budget (Koven et al 2009). As seen in figure 2, carbon at the tussock site is confined to the top 40 cm of soil, which is within the seasonal active layer. Carbon transported deeper into the soil column would decompose more slowly thereby lowering the total annual HR, and also reducing the seasonal amplitude of HR because decomposition at depth extends into the fall and winter (figures 3 and 4). The effects of cryoturbation were not explicitly included in this study because rates of vertical mixing arising from cryoturbation are poorly constrained and spatially variable (Koven et al 2013). However, this is an avenue for future development that may improve the model representation of subsurface carbon abundance.

Previous studies have concluded that GPP dominates the variability in the growing season carbon budget (Shaver et al 2007, López-Blanco et al 2017, Schädel et al 2018). Simple models that only include LAI, light, and temperature can explain much of the NEP variability across tundra sites (Shaver et al 2007, Loranty et al 2011). However, these studies do not consider soil carbon variability as a predictive factor, and recent year-round measurements are challenging the notion that GPP is the primary control on annual NEP, i.e. whether the Arctic ecosystems are a net source or sink of carbon to the atmosphere (Euskirchen et al 2014, 2017, Zona et al 2016). GPP may be the largest source of interannual variability in areas of low soil carbon such as some upland tundra ecosystems. However, much of the Arctic has large stores of soil carbon (Crowther et al 2019), and in these areas the large soil carbon reservoir means that there is far larger potential for large changes in rates of HR. Small changes in the rate of respiration, due to a warming climate for example, applied to a large reservoir of carbon can greatly affect the net carbon flux. Both observations and model predictions show increases in wintertime respiration during the latter half of the study period (2013 onwards), without corresponding increases in NPP (Euskirchen et al 2017). We find that changes in HR predominantly determined whether an ecosystem was a net source or sink of carbon during this period (figure 5), and sensitivity tests illustrate that the leading driver of variation in HR and NEP between the tower sites is the soil carbon abundance (figures 3, S14–S17).

The relationship between soil carbon abundance and net CO₂ flux is often not emphasized in local scale studies; however, it is an underlying assumption in papers discussing long-term carbon release from permafrost loss (Schaefer et al 2014, Koven et al 2015, Schuur et al 2015, McGuire et al 2018, Wieder et al 2019). A recent statistical upscaling of net CO₂ fluxes across the panarctic found that SOC is the most important variable for predicting annual NEE (Virkkala et al 2021). We also find that the soil carbon abundance influences GPP (and NEP) through hydraulic feedbacks; in particular, the simulated over-productivity of the heath site (figure 5) results, at least in part, from differences in the hydraulic properties between organic and mineral soils. Reflecting observations, the heath site was initialized with approximately 30% less soil carbon, and our simulations show that springtime stomatal opening is substantially enhanced in soils with shallower organic layers that results in higher GPP (figures S12 and S13 panels (a) and (d)).

More generally our findings suggest that more accurate quantification of spatial variation in below-ground carbon stocks will be critical for improving predictions of how the carbon budgets of the Arctic ecosystems are responding to ongoing changes in temperature and precipitation. As illustrated in figure 3, the amount of soil carbon initialized in this study greatly affects the net CO₂ fluxes. The model sensitivity to the soil carbon stock is not surprising, as it is a first order term in the equation for HR (EQ S6). Model assumptions about soil carbon can determine whether the model predicts a net source or sink in the future (Wieder et al 2019); however, current models largely fail to reproduce observed values of soil carbon abundances (Huntzinger et al 2020). We circumvent this problem by initializing the model with observed soil carbon abundances. Patterns of interannual variability in the carbon budget are predicted to vary spatially, and differences in soil carbon stocks—rather than differences in vegetation—may be the leading driver of spatial variation across tundra landscapes (figure S17). This largely depends on the relative magnitude of the heterogeneity in soil carbon and vegetation abundance.

In summary, we find that tundra ecosystems that are both net carbon sinks (heath) and sources (tussock and wet sedge) to the atmosphere are strongly sensitive to changes in temperature and precipitation, and the magnitude of the net CO₂ flux increases with temperature. The effects of changes in precipitation on the tundra carbon budget are more complex however, and depend on both the season and form (rain vs. snow) of precipitation, and the existing moisture conditions of the soil. Between 2008 and 2017, we see a transition from a net CO₂ sink to a source in ecosystems with deeper organic soil layers during warmer and wetter years. This change is driven by a longer zero-curtain period and deeper active layer resulting in increased HR. Our results suggest that as the climate warms and wets in the future, the active layer will deepen and the zero-curtain period will extend, enabling more wintertime respiration (McGuire et al 2018, Natali et al 2019). As a consequence, ecosystems with shallower organic soils may continue to be net sinks while ecosystems with deep, organic rich soils may become carbon sources to the atmosphere.