Reconciling the Carbon Balance of Northern Sweden Through Integration of Observations and Modelling

The boreal biome plays an important role in the global carbon cycle. However, current estimates of its sink‐source strength and responses to changes in climate are primarily derived from models and thus remain uncertain. A major challenge is the validation of these models at a regional scale since empirical flux estimates are typically confined to ecosystem or continental scales. The Integrated Carbon Observation System (ICOS)‐Svartberget atmospheric station (SVB) provides observations including tall tower eddy covariance (EC) and atmospheric concentration measurements that can contribute to such validation in Northern Sweden. Thus, the overall aim of this study was to quantify the carbon balance in Northern Sweden region by integrating land‐atmosphere fluxes and atmospheric carbon dioxide (CO2) concentrations. There were three specific objectives. First, to compare flux estimates from four models (VPRM, LPJ‐GUESS, ORCHIDEE, and SiBCASA) to tall tower EC measurements at SVB during the years 2016–2018. Second to assess the fluxes' impact on atmospheric CO2 concentrations using a regional transport model. Third, to assess the impact of the drought in 2018. The comparison of estimated concentrations with ICOS observations helped the evaluation of the models' regional scale performance. Both the simulations and observations indicate there were similar reductions in the net CO2 uptake during drought. All the models (except for SiBCASA) and observations indicated the region was a net carbon sink during the 3‐year study period. Our study highlights a need to improve vegetation models through comparisons with empirical data and demonstrate the ICOS network's potential utility for constraining CO2 fluxes in the region.

variable (Smith et al., 2020). There are urgent needs to address these uncertainties and rigorously characterize the variations, for several reasons. First, it is well known that the increasing atmospheric CO 2 concentration (Mittler & Blumwald, 2010;Pachauri et al., 2014) has resulted in global warming and unprecedented climatic changes (Barber et al., 2000;Serreze et al., 2000). Moreover, "Arctic amplification" has contributed to dramatic melting of Arctic Sea ice and spring snow cover (Cohen et al., 2014), but clarification of the extent to which these climatic changes will modify the boreal region's role in the global carbon cycle is needed (Bradshaw & Warkentin, 2015;Pan et al., 2011;Tagesson et al., 2020). Second, the Nordic countries are striving to develop more sustainable low-carbon societies (Urban & Nordensvärd, 2018), and Sweden in particular aims to become the first "fossil (fuel)-free welfare nation" in the world, with zero net emissions of greenhouse gases by 2045 (Xylia & Silveira, 2016). Efforts to achieve this ambitious goal would be helped by thoroughly verified negative emissions that could offset some anthropogenic greenhouse gas emissions (Fauré et al., 2019;Rodriguez et al., 2020). Thus, it is essential to obtain robust and precise estimates of the national carbon balance including the sink strength of the boreal forest region.
Climate change and the associated temperature and precipitation trends trigger both positive and negative carbon cycling feedbacks (Braswell et al., 1997;McGuire et al., 2009). Understanding of biosphere-climate feedback at a regional scale is necessary to evaluate the boreal forest ecosystems' potential to support efforts in mitigating global climate change (Kondo et al., 2020). However, current understanding of these processes at a regional scale primarily relies on terrestrial biosphere models (Cramer et al., 2001;Friend et al., 2007;Haas et al., 2013). Hence, comprehensive evaluations of these models with observations covering multiple spatiotemporal scales is required to assess the robustness of their estimates of biospheric responses to changes in meteorological variables (Baldocchi et al., 2001;Jung et al., 2007;Mahecha et al., 2010;Moorcroft, 2006;Siqueira et al., 2006).
A key variable to consider is Net Ecosystem CO 2 Exchange (NEE). This is the net flux resulting from the major vertical CO 2 fluxes between an ecosystem and the atmosphere (Baldocchi, 2003) and a primary measure of the ecosystem's carbon sink-source strength (Kramer et al., 2002). It is generally determined by subtracting fluxes into the atmosphere generated by soil and plant respiration from fluxes in the other direction generated by photosynthesis. Other components (such as lateral transfers of carbon in and out of the ecosystem; emission of volatile organic carbon, methane, and carbon monoxide; release of soot and CO 2 from fire; and fluxes from terrestrial inland water bodies) also contribute to the "Net Ecosystem Carbon Balance"(NECB) (Chapin et al., 2006). Available models account for these flux components in various ways. Here, the term NEE refers to the net ecosystem land-atmosphere carbon exchange, defined as the net carbon flux resulting from gross primary production (GPP, regarded as the total carbon flux from atmosphere to biosphere) and ecosystem respiration (Reco, regarded as the total carbon flux from biosphere to atmosphere). Discrepancies in estimates of the NEE derived using different methods limit the understanding of interactions between Nordic ecosystems and atmospheric CO 2, as well as these ecosystems' likely responses to future climate change and associated meteorological events.
Current estimates of the various flux components are generated by one of three approaches: ground-based field measurements (Riederer et al., 2014;Ramonet et al., 2020), bottom-up (vegetation/ecosystem models) modeling (Peng, 2000;Bastos et al., 2020), and top-down (atmospheric inversion) modeling (Gurney et al., 2002;Monteil et al., 2020). In field studies, eddy covariance (EC) techniques are often used to measure exchange fluxes directly and estimate seasonal dynamics and interannual variations of NEE, GPP, and Reco (Law et al., 2000Stoy et al., 2006;Schmidt et al., 2012). The establishment of regional and global networks of ground-based micrometeorological EC flux towers has provided continuous and instantaneous data on the atmosphere-biosphere exchange of water, energy, and carbon via fluxes to and from all types of ecosystems (Baldocchi et al., 2001;Papale et al., 2006). When sensors are mounted on tall towers (at appropriate heights), EC measurements can provide indications of landscape-scale fluxes spanning several tens of km 2 (Chi et al., 2019(Chi et al., , 2020. In addition, terrestrial vegetation/ecosystem models have become powerful tools for filling measurement gaps, extrapolating observed trends to responses in anticipated future conditions, predicting climate change impacts, understanding processes in broader contexts, and retrospective analyses of ecosystem responses to previous climatic changes (Amthor et al., 2001). They can be used to estimate the growth of vegetation and changes in carbon storage between different carbon pools (leaves, wood, roots, and soils) at various spatial and temporal scales (Pan et al., 2014). Ecosystem models vary in complexity of structures, parameters, and processes that control photosynthetic carbon uptake and respiratory release. Thus, the models' estimates of terrestrial primary production depend on the underlying assumptions, their structure, and formulation. They range from purely empirical models to process-based ecosystem models that can be used to investigate responses of ecosystem processes to both environmental and anthropogenic factors at regional and global scales (Melillo et al., 1993;Tian et al., 2010). While these bottom-up models can provide detailed estimates of the fluxes and characterization of their spatial and temporal variability that can be applied in future climate simulations, large uncertainties may be propagated into the cumulative fluxes (Sitch et al., 2015).
Carbon fluxes can be highly variable even within a single biome (Schulze et al., 1999), so validating the model's estimate at one location does not necessarily guarantee its validity at a large scale. Since the land-atmosphere carbon fluxes clearly affect the atmospheric CO 2 content, they can also be estimated from observations of atmospheric CO 2 concentration. For this, an atmospheric transport model that computes the concentrations as functions of fluxes is needed. Because of atmospheric mixing, fluxes at one location influence CO 2 concentrations over a large area and conversely a single CO 2 observation is sensitive to flux changes over a large area. However, estimates from these approaches sometimes significantly differ, and the discrepancies are more pronounced for boreal forests (where carbon fluxes are relatively small) than in many other ecosystems (Campioli et al., 2016). Therefore, to improve understanding of the boreal biosphere's role in the climate system, it is essential to reconcile the carbon flux estimates from ecosystem to global scales (Baldocchi et al., 2001;Ballantyne et al., 2021;Running et al., 1999).
Valuable data for such purposes are provided by the ICOS SVB (Integrated Carbon Observation System-Svartberget) atmospheric site, including tall tower EC and atmospheric concentration measurements that collectively facilitate novel model-data fusion studies covering a homogenous boreal landscape that is very similar to landscapes in many parts of Northern Sweden (Lindström et al., 2002;Laudon et al., 2021). Another five ICOS sites across the Nordic region provide CO 2 concentration measurements that can be compared to modeled concentrations from an atmospheric transport model, driven by an estimate of the regional carbon fluxes. In this Nordic region, there was a severe drought in 2018, which was even stronger in some places than the drought that occurred in central Europe in 2003 . The aim of this study is to reconcile the carbon balance of Northern Sweden by integrating fluxes from models and observations over Scandinavia at site to regional scales. There were three specific objectives. First, to compare model estimates of the fluxes and concentrations with site-scale observations during the years 2016-2018. Second, to determine how well the atmospheric observation network can constrain the fluxes in the Nordic region. Third, to compare modeled and observed carbon fluxes and atmospheric concentrations during the drought in 2018 (one of the three warmest years on record for Europe; WMO, 2019, Figure S1 in Supporting Information S1) at the ICOS SVB site and the Northern Sweden region.

Models and Data
We used flux estimates from four gridded vegetation models. Two are dynamic global vegetation (DGVM) models: the ORCHIDEE (Organizing Carbon and Hydrology In Dynamic Ecosystems) and LPJ-GUESS (Lund-Potsdam-Jena General Ecosystem Simulator) models. The others are diagnostic models: the VPRM (Vegetation Photosynthesis and Respiration) and SiBCASA (Simple Biosphere/Carnegie-Ames-Stanford Approach) models. We assessed their performance by comparing results they yielded with empirical data at two spatial scales, as follows.
In the first approach, we compared modeled and observed fluxes of NEE and its components Reco and GPP at site level. The selected modeled fluxes are those generated for the nearest model grid cell to the SVB observation site (64. 25°N, 19.77°E) during the study period. The observed fluxes are from the tall tower EC data set, which are typically representative of the net flux at landscape scale (Chi et al., 2019). These direct comparisons reveal the vegetation models' ability to reproduce observed fluxes at site scale. Fluxes yielded by the four models for the Northern Sweden region are also compared.
In the second approach, we compared atmospheric concentration measurements at a large scale, which requires forward simulations from an atmospheric transport model. In this study, we used the FLEXPART 10.4 (Flexible Particle dispersion model) Lagrangian transport model (Pisso et al., 2019) to compute flux transport in a regional domain ranging from 0°W, 54°N to 20°E, 70°N ( Figure 1). In addition, background concentrations (boundary condition) were taken from a global TM5 inversion, following the LUMIA (Lund University Modular Inversion Algorithm) setup described by Monteil and Scholze, 2021. For the sake of readability, we refer to this transport model setup as LUMIA in the rest of the paper. LUMIA computes the CO 2 concentrations at the six ICOS sites in the domain (Figure 1) using estimates of the surface CO 2 fluxes (NEE and anthropogenic emissions), gridded on a 0.5° × 0.5°, 3-hourly resolution (fluxes were regridded by aggregation or re-binning when necessary). Monthly anthropogenic emissions are based on a pre-release of the EDGAR4.3 inventory (Janssens-Maenhout et al., 2019) for the base year 2010 (doi: 10.18160/Y9QV-S113). At all six sites, continuous CO 2 observations are made at several vertical levels. Only the upper level at each site was computed in LUMIA, as these are where the accuracy of the transport model is expected to be the highest.
The in-situ observations of fluxes and concentrations are described in detail in Section 2.1. Table 1 provides a summary of the observations of fluxes and concentrations and Table 2 provides the model estimates, their sources, original resolution, and geographical position of the sites. Section 2.2 details the structure and functioning of the models and differences between them, while Section 2.3 provides a detailed description of the LUMIA atmospheric transport model. The study period covers the years 2016-2018 (inclusive).  Vindeln, in the Krycklan catchment, Northern Sweden (Chi et al., 2020). In addition to the ICOS ecosystem-level EC system (mounted at 34.5 m), an additional EC system was installed on the same 150 m tall tower to measure CO 2 , water vapor, and energy exchange between the boreal landscape and atmosphere. This EC system was mounted at 70 m height above the ground from March 2016 until October 2017, then lowered to a measurement height of 60 m. These tall tower EC measurements are typically representative of the net fluxes at landscape-scale (covering ca. 3-4 km around the tower, giving a ∼25-50 km 2 footprint for the 60 m measurement level). Hereafter, we refer to data provided by this landscape-level EC system as "EC data". The EC measurement set-up consists of a CSAT3 3D sonic anemometer and EC155 closed-path gas analyzer (Campbell Scientific Inc., USA). Further information about the data processing, corrections, quality checks and controls, and final NEE estimations are reported in detail by Chi et al. (2019). The set of NEE measurements used here was gap-filled and source-partitioned into the two flux components, GPP and Reco, using the R-package REddyProc (v3.6.3, Wutzler et al., 2018). All data processing steps closely followed the ICOS protocol, but it should be noted that the EC system is not part of the ICOS network.

Atmospheric CO 2 Concentration Measurements
We collected continuous in situ observations of atmospheric CO 2 concentration from ICOS tall towers (ICOS, 2019) (https://meta.icos-cp.eu/collections/4H3RS8YtXIt_WTcjrSSkaQ-O; doi: 10.18160/RHKC-VP22). These included observations from six sites in our regional study domain, for 2017 and 2018 with data availability changing from site to site within that period. Data for 2016 are available at some stations that are not considered for maintaining the consistency across all stations. Most of these high frequency CO 2 concentration measurements were recorded at several vertical levels up to 150 m using tall towers that are now included in the European ICOS atmospheric network. We selected the observations at the highest level from each site. All ICOS sites in the domain where concentration measurements used in the study were marked in Figure 1. Hereafter the term "observations" refers to the measured CO 2 concentrations, unless otherwise specified.

LPJ-GUESS (Lund-Potsdam-Jena General Ecosystem Simulator)
LPJ-GUESS (Smith et al., 2014;Smith, 2001) is a process-based terrestrial vegetation model that can predict structural, compositional, and functional properties (including biomass and height) of native ecosystems at landscape, regional, and global scales. It is a dynamic global vegetation model that includes detailed description of the biochemistry (unlike diagnostic models) of the focal vegetation and representations of the energy absorption, carbon, and nitrogen cycling. The current version of LPJ-GUESS includes an interactive plant and soil nitrogen cycle and associated constraints on plant production and ecosystem carbon balance. Photosynthesis is modeled with a simplified Farquhar scheme (Collatz et al., 1991) and vegetation growth is simulated dynamically in a so-called cohort-mode representing different tree age classes in a grid cell. LPJ-GUESS simulations are driven by the 6-hourly Climatic Research Unit Japanese Reanalysis (CRU-JRA) meteorological data set (Weedon et al., 2014) as driving data. Fast (plant physiology) processes are implemented on the time step of the meteorological input (air temperature, precipitation, and radiation) data and slow processes ( LPJG_Conif LUMIA transport simulations FLEXPART regional transport model 0.5° × 0.5°LU, Sweden

Table 2 Details of Model Simulations (Including Both Gridded and Site-Specific LPJG Simulations) Used in the Study
two site-specific NEE simulations generated by the LPJ-GUESS model forced with local site meteorological data (surface air temperature, precipitation and incoming shortwave radiation measured at the site for the years 2015-2018) in the time series analysis at site scale ( Figure 2). For the other analysis presented here, we used flux estimates from the LPJ-GUESS gridded simulations. The site-level runs were carried out without site-specific calibration, using inputs of nitrogen deposition and soil properties from global datasets (Smith et al., 2014). One (denoted by LPJG_PNV) involved use of the PNV (Potential Natural Vegetation) vegetation type that establishes under the model's dynamic vegetation processes (at SVB: boreal evergreen coniferous trees, understory grass and some boreal summer-green trees). In the other (designated LPJG_Conif), the plant functional type that could establish was restricted to match the SVB site vegetation cover (boreal evergreen coniferous trees and understory grass), but neither site management history nor disturbances were considered. Finally, no wetland features were activated in the site simulations.

ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosystems)
ORCHIDEE is a global process-based terrestrial biosphere model developed as a land surface component of Institut Pierre-Simon Laplace (IPSL)-CM5 earth system model (Krinner et al., 2005 and http://www.ipsl.jussieu. fr/∼ssipsl/). It computes carbon, water, and energy fluxes between the land surface and atmosphere, as well as within the soil-plant continuum. The model computes GPP with CO 2 assimilation based on functions presented by Farquhar et al., (1980) for C3 plants. Land cover changes (including deforestation, regrowth, and cropland dynamics) were prescribed using annual land cover maps derived from the harmonized land use data set (Hurtt et al., 2011) in combination with the European Space Agency-Climate Change Initiative (ESA-CCI) land cover products. The carbon module estimates: growth and maintenance respiration, which sum up to autotrophic 7 of 20 respiration (Ruimy et al., 1996); Net Primary Production (NPP); heterotrophic respiration (Rh); and NEE (NEE = NPP-Rh) at 3 hourly temporal and 0.125° spatial scale using the CRU-HARMONIE meteorological input data (Petrescu et al., 2021). ORCHIDEE has been shown to simulate responses of natural ecosystems in Europe to climatic anomalies, such as the 2003 heat wave, reasonably well (Ciais et al., 2005).

VPRM (Vegetation Photosynthesis and Respiration Model)
VPRM (Mahadevan et al., 2008), a diagnostic model, is a satellite-based assimilation scheme that estimates hourly NEE values based on a light-use efficiency approach and temperature-dependent estimates of ecosystem respiration using high-resolution data. It has an extremely simple mathematical structure, with minimal parameters. VPRM uses ECMWF (European Centre for Medium-Range Weather Forecast) operational meteorological data for radiation and temperature inputs, the SYNMAP land cover classification scheme (Jung et al., 2006), together with EVI (Enhanced Vegetation Index) and LSWI (Land Surface Water Index) values derived from Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance data. Model parameters were optimized for Europe using EC measurements collected during 2007 from 47 sites (Kountouris et al., 2015). The VPRM simulations used here has been produced at 0.5° spatial and hourly temporal resolution. VPRM also provides partitioning of NEE into its two flux components.

SiBCASA (Simple Biosphere/Carnegie-Ames-Stanford Approach)
SiBCASA (Schaefer et al., 2008;Van der Velde et al., 2014) combines parameterization of the Simple Biosphere model (SiB) with the biogeochemistry of the Carnegie-Ames Stanford Approach (CASA) and additional features. It can estimate terrestrial carbon fluxes from diurnal to decadal timescales by calculating water, carbon, and energy exchanges between 25 soil layers, plants, and the atmosphere. It derives photosynthesis rate using the Ball-Berry-Woodrow model of stomatal conductance (Ball et al., 1987), and C3 and C4 vegetation types are treated separately in the kinetic enzyme model of Farquhar et al. (1980). The simulation used here is forced using meteorological inputs from ERA-Interim, and run with 10-min time steps, 1 × 1 degrees spatial resolution, and hourly temporal resolution. SiBCASA uses observed LAI (Leaf Area Index) measurements to guide semi-prognostic canopy development; and in this study, we used the particular SiBCASA simulation focused on the 2018 drought.

Atmospheric Transport Simulations-LUMIA
CO 2 is a stable tracer in the atmosphere: its sources and sinks are almost exclusively at the surface (small quantities are produced directly in the atmosphere by CO oxidation and emitted by aircraft, but these are negligible in comparison to CO 2 exchanges at the surface). Thus, the atmospheric distribution of CO 2 largely reflects the strengths and spatio-temporal distributions of surface CO 2 fluxes and subsequent mixing (convective and advective) of air. Thus, provided that the atmospheric mixing can be properly simulated, it is possible to relate surface CO 2 fluxes to observed CO 2 concentrations.
We implemented a regional transport model, following the approach used in the LUMIA system (Monteil & Scholze, 2021). The modeled mixing ratios y m are defined as the sum of two contributions: a "foreground" (near-field) contribution, which is the impact of fluxes on CO 2 concentrations in the direct vicinity, spatial and temporal, of the observations, and a "background" (far-field) contribution, which is the contribution of fluxes further away in time and/or space. In LUMIA, the background contribution is computed using the global, coarse resolution (6° × 4°) transport model TM5, and the foreground using observation footprints of the FLEXPART Lagrangian particle dispersion model. The observation footprints were computed using the FLEXPART 10.0 model, which simulates the dispersion backward in time of numerous virtual air particles (10,000 in our setup), based on meteorological input data. The aggregated residence time of the particles in a layer below 100 m above ground is taken as a proxy for the observations' sensitivity to the surface CO 2 exchanges. The footprints were computed at half degree spatial and three-hourly temporal resolution, and the FLEXPART model was driven by ECMWF ERA-Interim reanalysis data (Pisso et al., 2019). The footprints cover a regional domain ranging from 0°E, 54°N to 30°E, 70°N (Figure 1), and each covers a period up to 14 days before the observation time (although, in most cases, the sensitivity drops to zero in just a few days, as the particles all exit the regional domain).
The background concentrations ( bg E y ) were computed using the global, coarse resolution TM5-4DVAR inverse model, following the approach described in further detail by Monteil and Scholze (2021). Briefly, TM5 not only provides a boundary condition at the edge of the domain, but also computes the transport of that boundary condition toward the observation sites. There is therefore no exchange of masses between the two transport models; each calculates a partial mixing ratio. The TM5 simulation was driven at global, 6° × 4°, 3-hourly resolution, by ECMWF ERA-Interim meteorological data and transporting CO 2 fluxes. The latter were obtained from global TM5-4DVAR inversion assimilating flask observations from the NOAA ESRL Carbon Cycle Cooperative Air Sampling Network (Dlugokencky et al., 2013) and, within Europe, from the ICOS network of tall towers, in efforts to ensure that the background concentrations were as unbiased as possible. See Monteil and Scholze (2021) for more complete details.

Results
We compared the carbon flux estimates of the Northern Sweden region using vegetation models at regional scale and empirical observations of fluxes and concentration at site scales. Section 3.1 summarizes results of the comparison of fluxes from vegetation models with tall tower EC measurements to illustrate how well the vegetation models reproduced the observed ecosystem fluxes at the study site. In Section 3.2, we compare the atmospheric CO 2 concentrations simulated by LUMIA using the terrestrial fluxes simulated by the vegetation models with the observed atmospheric concentrations to assess the applicability of site-scale conclusions to regional scale patterns and thus facilitate evaluation of ecosystem models' performance.

Analysis of the In Situ Data
Seasonal variations in weekly sums of NEE, GPP, and Reco simulated by the four vegetation models are compared to tall tower EC measurements at Svartberget during the 2016-2018 study period in Figure 2. Flux estimates obtained using the LPJG_Conif and LPJG_PNV (see Section 2.2.1) are shown in this comparison, but other LPJ-GUESS modeling results presented in this paper are based on gridded simulations. Results of the two LPJ-GUESS site simulations and gridded LPJ-GUESS simulations at site and region levels are presented separately in Figure S2 in Supporting Information S1. Briefly, the fluxes obtained from the site simulations are similar to the one obtained from the gridded simulation, but have larger variability on daily timescales. We suspect this is because the landscape around the Svartberget site is fairly homogeneous. However, the locally measured meteorological parameters are more variable than the corresponding parameters in the gridded CRU-JRA meteo-data. Interestingly, the net fluxes from the gridded simulation are closer to the observed fluxes, possibly because the site simulations were not calibrated to the specific site conditions as explained above. According to both in the EC data and models, CO 2 was emitted (net release of CO 2 occurred) from October until the end of the following April, regarded here as winter period. This offset uptake during the spring and summer, from March until the following September, with a maximum rate of 50-100 gC m −2 week −1 . Simulations by all models for all years showed the same general pattern, but emissions and uptake were stronger according to LPJ-GUESS and ORCHIDEE than the other models ( Figure 2b). Modeled fluxes from VPRM appeared to be underestimated with lower CO 2 release in winter. Both VPRM and SiBCASA showed lower net CO 2 uptake in the summer of 2016, and it persisted for most of 2017 and 2018 for SiBCASA compared to the tall tower EC measurements. The differences between modeled and observed NEE were larger during winter periods (weekly model-data mismatch 20 gC m −2 ), than during summer (weekly model-data mismatch 10 gC m −2 ), and largest in the VPRM and LPJ-GUESS simulations. While VPRM simulated zero fluxes throughout the winter, fluxes were only non-zero in LPJ-GUESS simulations in the beginning of each year (January to March). All site-scale LPJ-GUESS simulations show an emission peak in August 2018 with a mean model-data mismatch of 30 gC m −2 with respect to the EC data. This summer emission peak in LPJ-GUESS simulations, which conflicts with simulations by the other models and the EC-based NEE observations, is largely due to substantially higher Reco values. In contrast, the EC data show a sharp increase in net CO 2 uptake, compared to the seasonal mean, in June 2018, due to an increase in GPP and decrease in Reco during this period (Figures 2b and 2c).
For the GPP and Reco component fluxes, differences between the simulations and EC data were more pronounced during the summer periods in the three study years. Like the LPJ-GUESS site-scale simulations, OR-CHIDEE overestimated GPP and Reco substantially during the summers. Essentially, the two process-based models (LPJ-GUESS and ORCHIDEE) overestimated both component fluxes with respect to the EC fluxes. Nevertheless simulated and recorded NEE fluxes are more similar. In stark contrast, VPRM underestimated the two component fluxes, but the simulated NEE matches the observed NEE well, except for zero winter fluxes. In addition, ORCHIDEE and SiBCASA simulated winter respiration and photosynthesis well, but not the other vegetation models. However, all vegetation models reproduced the observed fluxes and their seasonal dynamics within the model limitations during the entire annual period.
To explore the disparities in intra and interannual results of the applied models further, it is important to examine the net CO 2 exchange at landscape scale, the total annual sums and differences in spring onset as well as periods when gross CO 2 uptake via photosynthesis exceeds respiratory CO 2 release (or vice versa). These features can be explored with the help of cumulated curves of ecosystem fluxes (Figure 3) generated by the four vegetation models and tall tower EC measurements at SVB for the years 2016, 2017, and 2018. Both models and observations showed the same general seasonal patterns, which switch between winter net emission to summer net uptake in May and October, respectively. However, there are some notable differences among the simulations generated by the four models and the EC data set. The interannual difference in NEE curves generated by all four models started around May. This marks the onset of the spring period and corresponds to the start of increased CO 2 uptake. According to the EC data, this onset is different between years, occurring in March in 2016 but later in subsequent years (Figure 3e). The simulations do not display this difference, so the models apparently have limited ability to capture the interannual variation in the onset of the vegetation period manifested in the EC data. In addition, the cumulated NEE curve based on EC data showed a sharp shift upward in September 2017 (Figure 3e) due to an increase in Reco, which is captured by the VPRM and SiBCASA models, but not the others (Figures 3b and 3c). Meanwhile, the annual sums obtained from the EC data and simulations agreed well considering the representation error of the vegetation models even though their validity cannot be compared, because of the use of data pertaining to a grid point in a regional simulation in the analysis. The summer and annual net CO 2 uptake were both lowest in 2017 (Figures 3b-3e) rather than during the drought year 2018 according to the EC data and most of the models, and highest in 2016, with VPRM yielding higher annual sums for all years than the other models. However, the LPJ-GUESS simulation suggested that annual net CO 2 uptake was lowest in 2018.
The seasonal and annual patterns of NEE, GPP, and Reco are more precisely captured (in terms of similarity to the EC data) by SiBCASA and ORCHIDEE (Figures 2 and 3). The magnitudes of daily NEE fluxes and annual sums simulated by LPJ-GUESS and VPRM are considerably higher and lower, respectively. The ecosystem model LPJ-GUESS and diagnostic model VPRM did not simulate the respiratory loss during winter indicated by the EC data. However, LPJ-GUESS generally reproduced the summer NEE uptake shown in the EC data well, except for the 2018 summer. The VPRM simulation did not capture the distinctive seasonal winter emission to summer uptake transition in any of the three years referring to the observations ( Figure 3b). As shown in Figure 2, the component fluxes were respectively over-and underestimated by LPJ-GUESS and VPRM, but the balance between the two flux components yielded by LPJ-GUESS still resulted in comparable NEE fluxes to the EC estimates. The annual simulated NEE of 2017 differed from the EC data, but the modeled component fluxes generally agreed well, suggesting that GPP was higher and Reco lower in 2017 than in the other years. It is mentioned relative to the year-to-year difference, not referring to the absolute flux magnitude. In 2016, the GPP did not significantly differ from the GPP in other years and Reco was lower (which may account for the early transition to uptake in 2016). The daily NEE simulated by all models correlated well with the corresponding EC observations (r 2 ≥ 0.73, RMSE ≤ 1.22 gC m −2 day −1 ) with the lowest linear Pearson correlation and highest RMSE for LPJ-GUESS (r 2 = 0.73, RMSE = 1.22 gC m −2 day −1 ) and highest correlation with lowest RMSE for VPRM and SiBCASA (r 2 = 0.84, RMSE = 1.09, and 1.04 gC m −2 day −1 ) ( Table S1 in Supporting Information S1). Overall, the SiBCASA simulation provided the highest correlation with daily values of the observed EC fluxes (NEE, GPP and Reco) and lowest RMSE (r 2 > 0.84 and RMSE < 1.14 gC m −2 day −1 ; Tables S1 and S2 in Supporting Information S1) over the three study years. The annual GPP sums generated by the four models differ (Figure 3), but the linear correlation coefficients between modeled and EC estimates of daily GPP values are >0.88 with RMSE, which is <2.47 gC m −2 day −1 . The modeled Reco values are less consistent with the observations (correlation coefficients ranging from 0.72 to 0.79 and RMSE 0.92-3.01 gC m −2 day −1 ). In addition, annual GPP was lower in 2017 and annual Reco higher in 2018 than in the other years (Figures 3j and 3o) according to both the EC data and all models (Figures 3f-3i and 3k-3n).

Regional-Scale Analysis
In this section, we assess the potential utility of atmospheric CO 2 observations for validating bottom-up estimates of CO 2 fluxes at regional scale, and whether the site-scale features observed at SVB are representative of CO2 fluxes in Northern Sweden. We focus specifically on the "Northern Sweden" region (within the blue solid line in Figure 6), covering 320,000 km 2 around SVB. This region is regarded as fairly typical for much of the rest of Northern Sweden, and the measured concentrations at SVB are most sensitive to fluxes in this region. In Section 3.2.1, we considered the fit of results generated by the atmospheric transport model described in Section 2.3, driven by modeled CO 2 flux estimates, to CO 2 concentrations measured at ICOS sites in the Nordic study domain. It should be noted that the background concentration is always the major component of the total measured concentration at any location in the domain. The input ecosystem-driven fluxes only contribute significantly to changes in the relatively minor foreground concentration.

Atmospheric CO 2 Concentration
The magnitude of weekly variation in daytime (11:00-14:00 UTC) CO 2 concentrations at SVB (Figure 4) was similar according to both the simulations and observations. We focused on daytime values because the transport model is expected to be most accurate when the planetary boundary layer (PBL) is well developed. We considered weekly means partly for consistency across the flux and concentration analysis and partly because the weekly averages exhibit the seasonal pattern more clearly with less noise compared to the daily averages. Hence, they give a better appreciation of the overall model performance whereas shorter-time averages are relatively more sensitive to transport model errors and of the fine-scale performance of the vegetation models, especially for the 3-year study period (as illustrated by daily concentration in Figure S3 in Supporting Information S1). The Pearson correlation coefficients between observed and simulated CO 2 concentrations always exceeded 0.93 with RMSE below 3.42 ppm for daytime hours at SVB (Table S3 in Supporting Information S1). However, the correlation coefficients were ca. 0.2 lower for the winters than for the summers. The high Pearson correlations between the simulations by all models and observations confirm the quality of the transport model and background concentration data, especially during daytime and summer.
The observed and modeled concentrations were both lower than the background concentration during summer, owing to CO 2 uptake, and vice versa during winter due to CO 2 release. On average, over the entire study period (2016)(2017)(2018), VPRM-based estimates of CO 2 concentrations were closest to the observations, with less than ±5 ppm bias and higher than 0.92 Pearson correlation for all sites in the study domain. SiBCASA and ORCHIDEE occasionally over-/underestimated concentrations, relative to observations, by more than ±5 ppm sometimes. The unusual emission spike in the LPJ-GUESS simulation during August 2018 (Figure 2) at the SVB site was also manifested in the LPJ-GUESS-driven LUMIA simulation (2-4 ppm higher concentration than the observations during 2018 summer) while other LUMIA simulations closely followed the observed concentration throughout the study period. Apart from this, all models provided better fits to the observations during 2018 than the previous year, 2017. The near-zero wintertime fluxes in LPJ-GUESS and VPRM simulations led to underestimation of winter CO 2 concentration during the period 2017-09 to 2018-04 ( Figure 4). As the background concentration is the major component of the total measured concentration, it accounts for some of the concentration peaks, such as the one in November 2018. Weekly mean model-data mismatches (modeled minus observed concentrations) at all ICOS atmospheric sites in the study domain are shown in Figure 5. The weekly biases at SVB in simulations by all the models stayed within −10 to 6 ppm (Figure 5a), while at other sites, it varied between −10 and 10 ppm. On average across all sites, there was a weekly bias of 3 ppm in LPJ-GUESS and VPRM simulations, but in some cases (e.g., the summer of 2018 at SVB, PAL, and PUI) LPJ-GUESS estimates were up to 5 ppm higher than the observed concentrations. Estimates by the other models had a negative bias of up to −8 ppm at SVB. ORCHIDEE and SiBCASA simulations had a more positive bias than the others throughout the study period and overestimated the weekly concentrations by up to 10 ppm in winter. It should also be noted that concentration estimates across all the sites had higher bias in winter. This comparison of the transport models' correlations with atmospheric observations suggests that flux estimates by the diagnostic models, especially VPRM, are the most realistic at the regional scale (The correlations are shown in Table S3 in Supporting Information S1).

Regional CO 2 Fluxes
We investigated the strength and spatial distribution of annual NEE over the study domain to explore the representativeness of flux estimates at SVB for fluxes in the entire Northern Sweden region. Figure 6 illustrates the annual NEE estimated by the four vegetation models for the drought year 2018 over the Nordic domain. As already noted, VPRM indicated that both the site and the entire study region were larger sinks in 2018 than the other models. All of the models except ORCHIDEE agreed that the southern part of the domain was a larger sink in 2018 than the northern part. Although there were differences in the annual NEE simulated by the models, the distribution of NEE in each simulation within the Northern Sweden region was homogenous.
The relative differences in annual NEE at the site (landscape) and regional (Northern Sweden) scales for each of the three study years are shown in Figure 7. All models except SiBCASA suggested that both the site and region were a carbon sink during the three years, although the estimated sink strength varied from −275 to −25 gC m −2 year −1 . The site and region were smaller sinks according to all the models and EC data (or even a small source according to the SiBCASA simulation) in 2017 than in 2018 and 2016, except the LPJ-GUESS simulation, in which the sink was smallest in 2018. NEE budgets obtained from ORCHIDEE and VPRM were quite constant from year to year with standard deviation ±6.6 and ±19.3 gC m −2 year −1 , respectively. The substantial discrepancy in annual NEE for 2018 at site scale between the LPJ-GUESS and other simulations resulted from the summer emission peak it generated at SVB site ( Figure 2). However, this site-scale feature was not prominent at the regional scale (see also Figure S2 in Supporting Information S1), and the annual NEE in LPJ-GUESS simulations is very similar to that of ORCHIDEE and SiBCASA simulations for 2018 at the regional scale. Except for this peak in the LPJ-GUESS simulations, the NEE dynamics and annual budgets at the SVB site are similar to those for the rest of the Northern Sweden domain in the simulations by all models and the model-to-model differences are larger than the year-to-year differences within each simulation.
We included all available hourly daytime observations when calculating the weekly averages at SVB presented in Figure 4a. However, CO 2 observations at SVB are not equally sensitive to fluxes from within the Northern Sweden region. Depending on the wind intensity and direction, they can be more sensitive to other parts of the domain or the background concentration. Thus, we used the FLEXPART footprints (see Section 2.3) to filter  out the observations that were least sensitive to fluxes in the Northern Sweden region. For each observation, we computed a "regional sensitivity index": We arbitrarily used a threshold value of FP cum <= 70% rejecting observations, set as a compromise between selecting observations with close to maxima local influence and not rejecting too many observations. The interpretation is therefore mostly qualitative. We show weekly averages of the modeled and observed concentrations in Figure 4b, for the observations that met this criterion. The influence of regional fluxes was most prominent during some peaks, particularly in August 2017 and March 2018.
To determine how a network of observation sites could help us to further constrain the fluxes in Northern Sweden, we repeated the same "regional sensitivity index" calculation for the other ICOS sites in our domain. Figure 8 shows time courses of that fractional sensitivity metric for all the sites. Apart from SVB, the sites that were most sensitive to fluxes in the region are Pallas (PAL) and Norunda (NOR), and they were particularly sensitive during the times of the unusually high emissions in the LPJ-GUESS simulation at SVB in August 2018 and August 2017. However, these observations were also very sensitive to fluxes from their own region and other regions, which are uncertain. This complicates use of these observations for direct validation of the fluxes, but indicates that they could potentially provide important constraints in a formal inverse modeling approach.

Discussion
The carbon balance of boreal regions is of particular interest due to their high sensitivity to warming and high rates of projected warming (Ito et al., 2020;Vincent, 2020). Thus, we have explored in detail the agreement between CO 2 fluxes (NEE and its component fluxes GPP and Reco) obtained from vegetation models and observations, including their seasonal patterns and interannual variability, over Northern Sweden. We compared gross and net fluxes obtained using four well-established vegetation models (LPJ-GUESS, ORCHIDEE, VPRM and SiBCASA) to constrain the ecosystem carbon budgets at both the site and regional scale. All models unequivocally indicate that the region was a carbon sink during the 3-year study period (2016)(2017)(2018). Here, we focus on three aspects of the study: the difference in estimates among models and observation, the potential utility of ICOS observations for constraining fluxes in our domain, and the accuracy of the models' predictions for fluxes during the drought year 2018.

Disparities Among Estimates of Carbon Exchange in Northern Sweden Obtained From the Vegetation Models and Observations
There are substantial differences in the annual budgets estimated by the four models, largely due to differences in their internal structures and dynamics. The relative contributions of uptake and emission in the simulations differ from the observations at landscape scale, resulting in differences of 100-200 gC m −2 year −1 in NEE annual budgets. The VPRM and LPJ-GUESS models simulated the highest and second highest annual ecosystem carbon budgets generally. High carbon uptake is a known feature of VPRM, resulting from the model parameters' calibration against EC observations (Oney et al., 2017). The near-zero winter flux simulated by VPRM is not consistent with either the other models or tall tower EC fluxes. Since photosynthetic rates are extremely low in winter (so GPP is very close to zero), ecosystem respiration (Reco) continues and NEE is positive. The winter fluxes are more realistically represented by SiBCASA and ORCHIDEE, and the fluxes simulated by these models are closest to the empirically determined EC fluxes. This holds true for the component fluxes also from SiBCA-SA which compares best with EC fluxes for the whole study period despite estimating the lowest annual NEE budget among other models. In high-latitude forests, carbon uptake may generally be promoted more strongly by the higher radiation levels associated with drought than it is impaired by the decline in water availability (Flach et al., 2018), unless the drought is extremely severe. Our findings suggest that the considered model simulations of NEE fluxes differ significantly from those manifested in the EC data, although both the interannual and seasonal variations were quite similar. On the other hand, the EC flux components, GPP and Reco, are also estimates generated by flux partitioning algorithms from direct measures of NEE obtained using the EC technique. Thus, use of a source partitioning model in the prediction of fluxes may also contribute to uncertainties. We found that the sharp reduction in net CO 2 uptake in September 2017 recorded only in EC observations was caused by a rainy period (with 80 mm precipitation between 7 and 13th September 2017), in which global radiation (and hence GPP) was considerably reduced, while the rewetting likely increased Reco (see, Figure 3o). Since we only used the gridded simulations here with meteorological inputs from a regional-scale product rather than local observations of meteorological variables at SVB, it is not surprising that the models did not simulate any responses to such short-term weather events, which got smoothed out.

Can We Use an Atmospheric Observation Network to Constrain Fluxes in Northern Sweden?
As mentioned, one of the objectives of the study was to assess ICOS observations' potential utility for constraining surface CO 2 fluxes in Northern Europe, the study domain. The density of the observation network is limited, but the region has relatively homogeneous ecosystems (with north-south climate gradients), and fewer major urban areas than other parts of Europe. Thus, it is a good candidate area for estimating carbon fluxes at large scales using an inverse approach. The relative homogeneity of the flux patterns simulated by the four vegetation models confirms that the density of the network is not a fundamental obstacle (unless there is more variability in the fluxes than the terrestrial biosphere models indicate). Another potential obstacle is that the contribution of regional fluxes to the observed concentrations could be within the uncertainty of the background concentrations. However, the CO 2 transport simulations presented in Section 3.2 show that the atmospheric observations have good sensitivity to regional signals. The differences between the flux estimates lead to differences in concentrations of the order of 10-20 ppm ( Figure 5), which is significantly larger than the typical uncertainty of the observations and background concentrations (<1 ppm). The estimated regional sensitivities in Section 3.2.2 ( Figure 8) corroborate our conclusion that the observation network provides exploitable constraints, and we will implement a formal inversion approach in further research.

Performance of Vegetation Models During the 2018 Summer Drought
A negative effect of the 2018 drought on the NEE of Northern Sweden was observed, but it was less pronounced than its effect in the south of Sweden and other parts of Europe (Gharun et al., 2020;Rödenbeck et al., 2020;Peters et al., 2020;Lindroth et al., 2020). However, the NEE response to drought simulated by the models differed, with LPJ-GUESS simulation showing a net reduction in NEE uptake for 2018 at the SVB site, but not the observations or other models. The EC data showed an increase in net uptake in the early summer of 2018, which led to an annual budget similar to that of the "normal" reference year 2016, in which precipitation was near-normal across most of Europe . The observed annual carbon sink was lower in 2017, because in that year a late cold spring was followed by a cold and cloudy summer, resulting in lower than normal CO 2 uptake (Chi et al., 2019). Apart from the increased uptake shown in EC data and high emissions in LPJ-GUESS simulations in the summer of 2018, the general trend in the simulations indicated that annual net CO 2 uptake was lower in 2018 than the reference year except SiBCASA. The lower annual uptake reflects the limited impact of the dry summer, which was partly compensated by the colder start of the year and warmer spring (Gharun et al., 2020) than in 2016 at site and regional scale. Although 2018 was one of the warmest years on record Thompson et al., 2020; Figure S1 in Supporting Information S1), the drought did not have a strong impact on the annual carbon balance of Northern Sweden. Buras et al. (2019) obtained similar results and attributed the lesser impact of the 2018 drought to weaker sensitivity of the boreal mixed land-cover types in Northern Sweden. In addition, adverse effects of the 2018 drought were mediated more strongly by low soil water content (SWC) than high air temperature . Hence, it was likely to have relatively weak effects in the northern Nordic forests where the trees generally have relatively high rooting depths and access to deeper soil water resources than those in many other ecosystems (Teuling et al., 2010).
The spatial plots clearly show a sharp north-south gradient in ecosystem responses to the drought, as also reported in Bastos et al. (2020). The simulated NEE responses generally agreed well with previous reports based on NEE measurements , not only for Northern Sweden, but also for the whole of Northern Fennoscandia. However, the apparent discrepancies call for further efforts to improve model algorithms to enhance the accuracy of predicted carbon-cycle climate feedbacks under future climate change.

Conclusions
Knowledge of Nordic ecosystems' responses to climate change is crucial for assessing their potential contributions to mitigation efforts. To facilitate acquisition of such knowledge, we have established a framework to combine observations of fluxes and concentrations with data generated by vegetation and atmospheric transport models which enhance estimates of the carbon balance of Northern Sweden. We also assessed the SVB experimental site's representativeness for the region. We found that NEE estimates from four vegetation models and tall tower EC measurements generally agree that the region was a carbon sink during the three-year study period (2016)(2017)(2018). However, simulations by the four tested models differed in estimates of the sink strength and responses to the 2018 drought. They also did not perform well in reproducing the inter-and intra-annual variability and seasonal dynamics. Based on consistency of the simulations and observations of flux and CO 2 concentrations, we conclude that the 2018 drought had limited impact on the Northern Sweden region's capacity to sequester carbon, although it had severe effects in southern Sweden, Denmark, and Estonia . The high correlations between modeled and observed daytime concentrations demonstrate the high quality of the transport model and background concentration estimates used in the study, especially during daytime and summer. In addition, our regional sensitivity analysis confirmed that the source area/region of fluxes affect the measured concentrations. Thus, our study indicates that the ICOS atmospheric observation network has substantial potential utility for validating the modeled fluxes in the region in an inversion framework. However, we also noted discrepancies that may at least partly reflect a need for further refinement of these general vegetation models for use in Nordic countries, and their evaluation at other ICOS sites to improve their flux estimates for transport models. Overall, we conclude that the integration of vegetation and transport models with tall tower flux and concentration measurements is a promising approach for reconciling regional carbon balances and their responses to climate change.

Data Availability Statement
The modeled concentrations, cumulated sensitivities, background concentrations, etc as well as the regions definition are shared via the figshare link (https://figshare.com/articles/dataset/Data/16676656). The first author, Anusha Sathyanadh, acknowledges a Postdoctoral Scholarship funded by the Kempe Foundations (grant no: JCK-1815). A grant (no. 942-2015-49) from the Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning (FORMAS) funded the tall tower eddy covariance measurements at Svartberget. Financial support from the Swedish Research Council and research institutes contributing to the Swedish Integrated Carbon Observation System (ICOS-Sweden) and Research Infrastructure from the Swedish Infrastructure for Ecosystem Science (SITES) is also acknowledged. Marko Scholze and Guillaume Monteil acknowledge support from the three Swedish strategic research areas ModElling the Regional and Global earth system (MERGE), the e-science collaboration (eSSENCE), and Biodiversity and Ecosystems in a Changing Climate (BECC). The computations were enabled with resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC partially funded by the Swedish Research Council through grant agreement no. 2018-05973. We thank Greet Janssens-Meanhout (European Commission, Joint Research Center, 790 Ispra, Italy) for providing the fuel type and category specific version of the EDGAR v4.3 anthropogenic emission data and Hugo Denier van der Gon (TNO, The Netherlands) for making the temporal emission profiles available. We also thank the PIs of the atmospheric stations at ICOS Svartberget, Norunda, Hyltemossa, Puijo, Pallas, and Hyytiala (SE-SVB, SE-NOR, SE-HTM, FI-PUI, FI-PAL, and FI-SMR) for providing data on atmospheric CO 2 concentrations. The authors thank the staff at the SLU Unit for Field-Based Forest Research for technical and logistic support at the Svartberget station.