NEE estimates 2006-2019 over Europe from a pre-operational ensemble-inversion system

. 3-hourly Net Ecosystem Exchange (NEE) is estimated at spatial scales of 0.25 degrees over the European continent, 10 based on the pre-operational inverse modelling framework “CarboScope Regional'' (CSR) for the years 2006 to 2019. To assess the uncertainty originating from the choice of a-priori flux models and observational data, ensembles of inversions were produced using three terrestrial ecosystem flux models, two ocean flux models, and three sets of atmospheric stations. We find that the station set ensemble accounts for 61% of the total spread of the annually aggregated fluxes over the full domain when varying all these elements, while the biosphere and ocean ensembles resulted in much smaller contributions to the spread of 15 28% and 11%, respectively. These percentages differ over the specific regions of Europe, based on the availability of atmospheric data. For example, the spread of the biosphere ensemble is prone to be larger in regions that are less constrained by CO 2 measurements. We further investigate the unprecedented increase of temperature and simultaneous reduction of Soil Water Content (SWC) observed in 2018 and 2019. We find that NEE estimates during these two years suggest an impact of drought occurrences represented by the reduction of Net Primary Productivity (NPP), which in turn lead to less CO 2 uptake 20 across Europe in 2018 and 2019, resulting in anomalies up to 0.13 and 0.07 PgC yr -1 above the climatological mean, respectively. Annual temperature anomalies also exceeded the climatological mean by 0.46 °C in 2018 and by 0.56


Introduction
The atmospheric mole fractions of Greenhouse Gases (GHGs) like CO2, CH4, and N2O have drastically increased since the 30 industrial era began (Friedlingstein et al., 2019), primarily caused by anthropogenic GHG emissions. As a consequence, the globally averaged surface air temperature has risen by 0.87 °C from 1850 to 2015(Jia et al., 2019. Carbon dioxide is ranked as the most prominent anthropogenic GHG owing to its atmospheric abundance, resulting from a) the natural exchange through the biogeochemical interactions with the organic molecules in the biosphere and hydrosphere (represented by the Net Primary Productivity -NPP), b) significant anthropogenic emissions from burning of fossil carbon and from cement production, and c) 35 land use changes such as deforestation. The largest uptake of atmospheric CO2 is carried out through terrestrial Gross Primary Production (GPP) and thought to derive an uptake of about one-third of anthropogenic emissions owing to enhancement of photosynthetic CO2 uptake in the recent decades (Cai and Prentice, 2020). However, measurements of NEE cannot be easily achieved at finer spatial and temporal scales over the globe. Ancillary data from the atmosphere and the biosphere are thus applied in the inverse modelling setups to estimate the natural CO2 fluxes. Such a method of using atmospheric data to constrain 40 NEE obtained from the terrestrial biogenic models is also called a top-down method.
The continuous expansion of GHG in-situ measurement capabilities enabled atmospheric tracer inversion systems to better infer the sources and sinks of CO2 at global (Ciais et al., 2010;Enting et al., 1995;Kaminski et al., 1999;Rödenbeck et al., 2003), and regional scales Kountouris et al., 2018b;Lauvaux et al., 2016). Meanwhile, regional atmospheric inversions have employed atmospheric transport models at finer spatial resolution to deal with the complex 45 atmospheric circulation at continental measurement stations.
The observational site network across Europe has been markedly homogenized since the Integrated Carbon Observation System (ICOS) was established in 2015, allowing for better estimation of the regional budgets of CO2 over Europe (Monteil et al., 2020). Consequently, this has allowed for a better understanding of the impacts of climate extremes on the ecosystem productivity such as the drought episode that occurred in 2018 (Bastos et al., 2020;Rödenbeck et al., 2020;Thompson et al., 50 2020). The inversions typically assume anthropogenic emissions to be well known, and thus target the more uncertain biosphere-atmosphere fluxes.
The regional inversion framework encounters various sources of errors, such as 1) the uncertainty of a-priori knowledge (necessary in the Bayesian framework inversions to regularize the solution of the ill-posed inverse problem), and 2) the representation error resulting from the inaccuracies in simulating the atmospheric transport. The structure of prior error (e.g., 55 uncertainties in the prior biosphere flux estimates) is of particular importance as it determines the way in which the flux corrections calculated from the data information should be spread in space and time (Chevallier et al., 2012;Kountouris et al., 2015). Defining proper error covariance matrices in both flux and measurement space is therefore essential to obtain an optimal estimate of the true fluxes. Non-optimized flux components used as prescribed fluxes in the inverse frameworks should be provided with the highest achievable confidence, as any error in these components will directly modify the estimated 60 biosphere-atmosphere fluxes.
Here, we present NEE estimates from a pre-operational regional inversion system set up over Europe covering fourteen years since 2006. An ensemble is created by varying a) a-priori biogenic fluxes, b) a-priori ocean fluxes, and c) the number of available atmospheric observation sites in order to estimate their impact on a-posteriori optimized biogenic fluxes. We furthermore discuss the interannual variability (IAV) over this period, with special focus on the changes of NEE in 2018 and 65 2019, specifically in light of the water availability and temperature variations that occurred in the wake of anomalously warm and dry conditions over the continent. These changes are analysed using the seasonal and annual NEE fluxes aggregated over different subregions in Europe.
The inversion setup, observational dataset, and prior fluxes used are described in Section 2, including details on ensemble members configuration. A statistical analysis of uncertainty and spreads over the ensembles of inversions are presented and 70 discussed in Section 3.1. Section 3.2 presents the NEE estimated in the pre-operational inverse system based on several analysed cases. Finally, discussions and conclusions are summarized in Section 4.

Inversion framework
The CarboScope Regional inversion system (CSR) is used to infer NEE from observed atmospheric CO2 dry mole fractions at 75 high spatiotemporal resolution over Europe. The CSR makes use of Bayesian inference to regularise the solution of the underdetermined inverse problem (i.e., there are more unknown fluxes than atmospheric observations). For details about the mathematical concepts, we refer the reader to Rödenbeck (2005); the specifications of the set-up largely follow previous studies by Kountouris et al. (2018aKountouris et al. ( , 2018b. The inversion at the regional domain is embedded into the global atmosphere using the "two-step scheme" described in Rödenbeck et al. (2009). 80 Atmospheric transport is simulated by the Stochastic Time-Inverted Lagrangian Transport (STILT) model , which is utilized to calculate hourly surface influences at the stations (i.e., "footprints") at the spatial resolution of 0.25° (longitude) × 0.25° (latitude). The model is driven by meteorological fields from the high-resolution implementation of Integrated Forecasting System (IFS HRES) model of the European Centre for Medium Range Weather Forecasts (ECMWF).
The upstream influence is simulated over the past 10 days by releasing 100 virtual particles at the sampling heights of stations 85 (receptors). Additionally, we use the Eulerian global model TM3 (Heimann and Körner, 2003) within the CarboScope global inversion framework (Rödenbeck et al., 2003) to provide the far field contributions to the regional domain at a coarser spatial resolution of 5° (longitude) × 4° (latitude).

Atmospheric data
Since CO2-dry-mole-fractions are the main constraint of the inversion system, we have aimed in our study to maximize the 90 data coverage by using the observations available through the ICOS network as well as further atmospheric observation sites (both ICOS-associated and independent). All of the datasets are high-quality products of the level 2 ICOS Atmospheric Release, which underwent a strict filtration procedure described in Hazan et al. (2016). This homogenized data treatment makes the data suitable for inverse modelling. For measurement sites with multiple sampling levels, the top one is chosen, as this one is expected to be represented best in the STILT transport model. 95 The core of our observational dataset consists of data from 44 sites collected in the ICOS Carbon Portal under the 2018 drought initiative (https://doi.org/10.18160/ERE9-9D85), covering the period 2006-2018. This base dataset was extended into 2019 by the level 2 data (L2) released by ICOS in 2020, as well as included data from four new sites. From the total number of sites, 23 are currently ICOS-labelled and provide data since 2015, while the rest are non-labelled sites, providing datasets since 2006. Figure 1 shows the distribution of all sites throughout the domain of Europe. The figure also shows the division of the 100 domain into six sub-regions (North, South, West, East, South-east, and Central Europe) used for post-processing, to outline the impact of the observational constraint distribution on posterior fluxes.

South-eastern Europe (light red).
The representation error is assumed to be specific for different station types, which are categorized to 5 classes. Weekly representation errors are presented in Table 1, defining the measurement error covariance matrix in the cost function. The 105 observations are mostly provided at hourly frequency, especially in recent years. We also include measurements from flask sampling (mostly weekly) when available from the corresponding sites. To better represent the well-mixed boundary layer in the STILT model, we limit our analyses to measurements of 6-hour day-time for all stations, i.e., 11:00 to 16:00 (local time), except for mountain sites. For the latter, night-time hours (23:00 to 04:00, local time), are chosen, during which these sites regularly sample the residual layer. Particularly before establishing ICOS in 2015, the variability of station data coverage 110 across the period of interest was rather high, underlining the sparsity of available data (Fig. 2). Since then, the site network over Europe has been markedly expanded as new stations were installed. It should be noted that the variability in data coverage is expected to result in an inconsistency of annual flux variations. Therefore, we combined stations into 3 subsets: a) the full set of stations, referred to as "all sites", b) a subset of 15 stations that have consistent coverage from 2006 to 2019, to which we will refer as "core sites"), and c) a third subset of 16 stations that do not have gaps longer than 1 month during [2015][2016][2017][2018][2019]115 subsequently referred to as "recent sites". The far field contributions provided to the regional domain were calculated in the two-step inversion approach (Rödenbeck et al., 2009) using a global observational record from 75 sites (doi:10.17871/CarboScope-s10oc_v2020), which has best data coverage in the 2010-2019 period.

A-priori-fluxes
Terrestrial ecosystem flux models are utilized to provide prior knowledge of biogenic fluxes (NEE, defined as the net 120 ecosystem exchange). To appropriately represent the diurnal cycle in our modelling framework, NEE is obtained from the biosphere models at hourly temporal resolution. Three biosphere models were used as priors in the inversion runs. The first is the Vegetation Photosynthesis and Respiration Model, VPRM (Mahadevan et al., 2008). VPRM is a diagnostic model driven by shortwave radiation and temperature at 2-meter from the ECMWFs high resolution operational forecast product (IFS HRES). To calculate NEE and respiration fluxes, it uses MODIS (Moderate Resolution Imaging Spectroradiometer) indices 125 derived from surface reflectance, namely Enhanced Vegetation Index (EVI) and land surface water (LSWI) together with typespecific vegetation parameters optimized against the eddy covariance (EC) data. Parameter values for VPRM previously used by Kountouris et al. (2018b), were updated using the most recent EC data, and are available at https://www.bgcjena.mpg.de/bsi/index.php/Services/VPRMparam. The second biosphere prior is from the data driven modelling approach FLUXCOM, which combines eddy covariance measurements and satellite observations in several machine learning algorithms 130 to quantify the surface-atmosphere energy and carbon fluxes (Jung et al., 2020). Here, we use an extension of the modelling set-up described in Bodesheim et al. (2018), which employs daily and hourly surface meteorological information from ERA5 as well as a mean annual cycle of satellite data to produce NEE estimates at hourly temporal and 0.5-degree spatial resolution.
It should be noted that the magnitude of interannual changes in the data-driven flux estimates is generally found to be unrealistically small (Jung et al., 2012). The terrestrial biosphere model SiBCASA (Schaefer et al., 2008) is used as the third 135 biosphere model. SiBCASA is a combined framework based on the Simple Biosphere (SiB) model and the Carnegie-Ames Stanford Approach (CASA) model. As explained in Schaefer et al. (2008), Gross Primary Productivity (GPP) is calculated by SiB assuming that it is in balance with heterotrophic and autotrophic respiration (RH, RA), meaning that diurnal and seasonal variations are well represented; however, the long-term terrestrial carbon changes cannot be predicted by SiB component alone.
CASA fills this gap as it includes a light-use-efficiency model to estimate Net Primary Productivity (NPP, equal to GPP -RA). 140 In turn, CASA cannot calculate NEE during night time. Therefore, combining both models in a hybrid version of SiBCASA combines the advantages in their biophysical and biogeochemical aspects to calculate NEE from RH, RA (SiB) and NPP (CASA).
Following Kountouris et al. (2018b), we assume that the spatial correlation of the prior uncertainty follows a hyperbolic decay function, similar as in the inversion case nBVH described in that study. One notable difference in the current work is the 145 improved implementation of the directional dependence, with a twofold increase of decay distance in meridional direction.
Temporarily, prior errors are assumed to be correlated over 30 days, as found in Kountouris et al. (2015).
Ocean CO2 fluxes and anthropogenic emissions are considered as prescribed fluxes in the inversion system. Ocean fluxes are taken from two sources at coarse spatial resolution, (5 × 4 degrees): a climatological flux product with monthly fluxes of Mikaloff Fletcher et al. (2007), and the CarboScope pCO2-based ocean flux, providing fluxes at 6-hourly temporal resolution 150 (Rödenbeck et al., 2013). Fossil fuel emissions are taken from the category and fuel-type specific emission inventories of EDGAR-v4.3, and processed following the COFFEE approach (Steinbach et al., 2011) to include diurnal, weekly , seasonal, and annual variations. These emissions are updated annually according to national consumption data from the BP (British Petroleum) statistical review of world energy (BP, 2019), and are available from the ICOS Carbon Portal under https://doi.org/10.18160/Y9QV-S113. 155

Set-up of the inversion runs
We conduct three ensembles of inversion runs listed in Table 2 utilizing different setups of prior products (biosphere and ocean ensembles), as well as selected sets of observational data (station set ensemble). The inversion runs are labelled with unique codes for reference. B0 is defined as the base case of our analysis. It is configured using default settings of the inversion runs, with biogenic fluxes from VPRM, climatological ocean fluxes, and using all available atmospheric data as input. In the 160 biosphere ensemble (consisting of B0, B1, B2), FLUXCOM and SiBCASA replace the VPRM model in both B1 and B2, respectively, allowing for distinguishing the effect of using different prior flux models on posterior NEE. In the ocean ensemble, we replace the climatological ocean fluxes used in B0 with the pCO2-based CarboScope ocean fluxes in O1. The station set ensemble is formed by running the inversion with varying measurement station subsets: B0 -all sites, S1 -core sites, and S2 -recent sites, as explained in Section 2.2. For each of the three ensembles of inversions, its spread is calculated 165 as the standard deviation of the differences between each ensemble member and the ensemble mean over the respective overlapping period of time. The statistical uncertainty calculations are calculated in the inverse system based on model-data mismatch, performed for the base case inversion (B0) as they remain identical independent of the biosphere model used.

Statistical analysis of ensemble uncertainties 170
The annual NEE estimates among the biosphere ensemble ( Fig. 3) show good agreement across the three biosphere models but also across S1 and O1 inversions, yielding similar budgets of CO2 fluxes over the full domain. The findings suggest that atmospheric data constraints are more dominated in the posterior NEE fluxes in comparison with the prior constraint. are likely to entail this signal from the meteorological data used to force these models. However, the VPRM model overestimates the mean CO2 uptake compared to the a-posteriori fluxes, while SiBCASA underestimates the mean CO2 uptake, and this dissimilarity is persistent for all years as well.
The statistical uncertainty and spreads over the ensembles are evaluated, and affect our data (Fig.3). It is noticed that the spreads over the posterior fluxes and prior fluxes are comparable with the corresponding uncertainties over the full domain 185 (All Europe), Central Europe, and Northern Europe. There is a clear reduction of uncertainty and spread in posterior fluxes either over the full domain (All Europe) or in subregions like Central and Northern Europe. Unlike prior uncertainty, posterior uncertainty slightly differs from year to year following the number of atmospheric sites available (Fig 2). This gets even clearer when looking at the marked reduction of the posterior uncertainty in Central Europe, as well as in regions with high station density, resulting in a stronger observational constraint. In contrast to Central Europe, a smaller reduction of posterior spread 190 is found in Northern Europe as well as in other regions where there are few or no stations (e.g., Eastern Europe, not shown).
In this case, NEE estimates are not well constrained by atmospheric data. Instead, a-posteriori flux is driven by the inversion using biosphere models and their uncertainty, particularly for the distant areas that cannot be constrained by observations through the spatial correlation. Table 3 denotes the reduction of the biosphere ensemble spread in the a-posteriori relative to the a-priori over the full domain, Central and Northern Europe. It indicates less reduction in Northern Europe due to the 195 sparseness of observational sites. The large reduction of spread in Central Europe reflects a notable dependency of NEE estimates on the atmospheric measurements, substantially where the observation network is dense.
To analyse the seasonal variations, seasonal cycle from B0, B1, B2, S1, and O1 inversions is averaged over 13 years for the full domain of Europe together with the corresponding biogenic prior fluxes for VPRM, FLUXCOM, and SiBCASA (Fig.4).
Results show good agreement amid a-posteriori results of all inversions, while prior biosphere models show large differences, 200 a pattern similar to the one seen over the annual fluxes in Fig.3. Nevertheless, posterior NEE fluxes estimated in S1 inversion show differences during May-August when compared with the estimates of other runs, reflecting a larger sensitivity of IAV to summer fluxes when applying a different set of stations. In addition, the difference of posterior fluxes seen in Fig. 3 over the annually aggregated estimates computed from B1 inversion over the period 2014-2018 largely results from the estimates during May and June when comparing it to the rest of biosphere ensemble elements (Fig. 4). 205 Figure 5 illustrates the statistical uncertainty and the spread through the overall ensembles of inversions (listed in Table 2) calculated annually over three regions. As was discussed in the time series of NEE (Fig. 3) indicating large discrepancies between prior flux models. This confirms that the prior uncertainty assumed in the CSR system is realistic. The IAV of B0 was calculated separately for prior and posterior fluxes (blue bars) from the anomalies relative to the long-term mean to reveal the magnitudes of interannual deviation in comparison with the spread variability.
In terms of the spread of the biosphere ensemble, the standard deviation of posterior fluxes declines from 0.666 to 0.032 PgC 215 yr -1 over All Europe. This reduction is 95.1%, 96.0%, and 74.8% in All Europe, Central Europe, and Northern Europe, respectively. Spatial differences are expected as stations are not evenly distributed across the domain of Europe. This can be noticed from the spread over Central Europe (a large number of stations, 18 sites) and Northern Europe. (a lesser number of stations, 8 sites). As a result, lack of observations leads to inflating the spread over the biosphere ensembles.
The largest impact on NEE estimates in the ensembles are observed when the spread over station set ensemble is analysed. In 220 this regard, a robust analysis can be based on a subset from Central Europe, as the subsets of stations in this region are clearly contrasting in the two ensemble members (core sites and recent sites). The spread of the station set ensemble was found to be 0.11 PgC yr -1 -larger than those resulting from the biosphere and ocean ensembles (0.05 and 0.02 PgC yr -1 , respectively).
Noteworthy, the spread in the station set ensemble is slightly larger than the statistical uncertainty, highlighting the importance The spatially distributed spread of all ensembles is depicted in Fig. 6. In this instance, the standard deviation was calculated 230 for each grid cell rather than aggregating fluxes over regions first and then computing the spread (Fig. 5). The spatial spread here illustrates the deviations of the biosphere ensemble (Biosphere spread), the biosphere models (Prior biosphere spread), the station set ensemble (Stationset spread), and the ocean ensemble (Ocean spread). The maximum spread of 0.191×10 -4 (PgC yr -1 ) was observed over the a-priori terrestrial biosphere models, particularly concentrated in Central and Southern Europe.
The spread of a-posteriori biosphere ensemble is significantly reduced. In the station set ensemble, isolated stations like 235 Hegyhatsal in Hungary and Sierra de Gredos in Spain demonstrate relatively high impact on the NEE spatial patterns over broader areas, reflecting the inversion correlation length. However, such an impact is not clearly realized on the "Stationset spread" map amid dense clusters of sites due to the commutative constraint that compensates for the excluded sites within the subsets of stations. These results highlight the importance of defining a proper function of spatial correlation decay in the prior error structure. A quite small influence is seen through the spread over ocean where a slight impact emerges only in wider 240 coastal regions, being almost negligible inland (e.g., in Central and Eastern Europe).

NEE estimates of 2018 and 2019 in a preoperational system
In this section, we present CO2 fluxes for two selected years estimated in a preoperational system in the context of long-term estimates. The period of interest is chosen to start from 2006 in which a better coverage of observations exists within the 250 domain of Europe. Here, we give special attention to analysis of the drivers of spatiotemporal differences in line with climate disturbances that occurred in 2018 and 2019, during which inaccuracies of estimating the continental fluxes of CO2 have been reported (Friedlingstein et al., 2019). This is due to the sensitivity of ecosystem respiration and photosynthetic fluxes to extreme events like lasting droughts. The analysis here is based on two inversion runs using observational data only from the subset of core sites that have consistent measurements, i.e., as used in scenario S1 discussed in previous sections. The choice of using 255 consistent measurements is essential to study the IAV to diminish the uncertainty caused by gaps of data coverage over years. in the biosphere models from the meteorological data, as well as less information of radiation reflectance obtained from the 275 remote sensing data due to dominant cloudy scenes in winter, provided that the VPRM and FLUXCOM models use forcing data from meteorology and remote sensing.
To assess the temporal changes of NEE in response to such climate variations, we compare the seasonal anomalies of NEE (prior and posterior) to the anomalies of 2-m air temperature and Standardized Precipitation and Evapotranspiration Indices (SPEI) (Beguería et al., 2014) during spring, summer, fall, and winter, as well as the annual mean (Fig. 9). Here we show 280 estimated NEE integrated over the full domain. Monthly near-real-time data of SPEI (SPEI01) are obtained from https://spei.csic.es/map/maps.html at 1° spatial resolution and monthly 2-m air temperature accessed from https://psl.noaa.gov/data/gridded/index.html at 0.5 spatial resolution. The anomalies were normalized with the standard deviation of the interannual variations since 2006. In addition, Pearson correlation coefficients between posterior fluxes, prior fluxes, temperature and SPEI for the full year and calendar seasons were calculated. Of note, due to the fact that the biosphere 285 model VPRM utilizes temperature from meteorological fields and EVI data from the satellite sensor MODIS, it is anticipated to systematically correlate with temperature and SPEI. Consequently, we mostly devote our comparison to the posterior fluxes.
The findings of standardized anomalies in Fig. 9a show that the decrease of CO2 uptake in 2018 and 2019 summers (positive NEE anomalies) was exceptional and concurrent with a profound deficit of SWC (negative SPEI anomalies, dry conditions).
The reduction/very low SWC also coincided with an unprecedented rise of temperature (positive T anomalies, highest in 2018) 290 across Europe. Being an indicative factor of drought occurrences, SPEI links water availability in the surface including soil moisture (crucial limitation of GPP, especially in the temperate regions) and temperature to the precipitation and evapotranspiration rate. Hence, there is quite good agreement between posterior NEE and SPEI not only at spatial scales but also at temporal scales. The standard deviations of the interannual variations of posterior NEE, SPEI, and temperature over all Europe in the annual mean through the 14 years were equal to 0.17 PgC yr -1 , 0.12, and 0.45 K, respectively. When relating the 295 changes occurring during 2018 and 2019 to the context of the previous 12 years, the annual anomalies of SPEI found to decline to more than twice as much as the climatological deviation around -0.29 in 2018 and to around -0.08 in 2019. As a consequence, posterior NEE anomalies increased to 0.14 and 0.08 PgC yr -1 above the climatological mean in 2018 and 2019, respectively.
The excess of annually averaged temperature was predominant in 2018 and 2019, reaching around 0.40 and 0.47 °C above the climatological mean, respectively. Despite the fact that the impact of the 2018 and 2019 drought on NEE is realized from the 300 SPEI and temperature anomalies, there is a relatively moderate correlation between estimated NEE and SPEI and temperature at the annual scale over the full domain (Fig. 9b). However, the correlation coefficients largely vary between seasons. A high anticorrelation (-0.86) between estimated NEE and temperature is found to be consistent during spring. In contrast, anticorrelation drastically decreases and turns out to almost vanish in summer and fall. Nonetheless, the relatively moderate correlation of SPEI with posterior NEE during summer is adequate to deduce the lack of SWC under dry conditions through 305 the anticorrelation between SPEI and temperature. This implies, warm conditions accelerate the depletion of soil moisture content, in particular on the soil top layer that lacks water content in its deeper layers to compensate for the higher evaporation rate at the surface layer. This affects photosynthesis efficiency during the growing season by decreasing the gross primary productivity, but also increases the contribution of soil respiration, more pronounced in 2019.
During winter, water availability does not seem to be a limiting factor of NEE as we notice from low correlations between 310 posterior NEE and SPEI. Instead, temperature negatively correlates with posterior NEE indicating that the increase of temperature coincides with enhanced uptake of CO2 where photosynthesis can occur, e.g., in evergreen areas. But cold years have more and longer snow cover which decreases photosynthesis, whereas soil heterotrophic respiration contributes more to CO2 release since soil temperature is expected to be larger than air temperature owing to snow cover insulation. Figure 10 illustrates the seasonal contribution of NEE to IAV, which is dominated by summer and spring variability in comparison with 315 winter. We note that the posterior fluxes are in agreement with the biosphere model VPRM in summer and fall, while the variability of posterior fluxes is larger during winter and smaller during spring compared to prior fluxes; the opposite holds true for the prior fluxes. Temperature is, however, shown to largely vary during winter, while SPEI contribution does not show a significant variability between seasons.

Spatial differences of NEE estimates in 2018 and 2019 320
Using identical observations for the last 5 years, S2 inversion demonstrates the differences between NEE estimates in 2018 and 2019 as noticed from Fig. 11. Results emphasise the aftermath of drought episodes, showing a smaller uptake of CO2 in The smaller spread in the a-posteriori fluxes found through the ensembles of inversions is evident over All Europe reflecting the good performance of the inversion system. In the biosphere ensemble, flux estimates are not very sensitive to a-priori 355 terrestrial ecosystem fluxes. We deduce this from the small spread over the a-posteriori fluxes (Fig. 3), occurring despite major differences in a-priori fluxes. Likewise, different ocean flux models have the smallest effect on estimating NEE, in particular inland, where ocean-land exchange is dissipated. However, the spread in the station set ensemble is strongest, at 0.11 PgC yr -1 for the annually aggregated fluxes over the full domain. This points out a higher sensitivity of the inversion to the number of stations in comparison with 0.06 and 0.02 PgC yr -1 spreads in the biosphere and ocean ensembles, respectively. This effect is 360 most pronounced over Central Europe, where measurements of CO2-dry-mole-fractions are available from a large number of stations, and thus a contrasting number of sites manifests amongst station subsets, given that the station subsets were not selected based on geographical locations but on the long record of data coverage. Further, this finding corroborates the dominant influence of observational constraint on NEE estimates in the biosphere and ocean ensembles seen through interannual variations. Such influence is, however, subject to the availability of the atmospheric data which can otherwise be 365 altered by the prior constraint to maintain the posterior estimates according to Bayes' approach. For example, the relative spread over biosphere ensembles in Northern Europe, having less dense coverage, increases to 38.4%, compared to 23.9% in Central Europe. Conversely, in the stationset ensemble, the spread calculated for Northern Europe was found to be smaller than that in Central Europe (51.8% versus 71.7%), reflecting the biogenic flux constraints in this case, given the smaller differences of observations among the station subsets in this ensemble. The ocean ensemble spread remains at a minimum 370 percentage in both regions, albeit it elevates to 9.6% in Northern Europe in comparison with only 4.3% in Central Europe.
Although the impact of ocean fluxes on NEE estimates can be negligible over the full domain and inlands, the results point to a relatively higher influence in the coastal regions.
Our results denote a comparable reduction of posterior uncertainty and the spread relative to their a-priori (Fig. 5). Noteworthy, the indirect effect of statistical posterior uncertainty on the corresponding spread over the ensembles of inversions emerges 375 from the common dependency on observational data, which predominantly appear in the well-constrained areas in Germany, France, Benelux, and the U.K. Over such regions, the posterior uncertainty and spreads are greatly reduced and the inversions tend to converge regardless of which prior flux model is used. It is essential to consider the prior error assumption, as well as the prior error structure in the spatial and temporal aspects. This will determine to which extent the posterior fluxes are dependent on the prior biogenic fluxes, specifically in the regional inversions where the degrees of freedom can drastically 380 increase following the finer spatial and temporal resolution of biosphere flux models and atmospheric transport models.

Response of NEE estimates to climate variation
The linkage between NEE and climate variation has been examined via SPEI and temperature as proxy data of climate  (Ma et al., 2020). But GPP during 2019 spring showed a higher efficiency (larger uptake of CO2) than the spring of 2018, benefiting from the increment of temperature, SWC, and light availability. The finding is consistent with a study on seasonal NEE over North America implemented by Hu et al. (2019) and seems to hold for Northern regions where temperature is substantially considered as a limiting factor to NEE. Additionally, our results showed that Central 390 Europe experienced higher sources of CO2 during 2019, which can be impacted by an extended legacy of the drought of 2018, where forests were profoundly stressed and thus their growth was negatively impacted. MacKay et al. (2012) found about 17% reduction of drought plot growth relative to a reference plot and showed that growth reduction in the forests across Europe exceeded this value under drought conditions depending on tree species. Furthermore, the ecosystem respiration response to the temperature increment may contribute to such a positive anomaly, given that temperature anomalies during 2019 over 395 Central and Southern Europe were unprecedented in line with the 2018 anomaly (Hari et al., 2020). In agreement with Rödenbeck et al. (2020), we found that summer NEE anomalies were in agreement with the anomalies of temperature and SPEI, occurring in different summers including 2018.
In terms of estimated winter fluxes, the medium anticorrelation between temperature and NEE (also shown by the anomalies in Fig. 9) implies that an increase of CO2 release occurs in cold winters when thicker and lasting snow cover prevail. Snow 400 insulation prevents the soil temperatures from decreasing as strongly as air temperature. In comparatively warmer soils respiratory carbon dioxide emissions by the soil biota are sustained also during winter. Being in agreement with our results, Monson et al. (2006) found that soil temperature during winter increases nonlinearly with the depth of snowpack resulting in enhancement of soil respiration. Even though the study was carried out at site level in the Northern Hemisphere in the U.S, the agreement with our findings refers to a similar impact at a widespread scale. 405 Besides the previous explanation of T-NEE anticorrelation in winter, one should take into account a hypothesis of the vertical mixing height being systematically biased in the atmospheric transport models. In this case, a misrepresentation of mixing height during colder winters might occur, as was demonstrated in a study conducted by  devoted to characterizing the uncertainty of atmospheric transport models resulting from vertical mixing changes during day and night time. Given that, the shallow boundary layer developed during night holds similar characteristics during cold winters. A 410 separate study is required to investigate such an impact which will lead to improving vertical mixing in the atmospheric transport models, and thus reducing the uncertainty in atmospheric tracer inversions. Apart from that, the contribution of winter fluxes to the interannual variations is lower than that resulting from summer fluxes, denoting a lesser impact.
Overall, the NEE response to SWC and temperature varies depending on the temporal and spatial aspects of the region of interest, and is connected to the hydrological cycle and physical dynamics of soils and the ambient atmosphere. 415 In this paper, the NEE flux budgets of 2018 and 2019 are estimated in a pre-operational method to keep track of the changes of net terrestrial fluxes of CO2. Results still suggest the domain of Europe as a net sink of CO2, albeit in 2018 and 2019 CO2 uptake decreases to -0.22 ± 0.05 and -0.28 ± 0.06 PgC yr -1 , respectively, due to drought occurrences compared to a multi-year average uptake of -0.36±0.07 PgC yr -1 (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019). In contrast to the a-posteriori, the prior fluxes obtained from VPRM are by far overestimating CO2 uptake, especially in the growing season, resulting in -1.47 ± 0.43 and -1.37 ± 0.43 PgC yr -1 in 2018 420 and 2019, respectively. The posterior fluxes constrained in the regions or individual countries whose observational data are denser, such as France, Germany, and the U.K, are more reliable than those that have sparse observation network (e.g., Poland or Spain) or do not have any monitoring sites at all (e.g., Turkey). This implies that NEE in some countries within the domain of Europe is constrained either by weak observational signal from the neighbouring regions through the spatial correlation length of prior error, or mostly dominated by the prior flux constraint. This kind of flux constraints skews the total aggregated 425 fluxes at the expense of well-constrained regions and is expected to amplify the posterior uncertainty. For example, when comparing 2019 NEE budgets in France and Turkey, the a-posteriori fluxes are summed up to -0.032 ± 0.008 and 0.045 ± 0.020 PgC yr -1 , respectively. Despite the fact that they have different areas, the uncertainty in the latter is about 45%, notably greater than that realized in the first, up to 25%. These results emphasise the importance of using a wider coverage of CO2 observations in the regional inversions to better estimate the continental flux budgets, and thus understanding the biogenic flux 430 changes amid climate variation.

Code and Data availability
NEE estimates and the respective prior fluxes, as well as codes used in this study can be made available upon request.
Atmospheric dry mole fraction measurements of CO2 are available from the ICOS Carbon Portal and can be accessed from 435 http://doi.org/10.18160/GZ5S-4GPR. SPEI01 can be downloaded from the near-real-time data through https://spei.csic.es/map/maps.html.

Competing interests
Some authors are members of the editorial board of journal Atmospheric Chemistry and Physics. The peer-review process was 440 guided by an independent editor, and the authors have also no other competing interests to declare   BIR  BIS  LUT  PAL  CRP  LMP  MHD  WAO  BI5  BR5  CB4  GA5  HP4  HT3  HU4  KIT  KR3  LIN  LMU  NO3  OP3  PUI  RG2  SM3  ST5  SV3  TA3  TO3  TR4  HEI  GIC  LHW  ER2  FKL  MLH  SA3  UTO  BS3  IP3  OHP  CMN  JFJ  KAS  OX3  PD2  PTR  PUY  SI2   C  C+R  C+R