Warm-season net CO2 uptake outweighs cold-season emissions over Alaskan North Slope tundra under current and RCP8.5 climate

Arctic warming has increased vegetation growth and soil respiration during recent decades. The rate of Arctic warming will likely amplify over the 21st century. Previous studies have revealed that the most severe Arctic warming occurred during the cold season (September to May). The cold-season warming has posited significant CO2 emissions to the atmosphere via respiration, possibly offsetting warm-season (June to August) net CO2 uptake. However, prevailing Earth system land models poorly represent cold-season CO2 emissions, making estimates of Arctic tundra annual CO2 budgets highly uncertain. Here, we demonstrate that an improved version of the energy exascale Earth system model (E3SM) land model (ELMv1-ECA) captures the large amount of cold-season CO2 emissions over Alaskan Arctic tundra as reported by two independent, observationally-constrained datasets. We found that the recent seven-decades warming trend of cold-season soil temperature is three times that of the warm-season. The climate sensitivity of warm-season net CO2 uptake, however, is threefold higher than for the cold-season net CO2 loss, mainly due to stronger plant resilience than microbial resilience to hydroclimatic extremes. Consequently, the modeled warm-season net CO2 uptake has a larger positive trend (0.74 ± 0.14 gC m−2 yr−1) than that of cold-season CO2 emissions (0.64 ± 0.11 gC m−2 yr−1) from 1950 to 2017, supported by enhanced plant nutrient uptake and increased light- and water-use efficiency. With continued warming and elevated CO2 concentrations under the representative concentration pathway (RCP) 8.5 scenario, the increasing rate of warm-season net CO2 uptake is more than twice the rate of cold-season emissions (1.33 ± 0.32 gC m−2 yr−1 vs 0.50 ± 0.12 gC m−2 yr−1), making the modeled Alaskan Arctic tundra ecosystem a net CO2 sink by 2100. However, other geomorphological and ecological disturbances (e.g. abrupt permafrost thaw, thermokarst development, landscape-scale hydrological changes, wildfire, and insects) that are not considered here might alter our conclusion.


Introduction
Permafrost regions have undergone persistent warming during recent decades (Pithan and Mauritsen 2014, Huang et al 2017, Biskaborn et al 2019, and the Arctic tundra ecosystem has warmed more than any other biome (Bjorkman et al 2018). This warming, and likely the CO 2 fertilization effects, has increased ecosystem photosynthesis, productivity, and widespread expansion of tall shrubs in the Arctic tundra, leading to enhanced growing-season carbon uptake (Elmendorf et al 2012, Frost et al 2013, Zhang et al 2013, Martin et al 2017, Mekonnen et al 2018a, Berner et al 2020, Wang et al 2020. Meanwhile, permafrost thaw has made the organic carbon and nitrogen previously locked in permafrost soils available to microbial decomposition (Lawrence et al 2008, Koven et al 2011, Hugelius et al 2014, Mishra and Riley 2014, Schaefer et al 2014, thereby increasing both plant growth and ecosystem respiration (Waelbroeck et al 1997, McGuire et al 2012, Trucco et al 2012, Koven et al 2015, Schuur et al 2015, Parazoo et al 2018, Mekonnen et al 2018a, Gagnon et al 2019. The most severe Arctic warming has occurred during the cold season (Sturm et al 2005, Koenigk et al 2013, Cohen et al 2014, 2018, Huang et al 2017, Box et al 2019, which increases microbial decomposition of soil organic matter and thus enhances cold-season soil heterotrophic respiration (HR), releasing a significant amount of carbon via methane (CH 4 ) and carbon dioxide (CO 2 ) (Zona et al 2016, Natali et al 2019. Cold-season CO 2 emissions might largely offset the warm-season CO 2 net uptake, and become increasingly important in annual carbon budgets over Arctic tundra ecosystems (Commane et al 2017b, Jeong et al 2018, Campbell and Laudon 2019, Natali et al 2019.
Site-scale measurements have demonstrated large cold-season CO 2 losses over the pan-Arctic tundra ecosystem (e.g. Fahnestock et al 1998, Jones et al 1999, Kittler et al 2017. Some studies indicated that with warming, cold-season CO 2 emissions from the Alaskan Arctic tundra grew larger than summer net uptake, i.e. shifted from an annual sink to an annual source of CO 2 (Oechel et al 1993(Oechel et al , 2000(Oechel et al , 2014. Belshe et al (2013) also concluded that tundra sites were annual sources of CO 2 from the mid-1980s until the 2000s. Despite these conclusions at the site scale, long-term regional estimates of coldseason CO 2 emissions remain uncertain due to limited spatial representativeness of site-scale observations (Rustad et al 2001, McGuire et al 2012 and the lack of continuous measurements throughout the entire cold season (Oechel et al 1993, McGuire et al 2012, Parazoo et al 2016. Process-based terrestrial ecosystem models can facilitate the estimate of long-term, regional coldseason CO 2 budgets. However, most current land models integrated within Earth system models (ESMs) poorly represent cold-season carbon emissions (Zona et al 2016, Commane et al 2017b, Wang et al 2019, Natali et al 2019, making estimations and predictions of regional net ecosystem exchanges (NEE) of CO 2 over the Arctic tundra highly uncertain (McGuire et al 2012, Fisher et al 2014. The coupled model intercomparison project 5 (CMIP5) models, for instance, substantially underestimated early winter respiration rates over Alaska (Commane et al 2017b), partially due to inadequately represented zero-curtain periods (ZCP) (i.e. the period when soil temperatures linger around 0 • C during the freezing season) (Outcalt et al 1990). With a warming climate, the ZCP duration over the Arctic tundra is expected to increase (Zona et al 2016, Arndt et al 2019. However, different representations of ZCP and cold-season soil respirations over Arctic tundra might lead to contradictory conclusions on the annual CO 2 budget among models. For instance, Zhang et al (2014) predicted that under the representative concentration pathway (RCP) 8.5 scenario, the Arctic tundra will remain a carbon sink throughout the 21st century due to enhanced biogeophysical feedbacks. Similar conclusions were found by many other studies (e.g. Qian et al 2010, Zhang et al 2013, Mekonnen et al 2018a. Other modeling studies, in contrast, have concluded that cold-season carbon emissions might eventually offset summer net carbon uptake under 21st century warming climate, potentially transforming permafrost tundra ecosystems from a net sink to a net source of carbon (e.g. Piao et al 2008). Without close examination of the representation of ZCP and cold-season carbon emissions, it is hard to quantify the annual CO 2 budget with high confidence.
Relying on site-scale observations, Natali et al (2019) derived winter-time CO 2 emissions over highlatitude permafrost regions with a machine learning approach, and reported that 1662 TgC yr −1 were lost via CO 2 from their domain from October to April between 2003 and 2017. By combining the machine-learning cold-season CO 2 emission estimates and independent multiple-model estimates of warm-season net CO 2 uptake, they further concluded that the current permafrost region is a net annual CO 2 source with annual emissions ranging from 15 to 975 TgC yr −1 , and that, under a continued warming climate, increases in cold-season CO 2 emissions might exceed increases in warm-season net uptake. However, large uncertainties exist in the predicted future emissions because the data-driven model trained by current in-situ measurements might not be able to capture future ecosystem responses to a changing climate (Natali et al 2019). Therefore, reasonably representing the ZCP and overall cold-season carbon processes within process-based land models is urgently needed to better predict carbon budgets of Arctic tundra ecosystems (Commane et al 2017b, Campbell andLaudon 2019).
Recently, relying on observations from the Arcticboreal vulnerability experiment (ABoVE) and the carbon in Arctic reservoirs vulnerability experiment (CARVE), Tao et al (2020) improved the energy exascale Earth system model (E3SM) land model (ELMv1-ECA) in terms of simulating ZCPs and coldseason CO 2 emissions at several Alaska tundra sites (figure 1). Here we adopted the updated ELMv1-ECA to address the long-standing question of whether cold-season CO 2 emissions offset the warm-season net uptake over the Alaskan North Slope tundra (NST) currently and over the 21st century.

Study domain and data
We mainly focused on the Alaskan NST (figure 1), adopting the landcover map in Commane et al (2017b). We also examined results over the interior boreal forest (BF) and Southwest tundra (SWT) areas in supplementary section sup. 2 (available online at stacks.iop.org/ERL/16/055012/mmedia). We used two independent observationally-constrained spatial datasets to evaluate model results. Details about the two datasets and processing methodology are provided in the supplementary section sup. 1.
First, we used a spatial NEE dataset optimized with CARVE aircraft measurements (Commane et al 2017a(Commane et al , 2017b as a benchmark to evaluate the regional simulation results over Alaska from 2012 to 2014. The CARVE airborne measurements span April through November, and thus the optimized NEE are only available over these months (supplementary figure S1) 3 . Hereafter, we use 'C2017' to represent this observation-constrained, aircraft-optimized spatial NEE dataset.
We also used another observation-constrained spatial dataset of cold-season CO 2 emissions derived from ground flux observations with a boosted regression tree machine learning method (Watts et al 2019, Natali et al 2019. We denote this observationconstrained dataset as 'N2019' . The N2019 includes spatially-resolved September to April CO 2 emissions over northern permafrost regions from 2003 to 2017, and predicted cold-season CO 2 emissions from 2018 to 2100 using their method and CMIP5 outputs under the RCP8.5 scenario.

Model and experiment design
The ELMv1-ECA simulates terrestrial carbon cycles, energy, and water exchange fluxes, and accounts for the limitation of nutrient availability for plant growth (Riley et al 2018, Chen et al 2019, Golaz et al 2019. Here, we used the version of ELMv1-ECA updated by Tao et al (2020), denoted as ELMv1a hereafter. ELMv1a differs from ELMv1-ECA mainly in its phase-change scheme and soil organic carbon decomposition scheme, and coldseason methane transport mechanism. We have calibrated and evaluated the ELMv1a against eddy covariance observations at Alaska tundra sites (figure 1) (Tao et al 2020). Compared to ELMv1-ECA, ELMv1a has shown enhanced performance in simulated ZCPs of active-layer soils and cold-season CO 2 and emissions at the site scale. At the regional scale over the NST, we identified a generic decomposition scheme which uses a Q 10 value of 2.0 and a modified moisture-dependency function that prevents zero respiration in freezing and subfreezing soils (Tao et al 2020).
With ELMv1a, we conducted a transient simulation from 1901 to 2017 at 30 min temporal resolution over Alaska, driven by 0.5 • × 0.5 • climatic research unit Japanese reanalysis (CRU JRA) climate forcing (Harris 2019). The transient simulation started from properly spun-up simulations following Zhu et al (2019). We estimated model uncertainty associated with representations of temperature sensitivity to HR through simulations with a range of Q 10 values (1.6-2.2) (Tao et al 2020). We then ran three future simulations through 2100 under the RCP8.5 scenario (Friedlingstein et al 2006), including (a) with RCP8.5 climate forcing but CO 2 concentration fixed at 2011 levels (denoted as RCP8.5_climate), (b) with evolving elevated CO 2 concentration but repeatedly cycling 2011-2020 climate forcing (denoted as RCP8.5_eCO 2 ), and (c) with both predicted RCP8.5 climate and evolving CO 2 concentration (denoted as RCP8.5_climate + eCO 2 ). Then, we characterize the trends of warm-season net CO 2 uptake and coldseason CO 2 loss over the Alaskan NST and the whole Alaska domain.
Studies have shown that net primary productivity (NPP) enhancement resulting from elevated CO 2 might be attenuated by limited nutrient availability (Norby et al 2010, Zaehle et al 2014, Wang et al 2020. We thus disentangled the integrated ecosystem response into major nitrogen (N) and phosphorus (P) cycling processes, including plant N and P uptake, net N and P mineralization, N and P limitation, N-and P-use efficiency (NUE and PUE). We also examined water-and light-use efficiency (WUE and LUE). Details about these subcomponent processes are provided in the supplementary sup. 4.
We defined 'warm season' as June to August and 'cold season' as September to May. We further defined an early cold season period from September through October. We estimated the NEE trend for each period as the regression slope between annual or seasonal NEE and time. The computed trend is statistically significant if the p-value is less than 0.05. Further, similar to Ballantyne et al (2017), we calculated and discussed the climate sensitivity of warm-season total net CO 2 uptake and cold-season total net CO 2 emissions (supplementary sup. 3).

Evaluation of simulated NEE
We evaluated simulated NEE against the two spatial datasets (i.e. C2017 and N2019) at both pixel and regional scales. We used C2017 to evaluate results during the warm season and at the annual scale from 2012 to 2014. We also compared results with both N2019 and C2017 over their common period (i.e. 2012-2014). Then, we compared simulated coldseason CO 2 emissions against N2019 throughout 2003-2017.
Pixel-scale RMSE of ELMv1a-simulated monthly NEE (figure 2(b)) varied between 0.3 and 3.7 gC m −2 d −1 for the warm season and between 0.2 and 1.4 gC m −2 d −1 for the cold season (compared to C2017). There were large spatial patterns in the NEE RMSE, with the largest cold-season RMSE values clustering in the interior BF regions and largest warmseason values scattering across Alaska ( figure 2(b)). Reasons potentially contributing to higher NEE RMSE include inaccurate forcing and landscape parameters prescribed and model deficiencies, as summarized at the end of this section. Compared to ELMv1-ECA (figure 2(a)), ELMv1a improved simulated cold-season CO 2 emissions at the pixel scale, especially over the NST. The RMSE (compared with C2017) of ELMv1a-simulated regional mean NEE over the NST was reduced from 0.43 gC m −2 d −1 to 0.39 gC m −2 d −1 (figure S2). Evaluating against the N2019 CO 2 emissions over September to April, ELMv1a showed a significantly improved performance over ELMv1-ECA, with a 55% reduction in RMSE (0.14 vs 0.31 gC m −2 d −1 ) (figure S2). We, therefore, apply ELMv1a-simulated results for the remaining analyses.
ELMv1a overestimated warm-season CO 2 net uptake compared to C2017 (figure 3) by 0.23 gC m −2 d −1 (31%) and 0.44 gC m −2 d −1 (44%) on average for NST and the whole Alaska domain, respectively. Shoulder-season (i.e. September, October, and May) CO 2 release was underestimated compared to C2017 but not compared to N2019. The ELMv1a-simulated regional means agreed well with C2017 over the rest of the cold-season months for the NST and SWT regions, with slight overestimations for cold-season emissions over the BF region (figure 3). Compared to N2019, results demonstrate very good performance for all the major landcover types and entire Alaska. Because of a lack of coldseason aircraft measurements (figure 3) for C2017, we relied on N2019 for this comparison. Specifically, results show the best performance for the regional mean over entire Alaska with the smallest RMSE of 0.08 gC m −2 d −1 , followed by the RMSE of 0.14 gC m −2 d −1 for tundra ecosystems (table S1).
The ELMv1a simulated 2012-2014 averaged cold-season total CO 2 emissions over the NST show biases within 30% of the observed values, i.e. 24.7 gC m −2 yr −1 (24%) compared to N2019 and 31.9 gC m −2 yr −1 (29%) compared to C2017 (figure 4(a)). For the early cold season (September to October), when the ZCP of the active layer in Arctic tundra is most often present, the model simulated total CO 2 emissions match very well with N2019, with biases of 8.2 gC m −2 yr −1 (21%) for NST and 0.9 gC m −2 yr −1 (2%) for Alaska. The annual total NEE (positive) from C2017 indicates that the NST ecosystem is a CO 2 source on average from 2012 to 2014 ( figure 4(a)). The ELMv1a also estimates an annual CO 2 source for this tundra area, although it is lower by 53.5 gC m −2 yr −1 compared to C2017, mainly due to higher modeled summer CO 2 uptake ( figure 3(a)). The annual total difference is smaller (17.1 gC m −2 yr −1 ) for entire Alaska ( figure 4(b)), offsetting the larger summer CO 2 uptake by the slightly larger CO 2 emissions from the BF area (figure 3(c)).
Despite the larger summer uptake and thus lower annual emissions, ELMv1a better matched the C2017 and N2019 estimates compared to CMIP5 models reported by Commane et al (2017b). Specifically, compared to N2019, ELMv1a displays smaller differences in September to December total CO 2 emissions (figure 4(c)), i.e. 0.49 TgC yr −1 and 3.28 TgC yr −1 for ELMv1a and MIROC-ESM-CHEM (i.e. the best CMIP5 model compared to C2017), respectively. Indeed, MIROC-ESM-CHEM predicted a 1 month advanced beginning of the growing season compared to C2017 (Commane et al 2017b). In comparison, ELMv1a generally preserves reasonable estimations of peak summer NEE timing.
Compared to N2019 from 2003 to 2017, ELMv1a simulated slightly smaller total CO 2 emissions during the early cold season (September and October) (figure 5), i.e. −9.1 gC m −2 and −1.2 gC m −2 on average for NST and Alaska, respectively. The model predicted lower total CO 2 emissions during the cold season with a mean difference of −31.5 gC m −2 for the NST, but showed very good agreement for the regional mean over Alaska with a small mean difference of −7.9 gC m −2 . The modeled annual CO 2 budget for the NST indicates a net sink of −7.3 gC m −2 averagely from 2003 to 2017. Considering the underestimated cold-season total CO 2 emissions compared to N2019, however, would alter our conclusion from an annual sink to an annual source (+26.9 gC m −2 ). Similarly, the modeled SWT is a net sink of −15.7 gC m −2 and BF is a net source 8.2 gC m −2 , which would be changed to a net source (3.0 gC m −2 ) and a net sink (−19.6 gC m −2 ) for SWT and BF, respectively when replacing cold-season estimates with N2019 (sup. 2). Nonetheless, the spatial variability of CO 2 emissions across Alaska within N2019 is surprisingly small (figure S1). The large uncertainty in N2019 (about 0.93 gC m −2 d −1 for the NS tundra; personal communication with Dr Jennifer Watts) makes the annual CO 2 budget estimate highly elusive.
Overall, despite showing better performance than other CMIP5 model projections, ELMv1a predicts larger summer CO 2 uptake than C2017 and has smaller cold-season CO 2 emissions compared to C2017 and N2019 over the NST. Possible reasons for the  Watts) for the NST area, and thus is not shown here. Bluish shaded area indicates the model uncertainty associated with representations of temperature sensitivity to HR (section 2.2). Gray background highlight periods when CARVE aircraft measurements were available and used for optimization by C2017. Note, since few aircraft measurements were available in the cold season for optimization (supplementary figure S1), the C2017 estimates of NEE during cold-season months solely come from simulated results (supplementary section sup. 1). discrepancies between ELMv1a simulations and the observationally-constrained datasets include inaccurate driver and landscape parameters and model deficiencies. Specifically, unresolved grid-cell landscape heterogeneities (particularly in topography, vegetation, and soil properties) impact model estimates of carbon budget via influencing air-vegetationground water, carbon, and energy exchanges. The spatial resolution of the reanalysis forcing we used may also be coarser than exists in heterogeneous landscapes. Relevant model deficiencies include: (a) ELMv1a parameterizations for pan-Arctic plant functional types are known to need improvement (Sulman et al 2021); (b) ELMv1a lacks representations for dynamic vegetation (e.g. shrub expansion) and ecological disturbances from insects; (c) ELMv1a lacks an appropriate representation for geomorphological disturbances to permafrost landscapes caused by abrupt permafrost thaw and thermokarst formation that release a large amount of permafrost carbon via greenhouse gases (Nitzbon et al 2020, Turetsky et al 2020; and (d) ELMv1a simulations here did not incorporate the substantial ecosystem disturbance caused by wildfires which also initiate widespread abrupt permafrost thaw and thermokarst development (Holloway et al 2020).

Historical and predicted trends in NEE
The spatial pattern of trends in soil temperature reveals that the largest warming in the top 50 cm soil has occurred in the cold season over the NST (figures 6(a) and (d)). The warming trend of coldseason soil temperature is three times that of the warm-season soil temperature over the NST, i.e.    of regionally averaged total NEE over (a) the Alaskan NST area and (b) the whole Alaska domain. The total NEE during the early and entire cold season from N2019, an observation-based CO2 emission dataset derived with a boosted regression tree machine learning method, are also shown. Estimated negative trends of warm-season total NEE (gC m −2 yr −1 ) represent increasing net CO2 uptake, which outweigh the increasing rates of total cold season CO2 emissions. Shaded areas indicate the uncertainty bounds (i.e. 95% confidence interval) associated with the estimated trends. Figure 6. Spatial map of warm season (June to August) and cold season (September to May) trends  in mean soil temperature, liquid VSM, and total net CO2 uptake or emissions over Alaska. Pixels with a significant trend (p value less than 0.05) are marked by black circles. Both warm-season net CO2 uptake and cold season net CO2 emissions show widespread increasing trends. Warming trends in cold-season soil temperature generally are larger than the trends in warm-season soil temperature, especially over the North Slope. 0.06 • C yr −1 vs 0.02 • C yr −1 (table 1). The map of warm-season trends in top 50 cm liquid volumetric soil moisture (VSM) shows fewer pixels with statistically significant trends (indicated by black circles) than the cold-season map over the NST (figures 6(b) and (e)). The warm-season regional mean VSM over the NST shows the same (but statistically nonsignificant) increasing rate as for the cold-season (0.04% yr −1 ; table 1). Both warm-season net ecosystem CO 2 uptake and cold-season net CO 2 emissions show widespread increasing trends over the NST (figures 6(c) and (f)), but showing very different spatial patterns.
Integrated over the NST, we found that warmseason net CO 2 uptake has increased faster than that of cold-season net CO 2 emissions between 1950 and 2017 (figure S6(d) and table 1). Due to recent accelerated warming, the increases in warm-season NEE between 2000 and 2017 were greater than those for the long-term historical periods (i.e. 1950-2017 and 1980-2017) (figure S6(d)), supported by larger increases in warm-season monthly NPP (figure S6(e)) and enhanced plant N and P uptake ( figure  S7). The increases in HR between 1950 and 2017 were apparently smaller than those between 1980 and 2017 ( figure S6(c)), which is attributed to consecutive hydroclimatic disturbances from the late 1960s to early 1970s (Bieniek and Walsh 2017, Sulikowska et al 2019), i.e. dry extremes (late 1960s) followed by cold winter extremes (early 1970s) (figure S8). Due to soil moisture and temperature memories of hydroclimatic variabilities, microbial decomposition did not recover to the rate before the hydroclimatic disturbance until several years later (figures S8 and S9). Similar disturbances were also shown for net N and P mineralization (figure S7). In contrast, NPP (figure S6(b)), plant N and P uptake (figures S7 and S8), and leaf C, N and P (figure S10) did not show large reductions during these disturbance periods mainly due to plant nutrient storage and flexible plant stoichiometry for biomass growth (Riley et al 2018), demonstrating strong modeled plant resilience to these types of hydroclimatic extremes. In addition, the NST ecosystem sustained productivity under the dry summer extremes (late 1960s) also because of increased WUE (figure S11). However, due to the strong relationship between warm-season NPP and soil temperature (figure S9(e)), cold warm-season extremes would exert a larger impact on plant productivity than HR for a short-term period (e.g. 1998-2000 in figures S8(c.1) and (e.1)).
The climate sensitivity of warm-season net CO 2 uptake, i.e. changes in warm-season cumulative net CO 2 uptake in response to changes in warmseason soil temperature, is 25.1 gC m −2 • C −1 (figure 7(a)), which is three times greater than the climate sensitivity of cold-season CO 2 emissions (6.8 gC m −2 • C −1 ) ( figure 7(b)). The large Table 1. Historical trend (1950Historical trend ( -2017 in regional CO2 emission, 2 m air temperature, precipitation accumulation, soil temperature and liquid water content in top 50 cm from different periods over the Alaskan NST and Alaska. The negative trend in warm-season total NEE indicates increasing net CO2 uptake. The uncertainty bounds (i.e. 95% confidence intervals) associated with the estimated trends are indicated in parentheses.   . The change of last 10 yr mean relative to the first 10 yr mean for multiple variables including NPP, N-and P-based NPP (NPPN and NPPP), leaf N and P, plant N and P uptake, net N and P mineralization, N and P limitation factor, N-and P-use efficiency (NUE and PUE), light-and water-use efficiency (LUE and WUE) for (a) the transient run and (b) three RCP8.5 runs (see details in section 2.2) over the NST. climate sensitivity of warm-season net CO 2 uptake was reported earlier with tower observations by Yi et al (2010). It is mainly caused by enhanced plant productivity due to increased plant N and P uptake and increased LUE and WUE (figures S7-S11). As a result, despite the fact that cold-season soil temperature has warmed three-times as fast as warm-season temperature (table 1), the trend in warm-season net CO 2 uptake is expected to exceed the increasing coldseason CO 2 emissions due to the larger emergent climate sensitivity of warm-season net CO 2 uptake. Indeed, warm-season net CO 2 uptake had a larger increasing trend (0.74 ± 0.14 gC m −2 yr −1 ) than the positive trend for the cold-season total CO 2 emissions (0.64 ± 0.11 gC m −2 yr −1 ) ( figure 5 and table 1), leading to a statistically non-significant annual budget trend for the NST between 1950 and 2017 (table 1). The warming climate has increased the first 30 yr  cold-season CO 2 emissions from the NST by 44% (i.e. 86.4 gC m −2 vs 60.0 gC m −2 ) relative to the last 30 yr average , and increased the warm-season net CO 2 uptake by 47% (i.e. 92.3 gC m −2 vs 62.7 gC m −2 ).
With respect to future predictions, the ELMv1a results (RCP8.5_climate + eCO 2 ) of early cold season CO 2 emissions are similar to those predicted by N2019 for NST (with a mean difference as −10.9 gC m −2 ) and Alaska (with a mean difference as 13.3 gC m −2 ). Under this future scenario, both warm-season net CO 2 uptake and cold-season CO 2 emissions are projected to continue increasing (figure 8). Compared to the first 30 yr average , the NST cold-season total CO 2 emissions averaged over the last 30 yr (2071-2100) increased by 30% (i.e. 118.4 gC m −2 vs 90.8 gC m −2 ), while warm-season total net CO 2 uptake increased by 108% (i.e. 133.4 gC m −2 vs 64.2 gC m −2 ). Indeed, the increasing rate of warmseason net CO 2 uptake is more than twice that of cold-season CO 2 emissions from 2018 to 2100, i.e. 1.33 ± 0.32 gC m −2 yr −1 vs 0.50 ± 0.12 gC m −2 yr −1 , making this critical tundra ecosystem likely a CO 2 sink by 2100, with a significant increasing sink trend at 0.87 ± 0.34 gC m −2 yr −1 . However, the model predicted smaller cold-season CO 2 emissions than N2019 (figure 9) probably because it lacks representations for abrupt permafrost thaw and thermokarst dynamics; and it might also overestimate CO 2 uptake since it does not account for ecosystem disturbances by insects and wildfires. Also, although ELMv1a can simulate dynamics in plant growth and structure (e.g. leaf, stem, roots biomass) via simulating belowand above-ground biomass growth in response to light, water, and nutrient availability, currently it does not represent tundra plant community shifts (e.g. shrub expansion) which might increase CO 2 uptake (Mekonnen et al 2018a, 2018b). Representing these disturbances and vegetation community shifts is urgently needed within ESMs to better estimate current and future carbon budgets.
Results showed that enhanced plant N and P uptake and net N and P mineralization support the large CO 2 uptake during the warm season for both the transient (figures 9(a) and S11) and RCP8.5 runs (RCP8.5_climate and RCP8.5_climate + eCO 2 ) (figures 9 and S12). Although showing positive relative changes in plant productivity for the transient baseline run (figure 9(a)), results indicate potential declines in CO 2 fertilization effects on plant photosynthesis (figure S11) possibly due to nutrient limitations (as also shown by Wang et al (2020) for the NST). We evaluated this effect by comparing factorial experiments under the RCP8.5 scenario (section 2.2). Specifically, with the climate fixed at 2011-2020 conditions (RCP8.5_eCO 2 ), elevated CO 2 first causes large NPP enhancement (figure S13), but then increases in N and P limitation to plant growth from microbial N and P immobilization and decreases in net N and P mineralization (green in figures 9(b) and S13, and the N(P) limitation factor in figure S12). Consequently, the NPP enhancement due to elevated CO 2 is gradually attenuated due to limited N and P availability (RCP8.5_eCO 2 in figure S13). In contrast, the continued warming deepens active layer thickness (figure S14) and thus increases N and P availability, thereby somewhat relieving the N and P limitation, causing negative changes in the N limitation factor (red bar, figure 9(b)) and smaller increases in P limitation (red vs green bars, figure 9(b)). Indeed, both N and P limitation factors first increase and then decrease after their peaks with continued warming (RCP8.5_climate and RCP8.5_climate; red and yellow, figure S12), but persistently increases without climate warming (RCP8.5_eCO 2 ; green, figure S12). In addition, increased LUE and WUE (figure 9) also contribute to the large trend of warm-season CO 2 net uptake.

Summary and conclusions
We used an improved version of the E3SM land model (Tao et al 2020: ELMv1a) to investigate warm-season CO 2 net uptake vs cold-season CO 2 net emissions over Alaska under historical and the 21st century RCP8.5 scenario. The model results agree reasonably well with two independent observationally-constrained datasets derived with different methods: using atmospheric measurements and inversion (C2017) and using surface-flux measurements and machine learning (N2019). Results demonstrate that over the Alaskan NST, although cold-season soil warming was greater than that during the warm-season, the warm-season net CO 2 uptake increased faster (0.74 ± 0.14 gC m −2 yr −1 ) than that of cold-season CO 2 release (0.64 ± 0.11 gC m −2 yr −1 ) between 1950 and 2017. These patterns are held as well for future projections under the RCP8.5 scenario.
The response is attributed to a threefold larger climate sensitivity of warm-season net CO 2 uptake than that of cold-season CO 2 losses. This difference in climate sensitivity is partially (a) due to strong plant resilience to hydroclimatic disturbances, enhanced plant nutrient uptake, increased water-and light-use efficiency, and partially (b) due to weak microbial resilience, i.e. the slow recovery of microbial decomposition rates following consecutive dry summer and cold winter extremes that were built in soil moisture and temperature memories for extended periods. Based on ELMv1a simulations, we therefore conclude that the NST ecosystem will likely remain an annual CO 2 sink through the 21st century. We acknowledge that several factors, including fire, insects, abrupt permafrost thaw, thermokarst dynamics, and landscapescale hydrological changes that are not discussed here, might alter our conclusion.

Data availability statement
The modeling results that support the findings of this study are available upon reasonable request from the authors.