21st century tundra shrubification could enhance net carbon uptake of North America Arctic tundra under an RCP8.5 climate trajectory

Recent observed shifts in Arctic tundra shrub cover have uncertain impacts on 21st century net ecosystem carbon exchanges. Here we applied a well-tested ecosystem model, ecosys, to examine the effects of North America Arctic tundra plant dynamics on ecosystem carbon balances from 1980–2100 under the RCP8.5 scenario. Tundra productivity was modeled to increase from enhanced carbon fixation and N mineralization under recent and future climates. Between 1982 and 2100 and averaged across the region, predicted increases in relative dominance of woody versus non-woody plants increased ecosystem annual net primary productivity by 244 g C m−2 that offset concurrent increases in annual heterotrophic respiration (139 g C m−2), resulting in an increasing net carbon sink over the 21st century. However, smaller increases in seasonal carbon uptake during winter (1 g C m−2) and autumn (22 g C m−2) and greater increases in ecosystem respiration (winter (23 g C m−2) and autumn (47 g C m−2)) by 2100 versus 1982 resulted in larger carbon losses during these seasons that completely offset the gains in spring (13 g C m−2) and 25% of the gains in summer (140 g C m−2). Modeled soil temperatures were predicted to increase more slowly than air temperatures (~0.6 °C for every 1 °C increase in air temperature over the 21st century). This slower soil versus air warming, and thus greater increases in CO2 fixation versus soil respiration rates, also contributed to the tundra remaining a carbon sink through 2100. However, these higher gains versus losses of carbon may be a transient response and not sustainable under further soil warming beyond 2100. Our modeling analysis allows us to extend beyond results from short-term warming experiments, which cannot characterize effects associated with decadal-scale changes in plant communities.


Introduction
Climate change in northern ecosystems will affect the rate and duration of carbon fixation and ecosystem respiration (R e ) (Albert et al 2011, Klady et al 2011. The extent to which future climate change affects the net carbon exchange of these ecosystems is uncertain (McGuire et al 2009). Several warming experiments (Hill and Henry 2011, Hollister et al 2005, Klady et al 2011, Oberbauer et al 2007, Sistla et al 2013, Walker et al 2006 reported increases in productivity from enhanced nitrogen (N) mineralization (DeMarco et al 2014, Salmon et al 2015). Increases in ecosystem productivity from extended growing season length in the northern ecosystems were also reported by several studies (McManus et al 2012, Myneni et al 1997, Olthof et al 2008, Tucker et al 2001, Verbyla 2008, Zhang et al 2008. Concurrently with overall increases in ecosystem productivity, warming in higher latitudes increased permafrost thaw (Jorgenson et al 2001, Lantz and Kokelj 2008, Lawrence et al 2008, Nowinski et al 2010, Schuur et al 2008 Oechel et al 1993, Schuur et al 2008. Land surface models have reported contrasting responses of northern ecosystems to warming (e.g. Zhuang et al (2006) and Qian et al (2010) predicted these ecosystems to be a net source and sink, respectively) attributed to differences in the rates of carbon gains versus respiration losses.
In a meta-analysis of 61 tundra sites experimentally warmed for up to 20 years, Elmendorf et al (2012) showed that the responses of tundra plants to warming varied with site conditions such as soil moisture and plant fuctional types (PFTs). In a nine year warming experiment in three high Arctic tundra ecosystems, Welker et al (2004) observed contrasting responses of increasing carbon uptake versus ecosystem respiration (R e ) among sites with different soil moisture and soil organic matter stocks. In a meta-analysis of 32 northern ecosytem experiments with two-nine years of warming, Rustad et al (2001) found contrasting responses of net ecosystem C exchanges that varied with site conditions. A meta-analysis of International Tundra Experiment (ITEX) warming experiments (Bouskill et al 2014) indicated concurrent increases in carbon uptake and soil respiration with warming, and showed that two Earth System Model (ESM) land models were unable to accurately represent these responses.
Although these tundra warming experiments provide valuable warming scenarios, they were limited in numbers and time (Hollister et al 2005, Rustad et al 2001. The responses of these experiments were largely dependent on site conditions, thus the longterm responses across the broader spatial domain of Arctic tundra may not be captured , Hollister et al 2005, Lamb et al 2011. In particular, these experiments cannot fully represent the warming effects associated with relatively slower changes in species composition and abundance (Rustad et al 2001, Shaver et al 2000. Several studies have shown recent increases in shrub growth and abundance in the Arctic tundra from repeated photography (Tape et al 2006, Tremblay et al 2012 and long-term plot based warming experiments (Chapin et al 1995, Wahren et al 2005, Walker et al 2006. These increases in shrub growth may further enhance the ecosystem carbon sink due to increasing woody carbon stocks (Leffler et al 2016, Sistla et al 2013 with higher carbon to nitrogen (C:N) ratios and longer turnover times, which can result in greater carbon gains per N invested.
To characterize these complex interactions, which is difficult with experimental manipulations and observations alone, and predict future trajectories of North America Arctic tundra net ecosystem carbon exchange, we applied a well-tested ecosystem model, ecosys (Grant 2014, Grant et al 2015, to examine the effects of North America Arctic tundra plant dynamics on ecosystem carbon balances from 1980-2100 under the RCP 8.5 scenario. The model mechanistically represents key biological, physical, and chemical processes that control long-term carbon cycle dynamics such as internal resource allocation and remobilization; soil thermal and hydrological dynamics; and microbial soil carbon, nitrogen, and phosphorus transformations.

Model description
Ecosys is an hourly time-step model with fully coupled carbon, energy, water, and nutrient cycles (Grant 2001, Grant et al 2015. A detailed description of inputs, parameters, and algorithms used in ecosys can be found in the supplement I available at stacks.iop.org/ERL/13/054029/mmedia. We briefly review here the processes important for controlling terrestrial gross primary productivity, ecosystem respiration, and PFT dynamics in tundra ecosystems.

Gross primary productivity (GPP)
The CO 2 fixation rate is controlled by coupled schemes for gaseous diffusion and biochemical fixation as affected by plant water and nutrient status. The total CO 2 fixation from each leaf surface in the model results in GPP of each plant population. GPP is controlled by plant water status calculated from convergence solutions that equilibrate total root water uptake with transpiration (Grant et al 1999). Carboxylation is directly affected by canopy temperature (T ) modeled through the Arrhenius functions for light and dark reactions (Grant et al 2007). The model uses parameters for temperature sensitivity of key CO 2 fixation processes from Bernacchi et al (2003) for temperatures from 10 • C-40 • C and additional parameters for low and high temperatures inactivation by Kolari et al (2007). Carbon uptake is also strongly controlled by plant N and phosphorus (P) content. Increased soil temperature (T s ) enhances plant productivity by increasing soil N and P mineralization and root and mycorrhizal active uptake (Grant 2014). Changes in air temperature and precipitation affect GPP directly through its effects on carboxylation, oxygenation and indirectly thought its effect on soil-water-atmosphere water relations (Grant et al 2007). Leaf onset (leafout in deciduous, dehardening in evergreen) and termination (leafoff in deciduous, hardening in evergreen) was modeled from number of hours of canopy temperatures accumulated above or below set values during lengthening or shortening photoperiods, respectively.

Ecosystem respiration (R e )
Ecosystem respiration is calculated from autotrophic (R ) and heterotrophic (R h ) sources. Canopy and soil temperature and water contents are calculated from surface energy and water exchanges coupled with soil heat and water transfers through atmosphere-canopysnow-surface residue-soil profiles (Grant et al 2012).
Temperature-dependent oxidation of non-structural products of CO 2 fixation (R ) drives R by all branches, roots, and mycorrhizae. R is first used to meet maintenance respiration requirements (R ). Excess R over R is expended as growth respiration R , constrained by branch, root, or mycorrhizal turgor potential. When R exceeds R , the shortfall is met by the respiration of remobilizable C and translocation of associated N and P in leaves and twigs or roots and mycorrhizae and the loss of associated non-remobilizable (i.e. structural) C, N, and P as litterfall.
R h of each organic matter-microbe complex (coarse woody litter, fine non-woody litter, manure, particulate organic matter and humus) represented in ecosys is determined by the active biomass (M) of heterotrophic microbial populations and the substrate concentration (Grant et al 2006a). R h is controlled by T s through an Arrhenius function and by soil water content through its effect on aqueous microbial concentrations [M]. Decomposition generates dissolved organic carbon (DOC) that drives microbial growth. R h is also controlled by microbial N and P concentrations, DOC, T s , available O 2 , and soil water potential. Concentrations of C, N, and P in roots and mycorrhizae drive exudation of nonstructural C, N, and P to DOC, dissolved organic N (DON), and dissolved organic P (DOP) in soil. R h drives CO 2 emission from soil through diffusion and volatilization in aqueous and gaseous phases (Grant et al 2012). Changes in soil ice content are used to calculate active layer depth (ALD, defined as maximum annual thaw depth), modeled from the general heat flux equation driven by surface energy exchange and subsurface heat transfer (Grant and Pattey 1999). Net ecosystem productivity (NEP) is determined from the difference in ecosystem net primary productivity (NPP = GPP − R ) and R h .

PFT competition and dynamics
Ecosys represents multiple canopy and soil layers. The vertical profiles of canopy leaf area and root lengths are prognosed from allocations of plant nonstructural C, N, and P to each organ of each PFT (Grant 1994, Grant and Hesketh 1992, Grant et al 1989. Thus, each PFT competes for irradiance, water, and nutrients within each canopy and rooted soil layer depending on leaf area and root length. Light interception of incoming direct and diffuse radiation and back scattering is resolved across each canopy layer. Each PFT competes for nutrient and water uptake from common nutrient and water stocks held across multilayer soil profiles, calculated from algorithms for transformations and transfers of soil C, N, and P, and for transfers of soil water (Grant 2016, Grant et al 2003, Grant et al 2007. Nutrient and water uptake of each PFT depends on root length and density of primary and secondary root axes, driven by the allocation of non-structural C, N, and P to each axis (Grant 1993, Grant et al 1989. Allocation rates are determined by non-structural C, N, and P concentration gradients within each PFT arising from leaf CO 2 fixation, and by root N and P uptake versus consumption by R a in each axis determined by its nutrient and water status. Modeled differences in PFT functional traits determine the strategy of resource acquisition and allocation that drive growth, resource remobilization, and litterfall, and therefore each PFT's dynamic competitive capacity under changing growing conditions. There are several ecosys PFT-specific functional traits important for predicting high-latitude vegetation competition under a changing climate (i.e. CO 2 fixation kinetics, leaf optical properties, phenology, morphology, and root traits). These differences in plant traits result in emergent PFT variation in phenology, irradiance, CO 2 fixation rate, and water uptake and thereby each PFT's competitive ability.

Simulation design and testing
Plant establishment occurs from prescribed initial seed densities for coexisting PFTs (deciduous shrubs, evergreen shrubs, sedge, moss, lichen) across the simulation spatial domain. The model was initialized with soil attributes including layer depth, clay and sand fraction, pH, cation exchange capacity, bulk density, and soil organic carbon stocks (

Changes in tundra ecosystem productivity
The predicted spatial pattern of long-term mean annual  GPP agreed well with EC-upscaled values (Geographically weighed regression, R 2 = 0.78; figure 1), although we caution that the EC upscaling is based on very few observations in this region (Jung et al 2011). In particular, both the modeled and EC-upscaled tundra GPP have larger values in the south and southwest and smaller values in the high Arctic and Arctic cordillera in the northeast. Long-term (1982Long-term ( -2010 mean annual modeled GPP was 0.41 Pg C, similar to the estimates from EC-upscaled GPP with 0.44 Pg C. GPP, NPP, R h , and NEP of the North America Arctic tundra were all modeled to increase with climate change (figure 2). These increases in tundra carbon uptake were modeled from more rapid carboxylation kinetics (Bernacchi et al 2001) due to higher temperature sensitivity at lower temperatures (Sjögersten and Wookey 2002) and facilitated by enhanced N uptake from soils due to increasing ALD and more  Oberbauer et al 2007) that reported increases in northern ecosystem productivity. The modeled changes in NPP and NEP to recent warming were spatially heterogeneous. The spatial average modeled NPP increased by 27% and NEP by 62% from 1982-2010. Observations from long-term plots and experimental warming were also shown to have diverse and contrasting responses depending on the moisture regime and PTFs (Oberbauer et al 2007, Walker et al 2006. However, similar to our modeled results, most tundra plot based experiments have shown overall increases in productivity. For instance long-term biomass measurements from long-term plots in a Canadian High Arctic tundra site indicated 158% aboveground and 67% root biomass increase between 2005 and the 1980s (Hill and Henry 2011). Results from long-term ) experimental plots at Ellesmere Island, Nunavut, Canada (Hudson and Henry 2009) have also shown that aboveground biomass increased by ∼160%. Similarly, in a meta-analysis of 23 sites warmed from 2 • C-5 • C in green house warming experiments across the Arctic, Dormann and Woodin (2002) reported 120% average increases in biomass compared to the control plots.
Over the 21st century, the modeled growing season extended by ∼60 d, spatially averaged across the tundra and primarily in the spring, resulting in increased NEP (figure 3). Although increases in R e were also modeled with warming (figure 2), the carbon uptake rate was greater than the losses, resulting in an overall carbon sink in the North American Arctic tundra (figure 2).  varied spatially, with greater increases in much of the low Arctic and Alaska, despite localized declines in a few parts of the high Arctic and southwest Alaska (figure 4(a)). Consistent with these results, Sitch et al (2007), using a land surface model and atmospheric inversions, reported that the Arctic tundra has been a carbon sink in recent decades. Modeled annual NEP continued to rise over the 21st century, resulting in much of the North America tundra remaining a net carbon sink (greater in the low Arctic versus high Arctic; figure 4(b)).

Changes in modeled annual NEP in recent decades
Predicted 21st century northern ecosystem responses to climate change vary widely among land surface models. Consistent with our result, Qian et al (2010) predicted that northern ecosystems (above 60 • N) will be a carbon sink through the 21st century based on 10 models from the Coupled Carbon Cycle Climate Model Intercomparison Project. However, these 10 models had a wide range of prediction of carbon fluxes (e.g. NEP range from −0.2 to 1.0 Pg C yr −1 by 2100). In another carbon-nitrogen modeling study, Koven et al (2015) predicated that the permafrost regions of northern ecosystems will remain a carbon sink through 2100. On the other hand, Zhuang et al (2006) used a process based model and predicted that the northern ecosystem (above 50 • N) will become a carbon source by 2100.
The increase in woody shrub NPP was modeled to be the largest factor contributing to enhanced net ecosystem carbon uptake by 2100 ( figure 5(a)). Nonwoody plants were modeled to remain the dominant PFT over woody plants until ∼2045, after which woody plants dominated through 2100. The spatially averaged NPP:R h ratio of non-woody plants was greater than that of woody plants in recent decades and through ∼2045. The modeled gains in non-woody plant NPP offset ∼65% of R h from 1980-2045 ( figure 5(b)), suggesting that gains in non-woody plant carbon uptake were insufficient to offset increasing R h when the tundra was mainly dominated by non-woody plants. The non-woody plant NPP:R h ratio diminished to ∼0.35 by 2100 ( figure 5(b)). In contrast, woody plant NPP attained greater relative dominance after 2045, offsetting all R h carbon losses by 2085; the woody plant NPP:R h ratio increased by a factor of ∼3 from 1980-2100.
The modeled change in tundra woodiness ( figure  5(b)) is an important ecosystem response to climate change that will have large impacts on the tundra carbon cycle. Consistent with this model result, Sistla et al (2013) reported that a two decades warming experiment in Arctic tundra ecosystem in Alaska resulted in increased woodiness and net ecosystem carbon storage. In another warming experiment in Toolik Lake, a low Arctic site in Alaska, Leffler et al (2016) reported increased shrub growth that enhanced ecosystem net carbon uptake. As an N limited ecosystem, the higher C:N ratio of woody versus non-woody plants is an important functional trait that enhances ecosystem carbon uptake from greater gains of carbon per N invested and longer turnover times. Woody carbon stocks with higher C:N ratios also decompose more slowly than non-woody plants (Weintraub and Schimel 2005), which sustains nutrient availability and slows carbon losses from heterotrophic decomposition. Cornelissen et al (2007) compared decomposition rates of leaf litter from species in 33 northern biomes in an incubation experiment and reported that herbaceous plant litter decomposes 40% faster than shrub litter, from which they concluded that woody shrub expansion could result in a negative feedback to global warming.
The effect of warming on NPP was shown to vary with PFT composition in several manipulative experiments , Cornelissen et al 1999, Elmendorf et al 2012, Hollister et al 2005. However, our predictions of changes in tundra woodiness and the resulting net carbon exchanges may be different than those inferred from warming experiments because ecosystem responses vary with warming experiment duration (Henry and Molau 1997, Hollister et al 2005. Most warming experiments have been relatively short-term, implying that relatively slower tundra PFT changes that occur over decades (e.g. from progressive warming over 21st century) may not have been captured (Rustad et al 2001).

Annual and seasonal trends in tundra net carbon exchange
Between 1982 and 2100, annual NPP increased by 244 g C m −2 while R h increased by 139 g C m −2 , resulting in the tundra ecosystem becoming a greater net carbon sink (106 g C m −2 ) (table 1). The magnitude and trend of tundra carbon dynamics were modeled to vary across seasons (table 1). Despite GPP increases in all seasons, the relative gains in carbon uptake were greater during spring and summer from enhanced carbon fixation rates and extended growing seasons (figure 3). The spring and summer uptake were modeled to offset the annual R e carbon losses, resulting in a net increase of the tundra carbon sink over recent decades and through the 21st century (table 1; figure  3). However, smaller increases in carbon uptake during autumn (22 g C m −2 ) and greater increases in R e (winter = 23 g C m −2 , Autumn = 47 g C m −2 ) resulted in larger carbon losses during these seasons by 2100 versus 1982. Increases in total net ecosystem carbon lost during autumn and winter (48 g C m −2 ) offset the net carbon gains in spring (13 g C m −2 ) and 25% of the gains in summer (141 g C m −2 ). These predictions of greater increases in carbon loss during winter and autumn with permafrost thaw suggest that the increase in North America Arctic tundra net carbon sink may not persist much past 2100. The modeled net carbon loses in autumn was consistent with Piao et al (2008), who used long-term atmospheric CO 2 concentrations and EC flux data to estimate greater autumn carbon losses versus gains, offsetting 90% of the carbon gains in spring in northern ecosystems.
The modeled increase in carbon loss during nongrowing seasons (table 1) is an important ecosystem process that may control the future tundra carbon cycle. In a warming experiment in the northern foothills of the Alaska Range, Natali et al (2012) reported that carbon loss from enhanced respiration during winter completely offset gains in net ecosystem carbon uptake during the growing season. Although this warming experiment only represents selected sites and could Table 1. Changes in annual and seasonal modeled carbon fluxes (g C m −2 ) of the North America Arctic tundra over the 21st century (average (1982-1986) subtracted from average (2096-2100) not capture the impacts of decadal-scale changes in plant communities, the implied importance of nongrowing season carbon losses on the tundra carbon cycle was consistent with our modeled result. Carbon losses from winter soil warming was also shown in a long-term snow fence warming experiment (Walker et al 1999) in Toolik Lake, Alaska. Increases in snow depth at the Toolik Lake site were also shown to enhance CH 4 production from soil warming that increased soil wetness that reduced soil O 2 concentration and thaw depth (Blanc-Betes et al 2016). Increases in graminoids at the site were shown to enhance transport of CH 4 fluxes out of the soil, implying the importance of PFTs in carbon flux transport mechanisms in Arctic moist tussock tundra. Changes in tundra net carbon exchange may also depend on the relative changes in air and soil temperatures, which affects changes in modeled NPP and R h (figure 3) by 2100. In ecosys, increasing T s causes greater ALD and soil aeration and hence O 2 uptake by microbes, thereby increasing R h (Grant et al 2015). Despite substantial increases in air temperature by 2100 (supplement II, figure 2), which resulted in GPP gains (table 1) from enhanced CO 2 fixation, soil temperature increased at a much slower rate (figure 6). Soil temperature increased in all seasons, albeit with a slower rate during the growing season versus winter. The slower increase in soil versus air temperature during the growing season is attributed to increased GPP, thus leaf expansion that increased canopy shading of the soil surface. Increased insulation from the soil surface litter layer during the non-growing season resulted in higher soil versus air temperature (figure 6). Overall, spatially averaged 0-15 cm soil temperature increased by ∼0.6 • C for every 1 • C increase in air temperature over the 21st century. This slower soil versus air warming, and thus greater effects on CO 2 fixation rate versus R h , contributed to the tundra remaining a carbon sink by 2100. However, these higher gains versus losses of carbon may not be sustainable under further soil warming beyond 2100.
The modeled results were affected by increases in atmospheric CO 2 , precipitation, N deposition, and N 2 fixation. Other processes beyond those considered here need to be included in land modeling assessments of future tundra carbon cycling. For example, shrub expansion may increase snowpack depth locally, which can insulate the soil and increase winter soil temperatures , although effects on shrub productivity vary (Myers-Smith and Hik 2013). Warming and associated atmospheric feedbacks from albedo changes (Bonfils et al 2012, Chapin et al 2005 and earlier snowmelt could further warm the soil, thereby amplifying carbon losses (Cahoon et al 2012) from R h . Fine-scale (0.1-1 km) spatial heterogeneity can affect vegetation dynamics (Jorgenson et al 2013); our ∼25 km resolution simulations did not resolve those spatial scales. Mechanisms of surface and sub-surface flows of water, nutrients, and energy from differences in topography can affect the simulation of the carbon cycle in Arctic ecosystems (Grant et al 2017a). Arctic CO 2 and CH 4 fluxes associated with formation of thaw lakes (Sturtevant and Oechel 2013), are important processes that affect modeling Arctic ecosystems. Finally, projected increases in tundra fires can amplify carbon losses (Mack et al 2011), and much work remains to develop accurate fire models for tundra systems.

Conclusions
Extended growing seasons and enhanced N mineralization and plant uptake were modeled to increase carbon fixation rates in recent decades and over the 21st century across North America tundra. Tundra shrub expansion and woodiness increased carbon uptake and retention sufficiently to offset the large concurrent increases in heterotrophic respiration, resulting in a predicted net North American tundra carbon sink by 2100. However, by 2100, increases in autumn and winter heterotrophic respiration were increasing rapidly and offset the spring and 25% of summer carbon uptake. Slower soil versus air warming resulted in greater increases in CO 2 fixation rates versus respiration, which may be unsustainable under further soil warming beyond 2100.