The future of tundra carbon storage in Greenland-Sensitivity to climate and plant trait changes

The continuouschange in observed keyindicators such as increasing nitrogen deposition, temperatures andprecipita-tion will have marked but uncertain consequences for the ecosystem carbon (C) sink-source functioning of the Arctic. Here, weuse multiple in-situ data streamsmeasured bythe Greenland Ecosystem Monitoringprogramme in tight connection with the Soil-Plant-Atmosphere model and climate projections from the high-resolution HIRHAM5 regional model. We apply this modelling framework with focus on two climatically different tundra sites in Greenland (Zackenberg and Kobbefjord) to assess how sensitive the net C uptake will expectedly be under warmer and wetter conditions across the 21st century and pin down the relative contribution to the overall C sink strength from climate versus plant trait variability. Our results suggest that temperatures (5 – 7.7 °C), total precipitation (19 – 110 %) and vapour pressure de ﬁ cit will increase (32 – 36 %), while shortwave radiation will decline (6 – 9 %) at both sites by 2100 under the RCP8.5 scenario. Such a combined effect will, on average, intensify the net C uptake by 9 – 10 g C m − 2 year − 1 at both sites towards

Recent past and future 5-year annual mean net ecosystem exchange (NEE), gross primary production (GPP) and ecosystem respiration (R eco ) projected for Kobbefjord and Zackenberg using the SPA model forced with HIRHAM5 climate data.

Introduction
The Arctic is one of the regions on our planet with the strongest observed warming (IPCC, 2019), and is expected to warm strongly in the coming decades (AMAP, 2022).Increase in temperatures (Serreze and Barry, 2011), intensification of the hydrological cycle (Bintanja et al., 2020) and shortening of the spring snow cover duration (Derksen and Brown, 2012) will have marked consequences for ecosystem-atmosphere interactions, including both abiotic and biotic components at northern latitudes.Physical changes may trigger multiple ecosystem responses in carbon (C) sink-source functioning but also complex feedbacks that will further increase temperatures at high latitudes, a phenomenon called 'Arctic Amplification' (Previdi et al., 2021).Hence, it is important to address the gap in the understanding of Arctic C cycle sensitivity to climate change and so decrease the current large uncertainties (López-Blanco et al., 2019;Virkkala et al., 2021) in future projections.
The Arctic holds approximately half of the global soil organic carbon (SOC) (Hugelius et al., 2014;Jackson et al., 2017), and so plays an important role in the global C cycle.Yet, the terrestrial C cycle is currently the least constrained component of the global C budget (Bloom et al., 2016;Friedlingstein et al., 2006).In Arctic wetlands, particular attention has been paid to net ecosystem CO 2 exchange (NEE) as a key greenhouse gas and a fundamental part of biotic interactions with the atmosphere.The expected increase in temperature and precipitation (AMAP, 2022) may stimulate increased emissions and uptake in the Arctic C cycle.On the one hand, warming can increase CO 2 release via ecosystem respiration (R eco ) as a result of, e.g., thawed permafrost soils (Turetsky et al., 2020;Schädel et al., 2016), enhanced microbial turnover (Commane et al., 2017), enhanced heterotrophic respiration (Webb et al., 2016), extreme weather events (Christensen et al., 2020) or pest outbreaks (López-Blanco et al., 2017;Lund et al., 2017).On the other hand, there may be increases in gross primary production (GPP) caused by vegetation greening (Myers-Smith et al., 2020), shrub expansion (Myers-Smith et al., 2011), CO 2 fertilisation (Park et al., 2020), enhanced nutrient (N) availability (López-Blanco et al., 2020) or interactions with large herbivores (Mosbacher et al., 2019).These sensitive counteracting processes of uptake and loss, GPP and R eco , shape NEE and may feedback further with climate, thereby constantly rearranging this delicate balance, resulting in uncertain flux signs (C sink/source) and magnitudes.
The terrestrial C cycle of Greenland is generally understudied compared to other Arctic regions.Although in-situ C flux measuring techniques are widely used in ecosystem monitoring (Pastorello et al., 2020), their temporal and spatial coverage has been limited in Greenland.Logistic challenges in Greenland, mainly due to remoteness and harsh conditions, make it difficult to create long-term continuous datasets.During the last 25 years, a fundamental aim of the Greenland Ecosystem Monitoring programme (GEM; https://g-e-m.dk/) has been to ensure the continuity and integrity of its long-term time series of monitored variables and to improve knowledge of ecosystem processes and connections (Christensen et al., 2017).Additionally, as complements to in-situ long-term monitoring activities, ecosystem and regional climate models provide a unique opportunity 1) to understand in more detail these ecosystems and controlling processes beyond the temporal coverage of the available time series but also 2) to provide projections of ecosystem changes in response to different greenhouse gas emission scenarios (IPCC, 2019).Process-oriented models are powerful tools that can represent complex ecosystem processes determining the NEE of CO 2 controlled at the same time by meteorological forcing, vegetative phenology, snow dynamics, plant N traits and CN stocks from soils (Williams et al., 2000;López-Blanco et al., 2020).Moreover, coupled climate models are the best (and only) tools we have to explore and assess the changes in the recent past and future of the climate system (O'Neill et al., 2016) and the associated ecosystem responses (van den Hurk et al., 2016).Although large-scale climate models provide convincing numerical estimates at regional or global scales (Jones et al., 2016;Taylor et al., 2012), differences in model performance, emission scenarios and climate projections are far from perfection (Fisher et al., 2014;Nishina et al., 2015).Consequently, detailed comparisons with available in-situ data are critical to assess model biases and unavoidable uncertainties.
The aim of this paper is to quantify the relative sensitivity of Greenland's C balance to climate change based on regional variation in C cycling across a Greenland tundra climate and plant trait-nutrient gradient (Table 1).The key constraints to our understanding have been sparse time series of C fluxes and limited regional data.Now, with eddy flux and ecosystem trait observations over multiple years linked to ecosystem and climate models, we can project the variation of vegetative phenology, productivity and soil decomposition forward in time.We use proven model tools to analyse

Table 1
Meteorology and vegetation characterisation of the Kobbefjord and Zackenberg sites.Complementary pictures illustrating the tundra vegetation from both stations can be found in Fig. S1. a Note: The beginning and end of growing season are defined as the first and last day when there were three consecutive days with daily average net CO 2 uptake and release, respectively.The length of the growing season is the period between the beginning and the end of the growing season.
the underlying processes and links between the future climate, plant traits and terrestrial C cycling.Here, we apply a C cycle model, the Soil-Plant-Atmosphere model (Williams et al., 2000;Williams et al., 1996), to two GEM sites relying on previous substantiated efforts on source-code model implementation, model constraints, calibration and validation based on quality-controlled long-term data (López-Blanco et al., 2018;López-Blanco et al., 2020).In this study, our model framework is now forced forward in time with future projections from the high-resolution regional climate model HIRHAM5 developed by the Danish Meteorological Institute (Christensen et al., 2006;Boberg et al., 2018;Langen et al., 2015), which is specifically designed to characterise Greenland (and resolve its topography adequately) and follows different IPCC (2013) RCP greenhouse gas emission scenarios.With the complex geographical settings of the ice-free parts of the Greenlandic coastal margins, it is essential to use as high a horizontal resolution as possible.For that purpose, the HIRHAM5 model has here been operated at a grid resolution of 5.5 by 5.5 km, which depicts the main land sea contrasts and most coastal mountains adequately for a direct comparison with site-specific measurements.
We ask the ecological questions: "How sensitive is the C balance expected to be under warmer and wetter conditions projected for the 21 st century?" and "Is the C sink strength more sensitive to local characteristics such as plant nitrogen traits dependent on soil and geology conditions and other plant traits than to projected climate variability?"We hypothesise that the net C exchange will be exposed to higher temperatures and intensified precipitation levels increasing the future C sink strength, but plant nutrient traits will also play a significant role.The local nutrient status may be a more important regulator of the contemporary net C uptake than temperatures or growing season length (López-Blanco et al., 2020).Overall, we expect that higher photosynthetic rates will exceed respiratory losses following recent greening trends in response to warmer and longer snow-free periods (Myers-Smith et al., 2020;Berner et al., 2020).

Study sites
Multiple studies in Greenland have been conducted at the Zackenberg (Northeast Greenland, 74°28′N, 20°34′W) and the Nuuk-Kobbefjord fens (Southwest Greenland, 64°07′N; 51°21′W) since 2008 representing high and low Arctic tundra sites, respectively (Christensen et al., 2017).Extensive terrestrial CO 2 exchange measurements have been conducted at both sites under the auspices of the cross-disciplinary Greenland Ecosystem Monitoring (GEM) programme (https://g-e-m.dk/).The key environmental conditions and plant traits characterising Kobbefjord and Zackenberg are synthesised in Table 1 based on the site descriptions given in López-Blanco et al. (2020).

Data collection
2.2.1.In-situ CO 2 flux data and ancillary data NEE data consist of 30-minute temporal resolution measurements using eddy covariance (EC) (Baldocchi, 2003) during 2008-2016in Kobbefjord and 2008-2019in Zackenberg (López-Blanco et al., 2020).In Kobbefjord, the EC system consisted of a closed-path infrared gas analyser LI-7000 and a 3-D sonic anemometer Gill R3-50 at a height of 2.2 m (air intake at 2 m).In Zackenberg, the system was equipped with a closed-path infrared gas analyser LI-6262 and a 3-D sonic anemometer Gill R2 at a height of 3 m (and the air intake at the same level) until August 2012 when it was upgraded to a closed-path LI-7200 and Gill HS.Detailed information on EC system setup, flux pre-processing and quality checks of the systems can be found in López-Blanco et al. (2017).The quality-checked NEE data (with gaps) have been filled using the marginal distribution sampling method (Moffat et al., 2007), while the separation of NEE of CO 2 into its two key components GPP and R eco was done using established algorithms frequently utilised by the FLUXNET community (Pastorello et al., 2020).The data have been processed using 1) the night-time method (Reichstein et al., 2005) for Kobbefjord and 2) the daytime method (Lasslop et al., 2010) for Zackenberg due to the absence of true night-time during summer.The resulting fluxes follow the standard micrometeorological sign convention, i.e., C uptake is negative and C release is positive fluxes.
Additionally, a comprehensive suite of meteorological measurements was collected, processed and quality checked to complement the C flux data.The ancillary data, available at Zackenberg from 1995 to 2020 and at Kobbefjord from 2008 to 2020, include air temperature (°C), total precipitation (mm), relative humidity (%) and incoming shortwave radiation (W m −2 ).These datasets were used to validate the baseline regional climate model data used as forcing for the C cycle (SPA) model.Similarly, C stocks (from leaf, litter, stems, roots and soil organic matter; g C m −2 ) and foliar N stocks (g N m −2 leaf area) were used to define the initial conditions in the ecosystem model similar to López-Blanco et al. (2018) and López-Blanco et al. (2020).

Downscaled regional climate data
The simulated climate data used here were generated with HIRHAM5 (Christensen et al., 2006), a Regional Climate Model (RCM) based on the combination of the HIRLAM7 numerical short-range weather forecast model (Eerola, 2006) and the physical parameterisation schemes from the ECHAM5 general circulation model (Roeckner et al., 2003).HIRHAM5 has been applied in several studies focusing on Arctic and Greenland permafrost (Rinke et al., 2012), ice sheet surface mass balance (Boberg et al., 2018;Langen et al., 2015;Langen et al., 2017;Mottram et al., 2017), terrestrial CH 4 emissions (Geng et al., 2019) and ecosystem analyses (Schmidt et al., 2019).The model version utilised in this study (Boberg et al., 2018) incorporates a new dynamic snow/ice scheme together with an updated snow/ice albedo scheme (Langen et al., 2015;Langen et al., 2017).To adequately resolve the topography of Greenland, the spatial resolution of the HIRHAM5 RCM grid was set to 0.05°× 0.05°, corresponding to ~5.5 × 5.5 km grid cell sizes (Lucas-Picher et al., 2012).HIRHAM5 is considered to be among the best models for detailed climate modelling in Greenland (Westergaard-Nielsen et al., 2015).The complete model run covers the entire Greenland and Iceland domain (see Fig. 1).We employed HIRHAM5 experiments driven both by the ERA-Interim reanalysis (Dee et al., 2011) for the recent past as in Langen et al. (2017) and by the Global Climate Model (GCM) EC-Earth (Hazeleger et al., 2012) for future projections as in Boberg et al. (2018).Due to the extensive computational resources necessary for the simulations, the datasets available from Boberg et al. (2018) only include five 21-year time slices instead of a full transient coverage of the 21st century.The first year of each time slice was used as a spin-up and therefore not included, resulting in three 20-year time slices.This study employs a baseline time slice (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) and four future slices projecting the middle (2031-2050) and the end (2081-2100) of the century.Each future projection included here follows the IPCC (2013) featuring the Representative Concentration Pathways (RCP) 4.5 and 8.5 scenarios.These represent, respectively, the intermediate and business-as-usual future greenhouse gas emission scenarios without climate policy.We extracted air temperature, total precipitation (rain + snowfall), downward shortwave radiation, wind speed and relative humidity information from HIRHAM5 simulations for the Kobbefjord and Zackenberg fen coordinates using bilinear interpolation.HIRHAM5 air temperature data have been lapse-rate corrected to account for elevation biases similar to Langen et al. (2015).Additionally, we extracted the same variable information from the same locations but from the HIRHAM5 ERA-Interim (hereinafter, ERA-I) reanalysis-driven simulation between 1991 and 2018.The purpose of this was to compare observations with ERA-I-driven HIRHAM5 (which is closely tied to the observed climate evolution) in order to validate the GCM-driven HIRHAM5 (running freely) and thus improve the robustness and certainty of the overall study.

Terrestrial C cycling modelling
The year-round NEE of CO 2 is derived from modelling of gross C fluxes (GPP and R eco ) and their interactions using the Soil-Plant-Atmosphere (SPA) model (Williams et al., 2000;Williams et al., 1996) forced by downscaled climate data from HIRHAM5.SPA is a numerical process model that simulates coupled carbon, water and energy exchanges through ecophysiological principles in vertically resolved multi-level canopy and soil profiles.SPA simulates surface exchanges using I) canopy radiative transfer differentiating between sunlit and shaded foliage following Williams et al. (1998), II) plant photosynthesis simulated using the representation of electron transport and carboxylation reaction as described by Farquhar and von Caemmerer (1982) and sensitive to foliar N, III) surface energy and evaporation fluxes using the Penman-Monteith equation (Jones, 1992) and a stomatal conductance scheme having an explicit coupling between liquid and vapour phase water transport from soil to atmosphere, and IV) vertical exchanges of water and heat across soil layers linked to plant water use, soil evaporation and gravitational drainage.Based on the resulting estimates of photosynthesis, SPA then allocates C into autotrophic respiration and the remaining net primary production into living plant C pools (leaf, labile, stem, fine roots) using a phenology scheme sensitive to weather.Turnover of plant tissues generates inputs into dead organic matter (litter, soil organic matter) pools.Mineralisation of these pools is controlled by soil temperature (López-Blanco et al., 2018).The model has 29 key parameters and seven initial conditions, given the 36 variables (Fig. 6).These parameters can be divided into major groups, those that define: (i) allocation fractions for NPP; (ii) leaf traits, both structural and functional, including foliar N; (iii) turnover rates for live and dead C pools; (iv) hydraulic parameters that determine rates of water supply to leaf tissue; and (v) hydrological parameters that determine water flow through the plant-soil systems.
The SPA model runs at sub-diurnal temporal resolution (30 min), allowing direct comparison against in-situ eddy covariance flux timeseries data.SPA has been extensively calibrated and validated with in-situ C flux and CN stock data from both Kobbefjord (López-Blanco et al., 2018) and Zackenberg (López-Blanco et al., 2020) across the 2008-2018 period (Table S1).Field campaigns dedicated to retrieve in-situ plant trait data such as C and N stocks and plant quality contributed to explain the observed differences in plant C uptake between sites, and generate consistent model outputs with field measurements (López-Blanco et al., 2020).We achieve this by independently calculating maintenance respiration (Rm = Rm leaf + Rm root ) as a function of N based on the parameterisations introduced by Reich et al. (2008) instead of using an arbitrary fixed ratio.Rm in leaves is calculated based on air temperature, leaf N per leaf area, leaf C per area, and leaf area index, while Rm in roots is derived from soil temperatures at 10 cm depth, roots C:N ratio, and roots C stock (López-Blanco et al., 2018).Rm is calculated only when the air temperature > 0 °C.The C cycle runs, forced with HIRHAM5, cover five 20-year periods (one baseline and two future time slides featuring two different RCP scenarios, i.e.RCP4.5 and RCP8.5).The exact atmospheric CO 2 concentration pathway of the RCPs, ended at 538 and 936 ppm, respectively, at 2100.The first five years are regarded as spin-up to allow C stocks and fluxes to equilibrate between the end and start of each time slice.

Benchmarkingmodel evaluations against observations
To create robust ecological projections for these two sites, we first established a solid foundation with a comprehensive climate and ecosystem model benchmarking, including all possible and available in-situ data.Although HIRHAM5 GCM-driven simulations are dynamically self-consistent climate estimations and reconciled with atmospheric properties and physics, their variability is, as is generally the case for a freely-running GCM, out of phase with the actual climate evolution (Boberg et al., 2018) as only the radiative forcing from greenhouse gases and other anthropogenic drivers are specified as boundary conditions.Natural modes of variability are as the models simulate them and hence not to be expected to be in phase with observed modes.A freely-running climate model only takes information about the real-world time evolution of the climate through information on the overall external drivers (e.g., specified concentrations of greenhouse gases following a particular emission scenario) and therefore only represents the statistical properties of the weather (i.e., the climate) at any point in time, not the exact timing of the weather.This approach implies that naturally varying phenomena with an internal time scale from weeks to multiple years will not be in phase with that of the real-world climate system.We therefore also evaluated the performance of HIRHAM5 driven by ERA-I, a global product that combines vast amounts of historical observations using advanced modelling and data assimilation systems (Dee et al., 2011).In this set-up, the HIRHAM5 model is driven by the large-scale structures of atmospheric circulation as well as with the observed sea surface temperature and sea ice cover as they actually occurred.This allows for a direct comparison with the observed weather development and provides a benchmark test of the model's ability to reproduce observed conditions.
First, we compared the day-to-day variability between ERA-I-driven HIRHAM5 output and meteorological observations following the same method as Langen et al. (2015) to ensure quality assurance and consistency.We calculated the goodness-of-fit agreement (Pearson correlation and mean bias) in both Nuuk and Kobbefjord pixels using the exact same years and summer (AMJJAS) and winter (ONDJFM) periods.Second, and after verifying the consistency with Langen et al. (2015), we extended the benchmarking analysis covering the entire baseline HIRHAM5 ERA-Idriven time series used in Langen et al. (2017Langen et al. ( ) (1991Langen et al. ( -2016) ) and HIRHAM5 GCM-driven time series used in Boberg et al. (2018Boberg et al. ( ) (1991Boberg et al. ( -2010)), together with the available in-situ GEM climate data between 1995 and 2016 on Zackenberg and 2007 and 2020 on Kobbefjord.Since only the ERA-Idriven simulation, as argued above, can be meaningfully compared in a day-to-day evaluation against in-situ meteorology observations, the regional downscaled climate products demanded a comparison focusing instead on the probability distributions of the variables.To that end, we compared the two (observation vs model) by plotting their quantiles against each other using QQ-plots (Wilk and Gnanadesikan, 1968).We also generated climatological annual cycle plots averaged over multiple years to characterise the seasonal climate and flux variability of air temperature, incoming shortwave radiation and vapour pressure deficit as well as the different C fluxes (NEE, GPP and R eco ), respectively.

Uncertainty quantification and sensitivity analysis
To assess how critical each of the manually calibrated ecosystem parameters is to the estimated C fluxes in the future, we performed a sensitivity analysis over the 36 parameters tuned in SPA using the method described by López-Blanco et al. (2018).Such analysis helps us not only to identify model limitations in the C cycle (i.e.C allocation, C turnover and plant phenology) but also to evaluate the robustness of the model outputs in the face of uncertainty.Here we summed and propagated an uncertainty ensemble over the forecasted period of the combined variability introduced by 1) the two different emission scenarios (RCP4.5 and 8.5; U RCP Scenario ) and 2) the inherent model limitation of using stationary parameters (U Parameters ) that can (and likely will) change over time.The total uncertainty of the annual sums of NEE, GPP and R eco (U Flux ) was calculated using the error accumulation principle (Lund et al., 2012) until the end of the 21st century: To calculate the error emerging from the model parameters, we first used the 36 nominal parameters previously determined for Kobbefjord (López-Blanco et al., 2018) and Zackenberg (López-Blanco et al., 2020).We then modified each individual parameter ±10 % in sequence, ran the 1991-2010, 2031-2050 and 2081-2100 periods for both RCP scenarios and ranked from high to low sensitivity the percentage change in total annual C sink strength (NEE), our selected response variable.The ratio of the % change in NEE to the % change in each parameter is the so-called sensitivity index (SI).|SI| > 1 (|SI| = magnitude of SI) means that the parameter is sensitive to the response variable.We performed the sensitivity analysis across the entire time series, but we only compared the 1995-2010 and 2085-2100 periods including both RCP4.5 and 8.5 emission scenarios to account for the total difference between the baseline and the end of the century runs.
Although climate variability contributes to the C sink strength at both sites, a previous study has argued that local plant nutrient status may be a more important factor than climate for determining carbon flux magnitudes between 2008and 2018(López-Blanco et al., 2020).To further test the importance of climate together with local plant N status over the net C uptake, we arranged a simple experiment interchanging the richer value of foliar N measured in Zackenberg (2.26 g N m −2 leaf area) with that at Kobbefjord (1.61 g N m −2 leaf area).We ran the experiments across the entire dataset comparing the 1995-2010 and 2085-2100 periods using both future scenarios and compared the probability distribution function of annual mean NEE with nominal runs to assess the impact of foliar N changes and climate over the C sink strength.

Projected climate variability and uncertainties
HIRHAM5 projects a clear overall increase in air temperature, total precipitation and vapour pressure deficit (VPD) and a decrease in the incoming short-wave radiation by the end of the 21st century (Fig. 2).The degree of change differs significantly, however, depending on climate variable, site and RCP emission scenario.The average annual air temperature in Kobbefjord is projected to increase from −5.8 °C in the baseline period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) to −0.8 °C between 2081 and 2100, a warming of 5 °C, while Zackenberg will also experience an increase from −12.0 °C to −4.3 °C (a 7.7 °C temperature increase).On average, the RCP8.5 scenario estimates a 2.7 °C warmer temperature than the RCP4.5.The precipitation tendency per site is similar -Zackenberg will face a 110 % increase in total precipitation compared to a much lower (but still significant) increase in Kobbefjord of 19 %.On average, the RCP8.5 scenario estimates 10 % wetter conditions than the RCP4.5 scenario in both Kobbefjord and Zackenberg.Despite this, the magnitude of precipitation will continue to be ~6 times larger in Kobbefjord than in Zackenberg, which is similar to current patterns.Interestingly, the HIRHAM5 projections indicate a significant lengthening of the snow-free periods, on average occurring 9 days earlier in spring and 10 days later in autumn in Kobbefjord (totalling 19 days more without snow), while in Zackenberg the projections show snow-free conditions three days earlier in spring and onset of snow cover only two days later in autumn.According to HIRHAM5 runs, the incoming shortwave radiation used by plants will likely get reduced by 9 % and 6 % in Kobbefjord and Zackenberg, respectively.This behaviour is likely associated with the increase in precipitation and therefore increase in cloud cover and darker conditions.Finally, both sites will experience a similar increase in VPD (32 % in Kobbefjord and 46 % in Zackenberg), implying that both sites will likely hold more moisture in the air when it is saturated.Overall, the RCP8.5 scenario estimates, on average, 3 % darker conditions and 10 % higher saturation vapour pressure compared to the RCP4.5 projection, respectively.Fig. 2. Recent past and future 5-year mean air temperature (°C), precipitation (mm/ year), shortwave radiation (W m −2 ) and VPD (hPa) in Kobbefjord and Zackenberg estimated by HIRHAM5 Regional Climate model.The colored envelopes denote the modelled 5-year standard deviation.The RCP 4.5 scenario future projections are displayed as dashed lines, while the RCP 8.5 scenario is presented with solid lines.Note that these are results from the freely-running climate model, and the recent past (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) does therefore not correspond to the real-world time evolution.
The agreement between HIRHAM5 ERA-I-driven (assimilated from observations) and GCM-driven estimates (used for projection purposes) exhibits a reasonable predictive performance compared to 20 and 12 years in-situ observations from Zackenberg and Kobbefjord, respectively (Figs. 3 and S2, Table S2).HIRHAM5 ERA-I-driven simulations perform relatively well estimating contemporary air temperatures (Pearson correlation (corr) = 0.93 and mean bias (MB) = −0.94°C; Table S2) and shortwave radiation (corr =0.95 and MB = 10.9W m −2 ) but less so with precipitation (corr = 0.66 and MB = 0.53) and vapour pressure deficit (corr = 0.60 and MB = 0.52 hPa) (Table S2 and Fig. S2).As expected, the HIRHAM5 ERA-I-driven runs are in close agreement with the in-situ GEM field observations (Table S2).The ERA-I driven experiment follows the actual evolution of the weather (Dee et al., 2011), while the GCM driven experiments can be seen to represent the time evolution only in a statistical sense (see also Section 2.3.2).Yet, HIRHAM5 GCM is somewhat consistent with HIRHAM5 ERA-I simulations and in-situ observations (Table S2).HIRHAM5 GCM simulates significantly better temperature and radiation (often easier to predict and tied to largerscale processes and air masses) than precipitation and humidity (dependent on more localised processes that are inherently less predictable) (Fig. S2 and Table S2).The daily-aggregated seasonal variability of temperature, incoming radiation and humidity from both sites is also reasonably captured by both simulations (Fig. 3).There are, however, clear shifts in magnitude and timing (Table S2).For example, temperature shows a cold bias (−0.9 to −3.2 °C average range), especially in winter when simulated temperatures are colder than in-situ observations.Downward shortwave radiation is systematically higher compared to GEM observations (4-11 W m −2 ), suggesting an underestimation of cloud cover at both sites, especially in Kobbefjord.Humidity presents a clear lower degree of agreement, underestimating VPD in Kobbefjord by 0.9-1.0hPa and overestimating it in Zackenberg by 0.7-0.9hPa, which points to a divergent bias.

Projected C flux variability and uncertainties
The SPA model forced with HIRHAM5 climate data estimates a baseline net sink of C of −16 ± 14 g C m −2 year −1 in Kobbefjord and −52 ± 41 g C m −2 year −1 in Zackenberg (Fig. 4), values that are consistent with NEE fluxes measured in the field between 2008 and 2018 (~−18 and ~−50 g C m −2 year −1 , respectively; López-Blanco et al., 2020).The combined climate and ecosystem models suggest that the projected future increases in temperature and precipitation will intensify the net C sink strength of Kobbefjord by 10 ± 4 g C m −2 year −1 and Zackenberg by 9 ± 1 g C m −2 year −1 , on average, projecting the annual mean NEE to −26 ± 18 and −61 ± 41 g C m −2 year −1 by the end of the 21st century, respectively.In other words, Zackenberg is expected to have more than twice a C sink strength capacity compared to Kobbefjord, a difference in magnitude that is again consistent with recent estimates based on 20+ years of in-situ CO 2 measurements from both fens (López-Blanco et al., 2020).Photosynthesis will grow from the baseline period around 47 % in Kobbefjord and 62 % in Zackenberg, while respiration will also increase by 44 % and 77 % in the same two sites (Fig. 4).The estimated annual mean GPP and R eco for Kobbefjord at the end of 2100 are −160 ± 27 and 134 ± 21 g C m −2 year −1 , while in Zackenberg the same gross fluxes are expected to be slightly higher than twice this at −345 ± 52 and 284 ± 39 g C m −2 year −1 , respectively.Overall, by 2100, the plant uptake in Kobbefjord will increase more strongly than the ecosystem respiration, while Zackenberg will experience the opposite, i.e. stronger respiration than photosynthesis increases.Zackenberg has a more than two times stronger C sink strength capacity than Kobbefjord, thus respiration in Zackenberg can expect a higher increase than photosynthesis and not only still be a net C sink, but also a stronger C sink strength than Kobbefjord.However, the photosynthetic inputs in both sites will continue to dominate the respiratory outputs, and therefore the overall C sink strength is expected to become stronger.The fact that GPP will have a stronger influence than R eco in both sites is again consistent with previous findings based on field observations (López-Blanco et al., 2020).It is worth mentioning that the combination of the SPA model forced with HIRHAM5 climate also exposes a rather large inter-annual variability.For example, for a period of a few years, the ecosystems in both Kobbefjord and Zackenberg shifted from a consistent sink of C to an episodic C source (Fig. 4).This is not unexpected such sudden shifts in the C sink/source functioning due to climatic and biological episodic events have been previously reported for both sites (e.g.Christensen et al., 2020;López-Blanco et al., 2017;Lund et al., 2017).Additionally, the estimated seasonal C fluxes across the baseline period for both Kobbefjord and Zackenberg are, respectively, 3.4 ± 0.4 and 3.0 ± 0.6 g C m −2 in winter, 10.2 ± 1.4 and 16.7 ± 3.2 g C m −2 in spring, −38.1 ± 10.4 and −80.4 ± 35.1 g C m −2 in summer and 8.7 ± 1.5 and 8.4 ± 1.8 g C m −2 in winter (Fig. 5).The SPA model additionally quantifies an increase in C sink strength of −22.9 g C m −2 in Kobbefjord and −48.4 g C m −2 in Zackenberg only during summer (JJA) and an increase in the respiratory losses jointly in autumn (SON), winter (DJF), and spring (MAM) of about 54 % in Kobbefjord and 144 % in Zackenberg by the end of the century.
In general, the comparison between the SPA ecosystem model forced with HIRHAM5 (GCM-driven) climate data and the in-situ GEM observations reported in Kobbefjord by López-Blanco et al. (2017) and in Zackenberg by López-Blanco et al. (2020) shows that the seasonal C flux variability is captured by the simulations.The degree of correlation between measured and simulated climatological annual cycles of NEE (0.89-0.91),GPP (0.94-95) and R eco (0.86-0.94) is acceptable (Fig. S3).There are significant mismatches in magnitude (e.g.Kobbefjord GPP and R eco ; MB = 0.71, and −0.73 g C m −2 growing season −1 , respectively) and growing season start/end (e.g.Zackenberg GPP; MB = 0.38 g C m −2 growing season −1 ).In this implementation, snow cover information from the in-situ cameras is not used to constrain the SPA model simulation of vegetative phenology and the length of the growing season, decreasing the overall performance compared to López-Blanco et al. (2018) and López-Blanco et al. (2020).Furthermore, shortcomings related to the HIRHAM5 GCM forcing data further intensify the uncertainties from the SPA simulations on top of the intrinsic model limitations.Even with such obvious biases, SPA can estimate similar modelled C sink strength (−29.3 g C m −2 in Kobbefjord and −56.9 g C m −2 in Zackenberg) compared to measured fluxes in the field (−29.7 g C m −2 in Kobbefjord and −65.6 g C m −2 in Zackenberg) during the June-September period, ensuring that the SPA model forced with HIRHAM5 data is consistent with local measurements.

Climate and plant traits sensitivities driving the C sink strength
To better understand whether local C and N conditions or long-term climate variability have a bigger impact on the C sink strength, we calculated the relative sensitivity of all vegetation-related parameters used in the SPA model (Fig. 6).The sensitivity indices (SI) exposed significant differences between 1) model parameters, 2) the baseline (1995-2010) and future (2085-2100) periods and 3) sites.As an example, the leaf mass per area SI in the Zackenberg baseline run (SI-NEE = 2.25; Fig. 6) implies that if this specific parameter increases only by 1 %, the forecasted NEE flux will experience a 2.25 % change.Overall, parameters related to plant traits such as net primary production (NPP) allocated to roots, maximum foliar carbon stock, leaf mass per area, rate coefficient for J max (which is involved in the maximum rate of electron transport during photosynthesis) and average foliar N are the most sensitive parameters controlling the net C exchange in both sites.Our sensitivity analysis further revealed different sensitivity patterns between the beginning and the end of the centuryfor instance the 1995-2010 period showed sensitivities with larger variabilities (coefficient of variation (CV) = 56-76 % range) and responsive ranges between sites than the 2085-2100 period (CV = 47-52 %).Additionally, the future period was significantly less sensitive to the baseline period (51-54 % decrease, on average) and had similar ranked sensitivities between sites and emission scenarios.
By switching the richer input foliar N from Zackenberg in Kobbefjord, and vice versa, we further quantified the net effect of climate change vs local N conditions over the C sink strength of these contrasting sites.On the one hand, the SPA model suggests that using the N-richer conditions observed in Zackenberg (2.26 g N m −2 leaf area instead of 1.61 g N m −2 leaf area) could strengthen the net C sink in Kobbefjord by ~8 g C m −2 year −1 during the baseline run and ~10 g C m −2 year −1 at the end of the century (Fig. 7).Interestingly, a nutrient-richer Kobbefjord in the present day (annual mean −25.7 g C m −2 year −1 ) would not be as productive as the expected C sink strength at the end of the century influenced only by climate (−32.7 g C m −2 year −1 ).On the other hand, the experiment in Zackenberg shows a different hypothetical scenario.A weakened N input in Zackenberg will drastically reduce its C uptake strength by ~21 g C m −2 year −1 compared to the baseline N status and ~34 g C m −2 year −1 from the forecasted net C uptake driven only by climate at the end of the century (Fig. 7).It is worth mentioning that a theoretical decrease of C sink strength due to the poorer N conditions today could be counterbalanced and even slightly increased by the effect of climate change towards the end of 2100 (this is −52.1 g C m −2 year −1 vs the baseline run average of −48.8 g C m −2 year −1 ).Note that these values are consistent with the annual mean of −50 g C m −2 year −1 measured in-situ across 10 + years of observations (López-Blanco et al., 2020).

How sensitive is the C balance expected to be under warmer and wetter conditions forecasted in the 21st century?
The climate projected at Kobbefjord and Zackenberg by HIRHAM5 based on the RCP8.5 scenario suggests an overall increase in air temperature (5-7.7 °C), total precipitation (19-110 %) and vapour pressure deficit (32-46 %) by the end of 2100, while incoming shortwave radiation is expected to decrease (6-9 %) (Fig. 2).There is good evidence of the recent past and present warming and wetting of the Arctic (AMAP, 2022).First, air temperatures have increased three times more in the Arctic compared to the globe; this is a 2.8 °C increase between 1979 and 2019 based on the monthly EU Copernicus ERA5 reanalysis dataset (Hersbach et al., 2020).Second, the total precipitation likely increased by 9 % (mostly dominated by a 24 % rainfall increase) during the same period according to the same dataset.Third, a multi-snow-product ensemble created by Mortimer et al. (2020)  More importantly, the amplified Arctic warming will likely continue.The latest future projections using climate models from the Coupled Model Intercomparison Project Phase 6 (CMIP6) (O'Neill et al., 2016) point to intensified warmer and wetter conditions.On the one hand, 20+ CMIP6 models expect an average increase of 5.8 °C and 10.4 °C for the Arctic annual mean under the new Shared Socio-economic Pathways SSP2-45 and 5-85 scenarios (roughly equivalent to RCP4.5 and 8.5 used in this study), respectively, by the end of the century (Cai et al., 2021).Our estimates point to an annual increase of 2.7-5 °C in Kobbefjord and 4.5-7.7 °C in Zackenberg for RCP4.5 and 8.5, respectively (Fig. 2).This warming will increase the snow-free periods by 5-19 days in Zackenberg and Kobbefjord.On the other hand, larger (~40 %) and faster changes in the hydrological cycle are expected due to greater poleward moisture transport and increased sensitivity of precipitation to Arctic warming than previously anticipated (Bintanja et al., 2020;McCrystall et al., 2021).Our assessed sites showed a tendency towards wetter conditions, although with a much greater increase in Zackenberg (110 %) than in Kobbefjord (20 %) (Fig. 2).Directly linked with this process was the projected 6-9 % decrease of shortwave radiation in Zackenberg and Kobbefjord, which may be associated with a potential increase in cloud cover during the 21st century (Fig. 2).Clouds are an important control of the Arctic climate (Huang et al., 2021), and they are key regulators of energy balance processes (i.e.trapping warm temperatures) and other ecosystem processes influencing the amount of solar radiation received by the surface (Dong et al., 2010).The fact that air temperature will continue to increase globally will trigger an intensification of water evaporation and lead to a higher moisture content in the atmosphere, leading to an overall increase in precipitation.However, changes in local cloud cover should be seen as a combined effect of local forcing and long distance advection of water vapour as well as of clouds.Additionally, VPD is also expected to increase towards 2100 (32-46 %; Fig. 2), which agrees with empirical evidence that VPD is increasing globally (Yuan et al., 2019).VPD and incoming shortwave radiation are essential variables determining plant photosynthesis, and potential increases (Yuan et al., 2019) or decreases (Christensen et al., 2020) may weaken plant photosynthesis and the subsequent C sink-source functioning.According to our data, the fact that plants will have less light to photosynthesise combined with the presence of drier air above the surface seem to be compensated by the increase in temperatures and total precipitation (Figs. 2 and 4).
The likely warming and wetter future scenario (Fig. 2) will trigger substantive effects on the C cycle of the two climatically different tundra sites of Kobbefjord and Zackenberg according to the SPA ecosystem model (Fig. 4).Thus, the combined change in climate conditions will intensify the net C sink strength of Kobbefjord to −26 ± 18 g C m −2 year −1 and to −61 ± 41 g C m −2 year −1 in Zackenberg towards the end of the 21st century.This is an increase in C storage of 10 ± 4 and 9 ± 1 g C m −2 , respectively, compared to what is measured in the field today (López-Blanco et al., 2020).Our projections suggest that, even though Kobbefjord will increase its NEE more than Zackenberg in absolute numbers, the gap in C sink strength magnitude between Kobbefjord and Zackenberg will not be bridged as its difference (Δ) will remain within the 35-36 g C m −2 year −1 range.This coupling is associated with the compensatory effect between photosynthesis and ecosystem respiration (López-Blanco et al., 2017;Richardson et al., 2007;Wohlfahrt et al., 2008) and is consistent during the entire 21st century (Fig. S5).Although coupled, GPP and R eco fluxes are expected to be larger in Zackenberg (ΔGPP = 132 ± 55 and ΔR eco = 123 ± 39 g C m −2 year −1 ) than in Kobbefjord (50 ± 24 and 41 ± 18 g C m −2 year −1 ) (Fig. 4).It is therefore important to highlight that the future GPP will dominate over R eco in both sites, and controls on photosynthesis will thus play a greater role in driving the C storage in the coming years.This finding is in line with the widespread greening over the past few decades in response to warmer and longer summers (Myers-Smith et al., 2020;Myneni et al., 1997;Berner et al., 2020) but also with future productivity projections that account for change in plant trait responses to climate (Madani et al., 2018).Nevertheless, this general unbalance favouring the net C uptake will likely experience sporadic sudden shifts from net sink to source (see Fig. 4).Arctic warming will not necessarily lead to systematic increases in the C storage as there are several observed examples exemplifying different consequences of episodic climate-related events.Examples of such instabilities are those triggered by delayed onset of snowmelt and intense rain events combined with thermokarst erosion (Christensen et al., 2020), vegetation damage by extreme winter events (Parmentier et al., 2018) or biological disturbances facilitated by climatic conditions (Lund et al., 2017), among others.
Overall, the baseline Kobbefjord and Zackenberg sink power from HIRHAM5-SPA (−38.1 and −80.4 g C m −2 , respectively) is well aligned with the June to August range (−119.0 to −5 g C m −2 ) from 10 other Arctic wetlands at similar latitudes reported by Coffer and Hestir (2019) (Table 2).Additionally, our estimates from this modelling exercise are closely related to the values measured at both sites within the same period presented by López-Blanco et al. (2020).In a broader tundra context, our baseline annual NEE averages are stronger sinks than the −3.3 ± 44.2 g C m −2 year −1 annual average but within the monthly NEE range of −125-+75 g C m −2 estimated from an extensive 136-sites and 2217monthly dataset, only featuring in-situ observations from 1995 to 2020 recently reported by Virkkala et al. (2022).Virkkala and colleagues quantified the winter (DJF), spring (MAM), summer (JJA) and autumn (SON) periods in 9 ± 10, 6 ± 9, −26 ± 38 and 10 ± 21 g C m −2 , respectively.Our modelled baseline estimates in Kobbefjord (3.4 ± 0.4, 10.0 ± 1.3, −38.1 ± 10.5, and 8.7 ± 1.5 g C m −2 ) and Zackenberg (3.0 ± 0.6, 16.7 ± 3.2, −80.4 ± 35.1, and 8.4 ± 1.8 g C m −2 ) are consistent with the same climatological periods.Interestingly, our model additionally projects an intensification only during summer of the C sink strength up to −61.0 ± 12.4 g C m −2 and −128.9 ± 32.7 g C m −2 , respectively, but also an increase in the respiratory losses during the shoulder seasons (winter, spring, and autumn) of ~54 % in Kobbefjord and 144 % in Zackenberg by the end of the century.These numbers suggest that, even though an increase in respiration in the cold seasons is expected, summer processes will keep controlling the increases in total net C storage.

Is the C sink strength more sensitive to local carbon and nitrogen conditions than to climate variability?
Our modelling study shows a marked effect of the future climate change leading to increased terrestrial C storage.But the ecological question of whether future climate variability or local differences in N conditions will have a higher or lower contribution to the overall terrestrial C exchange is left open.
The Arctic is known for its well-documented N-limited conditions and its constraining effects on tundra plant productivity (Chapin III and Shaver, 1985;Shaver et al., 1992).Although there is a clear poleward decrease of biological N fixation and N mineralisation (Du et al., 2020), there is strong Present and future (based on RCP8.5 scenario) annual mean C sink strength from Kobbefjord (brown) and Zackenberg (green) comparing the control setups using the nominal parameter set and two experimental setups including enhanced and weakened foliar N inputs from the opposite site.Similar results using the RCP4.5 scenario are displayed in Fig. S4.
spatial N variability due to specific local differences such as climate boundary conditions, topography, geological substrate, flora and fauna.In Greenland, the potential limitations of higher latitudes featuring colder temperatures and shorter growing seasons controlling the net C uptake have been found to be counterbalanced by richer plant tissue N and soil nutrients (López-Blanco et al., 2020).This study based on 20+ years of in-situ CO 2 measurements assessed how the difference in nutrient availability, and not climate, can substantially enhance the C sink strength in Zackenberg fen compared to Kobbefjord.The larger foliar N stocks in Zackenberg enhance the overall plant productivity (Arndal et al., 2009).Plant N stocks are likely higher due to the presence of large grazers that further increase N concentrations and plant quality (Mosbacher et al., 2019) and the large abundance of mosses with high N content (Street et al., 2012).The higher availability of dissolved organic carbon and nutrients in the soil water is also associated with potential lateral flow from the adjacent slopes dominated by cation and nutrient-rich fast-weathering basalts and sediments (Cable et al., 2018).There is recent evidence that geomorphology may control nutrient and organic matter exports via flow paths, soil organic stocks and nutrient retention by plants (Pastor et al., 2021).Additionally, permafrost soils such as in Zackenberg might help to release more plant-available N but also retain N availability better than permafrost-free soils like Kobbefjord (Olefeldt et al., 2014;Reyes and Lougheed, 2015).Atmospheric N deposition in high Arctic systems is below 0.2 g N m −2 year −1 (van Cleve and Alexander, 1981), and ~0.03 g N m −2 year −1 specifically from Greenland (Lamarque et al., 2011), while biological N-fixation has been reported at ranges between 0.08 and 0.13 g N m −2 at another tundra site (Hobara et al., 2006).Since soil N mineralisation rates are ~0.35g N m −2 per growing season in a wet sedge tundra (Schmidt et al., 2002), we may expect that even having low deposition rates these tundra systems could undergo net gain of N. It is, however, unclear whether the Arctic is (and will become) (Chapin III et al., 1995;Martin et al., 2022) limited by N in response to elevated temperatures.It is also unclear what impact on future productivity and carbon storage this might have.The full implications and complete N budget partitioning (i.e.biological N-fixation, plant N uptake, mineralisation rates, microbial immobilization, volatilization, etc.) remain uncertain and unquantified as complex data are required to clarify the processes involved.This is not a trivial matter as there is recent evidence that models can underestimate the land C uptake by up to 20-21 % when nutrient limitation is accounted for in projected end-of-century estimates (Wieder et al., 2015;Meyerholt et al., 2020).Fully coupled C-N models can account for more complex ecosystem feedback loops such as the fate of N availability as response to warming.However, recent studies by Arora et al. (2020) and Spafford and MacDougall (2021) reviewed the validation of terrestrial biogeochemistry in the most modern CMIP6 Earth System Models (ESM) and noted that 5 of the 11 assessed models do not account for explicit N cycling.Despite the improvement compared to the previous CMIP5 (only two out of eight ESMs had terrestrial N coupled with their C cycle) there is large and unconstrained uncertainty in the magnitude of the regional N fluxes in the (Davies-Barnard et al., 2020).Our analysis cannot answer what the fate of N availability in the Arctic will be, but given that our modelling framework 1) takes into account future potential changes in plant N trait parameters by adding an uncertainty ensemble (Fig. 4) and 2) is consistent with the observed C cycle status in the field (and constrained by local foliar N; Fig. S3), we quantify the relative importance of climate change compared to the potential shifts in local N conditions in shaping the C uptake.Firstly, our sensitivity analysis revealed not only that plant traits are the most sensitive parameters controlling the net C exchange in both sites but also that different sensitivities exist between the beginning and the end of the century.For instance, the baseline 1995-2010 period displayed parameter sensitivities associated with larger variabilities and responsive ranges between sites compared to the 2085-2100 period (Fig. 6).The influence of climate change by the end of the century seems to have a similar impact on C sink strength to the plant trait sensitivity between sites.The fact that essential parameters in charge of photosynthetic efficiency, vegetative phenology and initial CN status are more sensitive in both baseline and future periods suggests that the forecasted increase in C sink strength will likely be as influenced by future changes in climate as the existing local N conditions exposed to a potential (although arbitrary) ±10 % change.We believe that trait sensitivity is an emergent ecosystem property, which highlights the potentially critical effect of changes in environmental conditions on trait shifts in tundra ecosystems (Bjorkman et al., 2018) and what impact changes in plant traits will have on terrestrial C uptake (Madani et al., 2018).
Secondly, we set up two experiments forcing realistic changes in plant nitrogen status based on in-situ information to study the implications of potential increase/decrease of plant N traits on NEE.The observed stronger carbon uptake and exchanges in Zackenberg were associated with additional nutrient-rich plant tissues (Table 1) and will likely prevail towards the end of the century even when we accounted for a 10 % change in average foliar N in SPA (Figs. 4 and 7).If the foliar N input decreases (29 %) from the observed 2.26 g N m −2 in Zackenberg (Arndal et al., 2009;Mosbacher et al., 2019;Street et al., 2012) to the 1.61 g N m −2 in Kobbefjord (López-Blanco et al., 2018), the net ecosystem exchange nearly halved from −48.8 to −28.0 g C m −2 .Interestingly, the potential booster effect of climate change under N-poorer conditions may overrule the net C uptake by surpassing the baseline average (Fig. 7).This finding suggests that realistic foliar N changes have impacts on C cycling of similar magnitude to climate change.Furthermore, the hypothetical addition of N-rich foliage content from Zackenberg increased the C sink strength of Kobbefjord from −18.1 g C m −2 to −25.7 g C m −2 compared to the forecasted −32.7 g C m −2 expected by the sole effect of climate change (Fig. 7).This finding suggests similar impacts from N changes versus climate change.Interestingly, Rasmussen et al. (2022) recently found that the overall net C

Table 2
Mean and standard deviation of summer (JJA) observations of net ecosystem exchange (NEE) in g C m −2 per month reported by Coffer and Hestir (2019), López-Blanco et al. (2020), and Virkkala et al. (2021) 4 and 7).Our modelling exercise demonstrates that the potential contribution from viable plant N trait variations on the future C sink strength is clearly important and could be as much as climate.

Present model uncertainties
There is a pressing need for efficient coordination and integration between communities in charge of long-term monitoring and ecosystem modelling to reduce uncertainties in current diagnoses, and future projections and associated feedbacks (Friedlingstein et al., 2014;Ahlström et al., 2012;Luo et al., 2015;IPCC, 2019).By coordinating, more informed and constrained projections will anticipate and guide adaptation to future conditions.Our study strongly emphasises the importance of integrating ecosystem models with long-term monitoring data in a typically understudied region of the Arctic.We use extensive and quality-checked in-situ datasets closely related to modern climate and ecosystem models to not only establish a solid foundation for model calibration and validation (Figs. 3, 5, and S2, Tables S1 and S2) but also to allow more credible ecological projections and sensitivity analyses (Figs. 2,4,6,7,and S2).We explored relationships between plant traits and the environment (Bjorkman et al., 2018) in order to improve estimates of ecosystem change (Wullschleger et al., 2014).It is, however, necessary to acknowledge inherent limitations on climate and C cycle projections generated by HIRHAM5 and SPA as they are associated with biases.
First, in order to evaluate future changes of C sink-source functioning, we can only rely on GCM-driven simulations.GCMs are dynamically selfconsistent climate models that mimic boundary conditions from actual atmospheric properties and physics although their variability, by design, is not in phase with reality (Boberg et al., 2018).In other words, GCMs cannot (and will not) predict the exact weather/climate that will occur tomorrow or in the next century.This is why it was essential to evaluate contemporary meteorological values using in-situ datasets and ERA-I-driven simulations, which are in phase with in-situ observations (Dee et al., 2011).Although the ERA-I HIRHAM5 runs adequately capture the climate variability observed in the field (Figs. 3 and S2, Table S2), ERA-I-driven estimates are not exempted from biases and uncertainties.For example, the HIRHAM5 ERA-I-driven simulation is known to have a cold bias in West Greenland, particularly in winter, but also a positive bias in shortwave radiation (Langen et al., 2015), which is consistent with the underestimation of cloud cover.Our bias numbers for both sites are consistent with these reported tendencies, showing a negative 0.94 °C temperature bias and a positive 10.9 W m −2 bias in downward shortwave radiation (Fig. 3 and Table S2).Even with such apparent magnitude mismatches, the correlations between observations and model (0.93 and 0.95, respectively) are reasonably captured, suggesting a fair degree of seasonal agreement.Langen et al. (2015) additionally reported large positive precipitation biases (~+200 %) near the ice sheet in Kapisillit (60 km from Kobbefjord) and significantly smaller biases closer to the sea (−20 % to +70 %) in the surroundings of Kobbefjord.This is again consistent with our findingsprecipitation modelled by HIRHAM5 in Kobbefjord and Zackenberg were, respectively, ~60 % and ~65 % higher than the annual observations (Table S2).This offset could partially be explained by the uncertainties of the observed precipitation.The issues of wind-induced undercatch of solid (but also liquid) precipitation in gauges (Yang et al., 1999) are well documented, and precipitation in West Greenland is typically adjusted with 1.4-1.6 (40-60 %) correction factors (Mernild et al., 2015).By applying such correction, our results suggest that precipitation might be reasonably well captured by HIRHAM5 in West Greenland too.The mismatch detected for VPD seems more complex as it overestimates in-situ observations in Zackenberg and underestimates those in Kobbefjord.The fact that precipitation is substantially larger in Kobbefjord (Fig. 2) might lead to wetter soil conditions in the model associated with higher local evaporation and lower near-surface VPD, while the opposite effect may occur in Zackenberg due to the much drier conditions.More dedicated analyses focusing on these non-trivial biases in humidity are encouraged to better understand spatiotemporal mismatches between data and models.
Second, the ecosystem model set-up used in this study gained knowledge from previous supporting studies.For instance, we implemented parameterisation and calculated maintenance respiration losses considering local nitrogen conditions based on a strong respiration-nitrogen relationship in both leaves and roots using a modified version of the Reich et al. (2008) calculation, which enabled us to avoid fixed ratios and therefore understand the implications of N in relation to C fluxes.We also relied on multiple time series of climate, C flux and CN stock data from the longterm monitoring GEM programme (https://data.g-e-m.dk/) to assist an exhaustive calibration and validation of the HIRHAM5-SPA modelling framework.Previously, as described in López-Blanco et al. (2018), we found that the model's most sensitive parameters to C fluxes were those related to leaf N traits and initial C stocks, most of which were derived from in-situ information that subsequently improved the certainty of the modelled runs.This study is now built upon the same uncertainty evaluation but also accounts for parameter uncertainties from both sites including two different future emission RCP scenarios.However, there are intrinsic limitations of this prognostic modelling exercise using HIRHAM5 GCM forcing compared to previous modelling exercises using in-situ climate data.Biases from the HIRHAM5 climate likely contribute to underestimating both GPP and R eco , particularly in spring and autumn.This limitation may intensify the already existing biases attributed to SPA alone forced with in-situ GEM climate where the C sink strength is underestimated by 1-2 g C m −2 month −1 in May and 3.5-5 g C m −2 month −1 in October.Additionally, the snowmelt time is of critical importance in the Arctic for vegetative and reproductive phenology (Assmann et al., 2019) and the net C balance (Lund et al., 2012).In the past, SPA has been improved to estimate the growing season start by resetting the growing degree day (GDD) based on snowmelt time data inferred from automatic cameras monitoring the snow cover percentage (Westergaard-Nielsen et al., 2017;López-Blanco et al., 2018).The SPA model version presented here does not reset GDD based on the snowmelt time but on the surface temperature, and therefore our estimation of the timing of the start of the growing season (and subsequent net C uptake timing) can be improved (e.g.Fig. S3).Finally, a thick snowpack may promote a warming effect on soil temperatures (Lund et al., 2012), enhancing heterotrophic respiration and decreasing the C sink strength (López-Blanco et al., 2018).The SPA snow cover subroutine to update the soil surface energy balance based on Essery (2015) was not used as data from the cameras were not available, and a potential transformation from HIRHAM5 snowfall into snow cover could create additional biases.

Conclusions
In this paper, we projected the future terrestrial C sink strength of two ecosystem stations in Greenland in connection with climate and nutrient interactions.We used multiple in-situ data streams measured by the Greenland Ecosystem Monitoring (GEM) programme integrated with the Soil-Plant-Atmosphere (SPA) ecosystem model and climate projections from the highresolution HIRHAM5 model.Based on our findings, we conclude that: 1.There is evidence suggesting that temperatures, total precipitation and vapour pressure deficit will increase, while incoming shortwave radiation will decline at both sites by the end of the 21st century, with a larger amplitude the higher the level of global warming.2. Such a combined effect will, on average, intensify the net C sink strength of Kobbefjord by 10 ± 4 g C m −2 year −1 and 9 ± 1 g C m −2 year −1 in Zackenberg, consequently reaching −26 ± 18 g C m −2 year −1 and −61 ± 41 g C m −2 year −1 towards 2100, respectively.
3. In relative terms, Kobbefjord will increase plant C uptake more than respiration losses, while Zackenberg will experience the opposite.However, the overall photosynthetic inputs in both sites will continue to dominate the respiratory outputs, and the overall C sink strength is therefore expected to increase.4. Our sensitivity analysis reveals that plant traits are the most sensitive parameters controlling the net C uptake in both sites.Furthermore, it shows that the forecasted increase in C sink strength will likely be as influenced by future changes in climate as the potential adjustments in plant traits.This emergent property indicates that the forecasted strengthening of the net C sink will be closely controlled by the future climate change and the existing local conditions.5.A series of experiments forcing larger but realistic changes in plant nitrogen traits at both sites corroborated this hypothesis.6.Multiple challenges and uncertainties remain.This work establishes robust baselines for 1) model calibration and validation, 2) pinning down uncertainties and 3) future upscaling exercises at different spatial scales focusing on Greenland.
Overall, we believe that it is vital to find and explore the common ground between the local-scale field observations and the coarse resolution climate projections focusing on Greenland, which are typically omitted in global modelling analyses due to complexity and lack of data.Likewise, this work may benefit the next generations of the Greenlandic population as more constrained and informed projections will guide adaptation to the future implications of feedback effects from climate change.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

•
Climate change will have unclear consequences on Arctic C sink-source functioning.• Climate and ecosystem models are the only tools able to predict climate change.• HIRHAM5 and SPA models anticipate stronger C storage in Greenland by 2100.• Foliar N changes will have impacts on net C uptake of similar magnitude to climate change.

Fig. 1 .
Fig. 1.Example of downscaled air temperature from the HIRHAM5 GCM-driven Regional Climate model (DMI) featuring Greenland on 1 August2000 at a 5 × 5 km spatial resolution.[Z] and [K] maps zoom in on the Zackenberg and Nuuk-Kobbefjord areas surrounding each research station (red dots) from where climate data were extracted to force the SPA ecosystem model between 1991-2010, 2031-2050 and 2081-2100.

Fig. 3 .
Fig.3.Daily-aggregated seasonal variability of temperature, radiation and vapour pressure deficit from in-situ GEM observations and HIRHAM5 ERA-I-and GCM-driven outputs in the Kobbefjord and Zackenberg sites.

Fig. 4 .
Fig.4.Recent past and future 5-year annual mean net ecosystem exchange (NEE), gross primary production (GPP) and ecosystem respiration (R eco ) projected for Kobbefjord and Zackenberg using the SPA model forced with HIRHAM5 climate data.The solid line represents the annual mean C flux, and the error bars account for the sum of the uncertainties propagated by 1) the sensitivity analysis averaging the ±10 % change in all vegetation parameters used in SPA and 2) the two different RCP emission scenarios from HIRHAM5.

Fig. 6 .
Fig.6.Ranked sensitivity indices influence the NEE variability between the baseline (grey) and future periods (brown and green for Kobbefjord and Zackenberg, respectively).The presented variability is subject to the average ±10 % change in each vegetation parameter in the SPA model.
Fig.7.Present and future (based on RCP8.5 scenario) annual mean C sink strength from Kobbefjord (brown) and Zackenberg (green) comparing the control setups using the nominal parameter set and two experimental setups including enhanced and weakened foliar N inputs from the opposite site.Similar results using the RCP4.5 scenario are displayed in Fig.S4.
López-Blanco et al. (2020)estimates calculated in this study.ofasimulatedtundraheath site in West Greenland was similarly strengthened by a 2 °C warming experiment than by the joint warminginduced increase in high N availability.The observation that NEE and GPP were not affected by high lateral N input as an isolated factor somewhat agrees with the fact that temperature (climate) will play a central role.There is no doubt that local N conditions are crucial for explaining the variability observed both in the field (Table1;López-Blanco et al. (2020)) and in SPA simulations (Figs. sink