Climate impact on mean annual cycle and interannual variability of CO₂ fluxes in European deciduous broadleaf and evergreen needleleaf forests: insights from observations and state-of-the-art data-driven and process-based models

The impact of climate on the annual cycle and interannual variability of CO₂ fluxes is assessed in European evergreen needleleaf forests (ENFs) and deciduous broadleaf forests (DBFs) using observations from 19 sites, alongside outputs from process-based and data-driven models. All models capture the temporal phasing of CO₂ fluxes, including a shorter sequestration period in northern rather than southern Europe, a more pronounced annual cycle for DBFs than for ENFs in central Europe and strong interannual variability across sites. However, they generally underestimate both the magnitude of CO₂ sequestration and its interannual variability compared to observations. Regarding the annual cycle, all datasets indicate enhanced CO₂ uptake from late spring to early fall, with a stronger climate–CO₂ flux coupling in northern and central Europe than in southern Europe, where seasonality is less pronounced. At the interannual timescale, the climate does not show a significant influence on observed and modelled net ecosystem exchange (NEE) when correlations are computed using monthly anomalies across all months combined. This apparent lack of relationship conceals meaningful seasonal patterns. In winter and fall, NEE tends to be positively correlated with temperature, soil moisture and vapour pressure deficit (VPD). In spring, NEE shows negative correlations with temperature and VPD but a positive correlation with soil moisture. The summer pattern is reversed compared to the spring pattern. In the observations, these relationships are noisy in both time and space, suggesting strong site-specific effects. In contrast, the models exhibit more structured and spatially coherent patterns with strong correlations, which may reflect an exaggerated response to climate forcing despite underestimated magnitude in CO₂ flux interannual variability.

1 Introduction

Forest ecosystems are the largest part of the land CO₂ sink (Lindeskog et al., 2021), with up to 20 %–50 % of anthropogenic CO₂ emissions (land-use changes excluded) sequestered for the 2000–2010 period (Le Quéré et al., 2018; Pugh et al., 2019; Pan et al., 2024). In Europe, recent estimations suggest a slight increase in CO₂ sequestration by forest ecosystems over the 2000–2021 period mainly due to the fertilization effect of increased atmospheric CO₂ concentration (Prentice et al., 2001; Piao et al., 2009; Schimel et al., 2015; Walker et al., 2021; Sitch et al., 2024). However, this trend remains weak, as the fertilization effect has almost been compensated for by a decrease in CO₂ sequestration induced by climate change (Sitch et al., 2024). The extent to which climate controls CO₂ flux exchanges between the atmosphere and European forest ecosystems is thus a burning question in the context of climate change.

Numerous studies have demonstrated the strong influence of climate on CO₂ exchanges between the atmosphere and forest ecosystems. The annual cycle and, to a lesser extent, interannual variability of these fluxes are driven by factors such as incident shortwave radiation, temperature, atmospheric evaporative demand and the water cycle, including soil moisture dynamics (Haszpra et al., 2005; Tang et al., 2014; von Buttlar et al., 2018; Kong et al., 2022; Sharma et al., 2022; Li et al., 2023; Xu et al., 2023). The dominant climate factor influencing CO₂ fluxes depends on the specific component considered. The variability in net ecosystem exchanges (NEEs) is a mixed response of its two components: gross primary production (GPP), which sequesters CO₂ into the ecosystem through photosynthesis, and ecosystem respiration (RECO), which releases CO₂ into the atmosphere from forest metabolism (autotroph respiration) and the decomposition of organic matter by fungi and bacteria (heterotrophic respiration). GPP is primarily driven by vapour pressure deficit (VPD), shortwave radiation, temperature and soil moisture, while RECO is mainly influenced by precipitation, soil moisture and temperature (Messori et al., 2019).

The influence of climate on CO₂ fluxes also depends on several additional factors, with seasonality playing a crucial role. Severe heat waves and droughts acted to reduce CO₂ sequestration in summer at the Europe-wide scale in 2003 (Ciais et al., 2005), in northern Europe in 2018 (Smith et al., 2020; Thompson et al., 2020), and in central and southeastern Europe in 2022 (van der Woude et al., 2023). On the other hand, anomalously high temperature under normal soil moisture conditions in spring set favourable growth conditions and hence increased CO₂ sequestration, such as in northern Europe in 2018 (Smith et al., 2020). The climatic zone under consideration is also a key factor. For instance, GPP is mostly influenced by soil moisture in the Mediterranean region, VPD over parts of central Europe, and temperature over Scandinavia, parts of eastern and southeastern Europe, and higher elevations (Seddon et al., 2016; Madani et al., 2017). The influence of climate on CO₂ flux exchanges is further shaped by various factors, including soil properties (Kurbatova et al., 2008; Besnard et al., 2018; Curtis and Gough, 2018; Martinez del Castillo et al., 2022), forest management practices (Carrara et al., 2003; Saunders et al., 2012), tree age (Kurbatova et al., 2008; Besnard et al., 2018) and tree species (Carrara et al., 2003, 2004; Welp et al., 2007; Kong et al., 2022).

Assessing the impact of climate on CO₂ flux exchanges remains challenging. The main reason involves the scarcity of multi-year CO₂ fluxes measured by eddy covariance above the canopy (Burba, 2021). At the European scale, the Integrated Carbon Observation System (ICOS) network provides standardized and open data from 98 ecosystem stations across 16 countries. The flux tower measurements remain limited in number and temporal depth and unevenly distributed spatially, making it difficult to assess the impact of climate on the interannual variability (and trends) in CO₂ flux exchanges and to map them. Process-based and data-driven models allow us to tackle the above limitations. Process-based models, such as dynamical vegetation models, are routinely used to assess CO₂ flux exchanges between the atmosphere and the biosphere (Friedlingstein et al., 2023). These are mechanistic models (Friedlingstein et al., 2006; Sitch et al., 2008) allowing for testing the response of CO₂ fluxes to individual and combined forcing (Sitch et al., 2024). Data-driven models rely on the identification of statistical relationships between flux tower measures by eddy covariance and corresponding land use, vegetation properties and climate characteristics. Based on these statistical relationships, empirical models are built and used for upscaling, i.e. for assessing CO₂ fluxes in regions where they are not measured (Tramontana et al., 2016; Jung et al., 2019, 2020; Zhuravlev et al., 2022). Both approaches have limitations. Estimations of CO₂ flux exchanges are highly sensitive to physical parameterizations (Cai and Prentice, 2020) and atmospheric forcing (Wu et al., 2017; Hardouin et al., 2022) in process-based models. The reliability of data-driven models is limited by the sparse and uneven distribution of flux tower measurements and by the underlying statistical methods used to build them (Jung et al., 2020). While not perfect, process-based and data-driven models provide satisfactory results for capturing large-scale patterns compared to, for example, satellite estimations (Wang et al., 2023). This makes them valuable complementary tools to observational data.

Most recent studies examining the influence of climate on the temporal dynamics of European forest CO₂ fluxes rely on case studies and primarily focus on spring and summer conditions (Smith et al., 2020; Thompson et al., 2020; van der Woude et al., 2023). However, a more comprehensive assessment is needed across the entire annual cycle, as CO₂ release during fall and winter is expected to increase under climate change. Additionally, climate conditions vary significantly between northern and southern Europe, necessitating a broader spatial perspective. These objectives are addressed at the monthly timescale, which is considered sufficiently fine to capture both the CO₂ flux annual cycle and its interannual variability.

This study addresses these gaps by investigating the impact of climate on both the annual cycle and interannual variability of CO₂ fluxes in European evergreen needleleaf forests (ENFs) and deciduous broadleaf forests (DBFs). Using ICOS network observations alongside state-of-the-art data-driven and process-based model estimates, we first characterize the observed annual cycle and interannual variability of CO₂ fluxes across Europe. We then evaluate model performance in capturing the temporal phasing and magnitude of these fluxes at the site scale. Finally, we assess the influence of climate on both seasonal and interannual CO₂ flux variations, leveraging the extended temporal coverage provided by models.

2 Materials and methods

2.1 Site description

Out of the 24 ENF and DBF sites from the ICOS network, we selected 19 sites (Fig. 1), 13 classified as ENF and 6 as DBF, for which observed CO₂ fluxes are available for at least 5 years (Table 1). These sites allow us to sample the different climatic zones of Europe. Three ENF sites (FR-Bil, FR-FBn and IT-SR2) are located in the northern region of southern Europe, close to 45° N, and ranging from sea level to 400 m in elevation. They are characterized by mild, wet winters and hot, dry summers, with annual mean temperature and precipitation of 12.9–13.90 °C and 700–960 mm, respectively. Four ENF sites (FI-Hyy, FI-Let, SE-Nor and SE-Svb) are located in northern Europe (60–65° N) at an elevation below 270 m. They are characterized by subarctic climate, with annual mean temperature and precipitation of 1.8–6.5 °C and 586–711 mm, respectively. The remaining 12 sites (six DBFs and six ENFs) are situated in central Europe within the 45–60° N, 2.5–20° E domain, encompassing a wide range of elevations (40–1730 m) and spanning temperate to continental climates. As a result, they exhibit substantial variability in annual mean temperature (4.3–11.4 °C) and precipitation (563–1338 mm).

Figure 1 Locations of the CO₂ flux measurement sites from the FLUXNET network across Europe selected for this study. The circles correspond to the six sites located in deciduous broadleaf forests (DBFs), while the triangles indicate the 13 sites located in the evergreen needleleaf forests (ENFs). The vertical colour scale represents the terrain height (above sea level, ASL) in metres. The horizontal colour scale indicates the length of the CO₂ flux record for each site. The elevation layer comes from the National Geophysical Data Center/NESDIS/NOAA/U.S. Department of Commerce (1995): TerrainBase, Global 5 Arc-minute Ocean Depth and Land Elevation from the U.S. National Geophysical Data Center (NGDC), https://doi.org/10.5065/E08M-448. Distributed under CC by 4.0. Publisher's remark: please note that the above figure contains disputed territories.

Table 1 European FLUXNET sites examined in this study.

Columns 1–5 provide information on the site, including its name, ID, latitude, longitude and elevation. Columns 6–8 describe the land cover classification, the period of available records and the record length for each FLUXNET site. Columns 9–10 present the site's mean annual climatic characteristics, including temperature (MAT) and precipitation (MAP).

2.2 Carbon flux data

2.2.1 Observations

Measured CO₂ fluxes come from Warm Winter 2020 (Warm Winter 2020 Team and ICOS Ecosystem Thematic Centre, 2022), an update of the FLUXNET2015 dataset (Pastorello et al., 2020) available on the ICOS platform (https://www.icos-cp.eu/data-products, last access: 6 January 2024). For each site, we selected daily time series of NEE (NEE_VUT_REF), GPP (GPP_DT_VUT_REF) and RECO (RECO_DT_VUT_REF), the latter two fluxes being derived from the daytime flux partitioning method (Lasslop et al., 2010). Preliminary analyses show the weak impact of the partitioning method (not shown).

The temporal coverage of the data varies by site (Table 1): less than 10 years for eight sites (SE-Svb, SE-Nor, SE-Htm, DE-HoH, DE-RuW, FR-Hes, FR-Bil and IT-SR2), between 10 and 20 years for five sites (FI-Let, DE-Hzd, CZ-BK1, FR-Fon and FR-FBn), and more than 20 years for six sites (FI-Hyy, DK-Sor, DE-Hai, DE-Tha, CH-Dav and IT-Ren). Given these limitations, the observational dataset likely lacks sufficient temporal depth to robustly assess the impact of climate on tower CO₂ flux interannual variability, highlighting the usefulness of models as complementary tools.

2.2.2 Data-driven models

Four data-driven models are used in this study (Table 2). The first data-driven model was developed by the CarbonSpace company to quantify carbon exchange at the site scale by integrating remote sensing data, meteorological variables and eddy covariance flux measurements. A Lagrangian particle dispersion model is used for footprint gas attribution. A machine learning model is used to solve the nonlinear regression problem of estimating fluxes from remote sensing and meteorological variables (Zhuravlev et al., 2022). For this study, the learning method is updated from the kernel method used in Zhuravlev et al. (2022) to an ensemble tree method (Chen and Guestrin, 2016). The key advantages of the CarbonSpace model include its scalability, high spatial resolution, and improved prediction accuracy through robust data quality control and advanced machine learning techniques. CarbonSpace provides monthly NEE only but at a very high spatial resolution (few hectares) from January 2000 to August 2023. This allows getting as close as possible to the 19 sites (around 1.8 ha centred on each tower) and their associated CO₂ flux measurement footprints.

The three other data-driven models come from the FLUXCOM initiative (Tramontana et al., 2016; Jung et al., 2019, 2020; Nelson et al., 2024). The first two are a three-member ensemble forced by both ERA5 reanalysis (Hersbach et al., 2020) and satellite data from the Moderate Resolution Imaging Spectroradiometer (MODIS) and a nine-member ensemble forced by MODIS only. The members differ by the machine learning method used to build each of the two models detailed in Jung et al. (2019). These models, named FLUXCOM-ERA5 and FLUXCOM-MODIS hereafter, provide global maps of monthly NEE, GPP and RECO derived with a daytime partitioning. The FLUXCOM-ERA5 model has a coarser resolution (0.5° × 0.5°) than FLUXCOM-MODIS (0.08° × 0.08°) but covers a longer period (1979–2018 versus 2001–2015). The last model, FLUXCOM-X, is a one-member model improving the coverage and quality of the training, as well as satellite data processing, and providing CO₂ fluxes at higher spatial resolution (0.05° × 0.05°) and for a longer period (2001–2020).

The four data-driven models include most, if not all, ICOS sites mobilized in this study. They accurately capture the mean annual and seasonal cycles of CO₂ fluxes (Tramontana et al., 2016; Jung et al., 2020; He et al., 2022; Zhuravlev et al., 2022) and are expected to outperform process-based models because the latter do not directly assimilate observed CO₂ fluxes. The methodological framework (e.g. machine learning model, forcing data and horizontal resolution) remains different between the data-driven models.

Table 2 Data-driven and process-based models used in this study.

Column 1 specifies the model type and name. Columns 2–4 detail the number of members in each dataset, the space-time resolution and the period of output availability. Columns 5–6 provide a brief description of each dataset along with references.

2.2.3 Process-based models

Two process-based models are considered (Table 2): the SMAP (Soil Moisture Active Passive) Level 4 Carbon model (SMAP-L4C hereafter) and an ensemble of dynamic global vegetation models (DGVMs) from the TRENDY project (https://sites.exeter.ac.uk/trendy, last access: 12 September 2024). The SMAP-L4C product is produced operationally by the NASA SMAP mission. It can be considered as a reanalysis product because it uses the Goddard Earth Observing System version 5 (GEOS-5) land model to assimilate SMAP L-band microwave observations and is forced with observed land cover and vegetation from the Moderate Resolution Imaging Spectroradiometer (MODIS) and Visible Infrared Imaging Radiometer Suite (VIIRS). The global processing is conducted on 1 km sub-grids using spatially aggregated MODIS plant functional types (PFTs) and VIIRS fraction of photosynthetically active radiation (fPAR) inputs, making it possible to distinguish up to eight individual PFTs within each 9 km × 9 km product grid cell. However, the model processing uses coarser-spatial-resolution (9 km and 0.25°) daily inputs from the SMAP L4 soil moisture (L4_SM) and GMAO Forward Processor (FP) surface meteorology. Among other variables, the SMAP-L4C outputs provide daily NEE and GPP (RECO deduced from the difference between NEE and GPP) in a consistent global grid from March 2015 to September 2023 for each PFT, including DBFs and ENFs (Jones et al., 2017; Kimball et al., 2022). The 1 km PFT subclass distinction allows the differentiation of ENFs and DBFs. The L4C product is derived using coupled photosynthetic light-use efficiency and soil organic matter decomposition models to estimate daily NEE and its component carbon fluxes, where GPP is reduced from PFT-specific optimal rates for unfavourable daily climate conditions, including cold temperatures, low light levels, excessive atmospheric vapour pressure deficits and low root zone (0–1 m depth) soil moisture levels defined from SMAP L4_SM and GMAO FP meteorology. The associated product quality assessment report gives details of the model algorithms and the calibration, validation and performance of the L4C version 7 product used in this study (Endsley et al., 2023).

In addition to SMAP-L4C, this study also uses outputs from 15 Dynamic Global Vegetation Models (DGVMs) of the Trends and Drivers of Regional-Scale Terrestrial Sources and Sinks of Carbon Dioxide (TRENDY version 12) project. These models are routinely mobilized to assess global carbon budget trends and for attributing changes to CO₂, climate and land use (Friedlingstein et al., 2023; Sitch et al., 2024). Appendix Table A1 provides a list of the 15 DGVMs used in this study. Here, we used outputs from the S3 scenario, with simulations starting in 1700 and forced by time-varying observed CO₂, climate and land-use change. All simulations have a horizontal resolution of 0.5° × 0.5° and monthly outputs.

2.2.4 Climate data

To investigate the impact of climate on CO₂ fluxes, we use the ERA5-Land dataset (Muñoz-Sabater et al., 2021) produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). This dataset results from the ECMWF land surface model (HTESSEL) operating at 0.1° spatial resolution and forced by the ERA5 reanalysis (Hersbach et al., 2020). This product provides hourly outputs for land surface, hydrological and meteorological variables from 1950 onwards. In this study, we use incident shortwave radiation, temperature at 2 m (T_2m) and averaged soil moisture (SM_AVG). We use the volumetric soil water content averaged across the four available soil layers (0–7, 7–28, 28–100 and 100–289 cm). Because it accounts for both liquid water and ice, this parameter remains above zero even when temperatures drop below freezing. Results obtained with incident shortwave radiation show no clear seasonality in correlation patterns with NEE, suggesting that greater light availability generally enhances CO₂ sequestration. For this reason, we do not include results for this variable in the main analysis. We also use relative humidity together with T_2m to compute the air vapour pressure deficit (VPD), an integrative metric accounting for both heat and water stress effects (Carrara et al., 2004; von Buttlar et al., 2018; Kong et al., 2022; van der Woude et al., 2023). The VPD is defined as the difference between the amount of moisture that is actually in the air and the amount of moisture that air could hold at saturation. The VPD is computed using the Tetens formula (Monteith and Unsworth, 2007) following Eq. (1):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathrm{VPD}=\left(\mathrm{1}-{\displaystyle \frac{\mathrm{HU}}{\mathrm{100}}}\right)\cdot \mathrm{saturation}\mathrm{vapor}\mathrm{pressure}=\\ & \left(\mathrm{1}-{\displaystyle \frac{\mathrm{HU}}{\mathrm{100}}}\right)\left(\mathrm{610.78}\cdot \exp \left({\displaystyle \frac{T}{T+\mathrm{237.3}}}\cdot \mathrm{17.2694}\right)\right).\end{array}\end{array}$$

In the end, three ERA5-Land climate variables (T_2m, SM_AVG and VPD) are used to assess the impact of climate on the annual cycle and interannual variability of CO₂ fluxes. These variables capture the influence of thermal, hydrological and atmospheric moisture demand conditions on CO₂ flux dynamics.

2.3 Methodology

For the gridded datasets (ERA5-Land, FLUXCOM, TRENDY and SMAP-L4C), we extracted the nearest grid point to each flux tower site. Note that SMAP-L4C simulates spurious CO₂ fluxes at the DE-Hzd site. Therefore, this site is not included in the analysis for this model. Because these datasets have varying temporal resolutions (Tables 1 and 2), all were aggregated to a monthly timescale. From these monthly values, we computed the mean annual cycle by averaging all available years in each dataset, along with interannual variability, defined by the standard deviation and coefficient of variation.

Model skill in capturing observed CO₂ flux variability is evaluated over overlapping periods between each model and observation. The number of overlapping years varies significantly across model–observation pairs (Fig. 2a). Two complementary metrics are used for model evaluation: the bias (model minus observation), which assesses errors in magnitude, and the Bravais–Pearson correlation coefficient (R), which evaluates temporal co-variability. These metrics capture distinct aspects of model performance and are not necessarily correlated. For the annual cycle, we computed monthly biases for each overlapping year and present the mean bias averaged across all months and years. Model skill in reproducing the seasonal timing of CO₂ fluxes is assessed by correlating the 12 monthly modelled and observed values within each overlapping year, with the multi-year mean R reported. Correlations are considered significant at the 95 % confidence level if the mean p value is below 0.05. For interannual variability, biases are calculated as the difference between modelled and observed standard deviations for each month. Co-variability between observed and modelled CO₂ fluxes is assessed only for model–observation pairs with at least 10 overlapping years, ensuring robust signal detection. Correlations were deemed significant when p<0.05.

The impact of climate on CO₂ fluxes is assessed for overlapping years between each CO₂ flux dataset and ERA5-Land. The number of overlapping years varies widely across datasets (Fig. 2b), ranging from low coverage in observations and SMAP-L4C to over 70 years in TRENDY. The correlation coefficient (R) was used to assess climate impacts, with the mean R reported for the annual cycle and individual R values for interannual variability.

To ensure results were not driven by long-term trends, analyses have been conducted using both raw and detrended climate time series (not shown), yielding similar outcomes. Additionally, for observed CO₂ fluxes, we verified ERA5-Land climate data reliability by comparing results with observed climate measurements from the FLUXNET database (not shown).

For conciseness, we primarily present results using the ensemble mean of the FLUXCOM (ERA5 and MODIS) and TRENDY models. However, uncertainties arising from machine learning methods and DGVM physical parameterizations are discussed in the model evaluation section.

Figure 2 Number of overlapping years between (a) observed and modelled CO₂ fluxes and (b) ERA5-Land climate and observed and modelled CO₂ fluxes. The overlapping periods in panel (a) are used to assess the model's ability to capture the annual cycle and interannual variability of CO₂ fluxes. In panel (b), the overlapping periods are used to evaluate the co-variability between CO₂ fluxes and the annual cycle and interannual variability of ERA5-Land climate. ENF sites are displayed in red text, while DBF sites are highlighted in green. Sites are ordered from north to south based on their latitude, with black vertical lines indicating the boundaries between northern (NE), central (CE) and southern (SE) Europe.

3 Results

3.1 Observed climate and CO₂ fluxes mean annual cycle and interannual variability

Figure 3 shows the mean annual cycle and interannual variability of T_2m, SM_AVG and VPD associated with each site. Overall, all sites depict higher T_2m and VPD and lower SM in summer than in winter (Fig. 3a–c). A south-to-north gradient is evident, with more marked annual amplitude and shorter summer in northern than in southern Europe. Few sites deviate from this pattern, including, for example, the Alpine site (CH-Dav), which depicts relatively cold and wet conditions, as well as low VPD, all year long. While the interannual variability of T_2m is the largest in winter regardless of the site, it increases markedly from south to north (Fig. 3d). The reverse is found for VPD, with higher interannual variability in summer than in winter, especially south of 60° N (Fig. 3f). The interannual variability in SM_AVG (Fig. 3e) is low all year long in the Alpine site, relatively low in northern Europe, and high from spring to summer in the mid-latitudes and in fall and winter in the Mediterranean region (FR-FBn).

Figure 3 (a–c) Mean annual cycle and (d–f) interannual variability (i.e. IV) in monthly 2 m air temperature (T_2m), soil moisture (SM_AVG) and vapour pressure deficit (VPD), respectively, for each study site for the 1979–2023 period. Climate conditions at each site are extracted from the nearest grid point of the 9 km × 9 km ERA5-Land product. ENF sites are displayed in red text, while DBF sites are highlighted in green. Sites are ordered from north to south based on their latitude, with black vertical lines indicating the boundaries between northern (NE), central (CE) and southern (SE) Europe.

Figure 4a displays the mean annual cycle of monthly NEE, GPP and RECO as provided by FLUXNET observations. The mean annual cycle in NEE is not necessarily phased on that of GPP and RECO, the two latter reaching their highest values from May to August in most sites. Significant differences are found between northern (SE-Nor, FI-Let, FI-Hyy and SE-Svb) and southern (FR-FBn, IT-SR2 and FR-Bil) Europe, where only ENF sites are available. The annual cycle is more marked, and the sequestration period (i.e. month associated with negative NEE values) is shorter in the former than in the latter region. Temperature conditions (and light availability) are the main drivers explaining these differences. In central Europe, where both ENF and DBF sites are available, there is a clear impact of land cover class. The DBF sites show a pronounced annual cycle with strong CO₂ uptake (below −3 tCO₂ ha⁻¹) from May to August, and up to some extent in September, and strong CO₂ release in the remaining months (above 2 tCO₂ ha⁻¹). Conversely, the ENF sites show a smoothed annual cycle: the summer peak of CO₂ sequestration barely exceeds −3 tCO₂ ha⁻¹, and the winter peak of CO₂ release rarely exceeds 1 tCO₂ ha⁻¹. Two sites deviate from the general pattern: the DE-Hzd DBF site, which acts as a CO₂ source nearly year-round, and the DE-RuW ENF site, which remains a consistent CO₂ sink on average.

The interannual variability of NEE, GPP and RECO, as defined by the standard deviation metric, tends to be stronger during the growing season (spring to fall) than in winter at almost all sites (Fig. 4b). The exact pattern depends on the CO₂ flux, site and land cover class considered. GPP is always close to zero during winter at DBF sites because trees are not photosynthetically active. This is not the case for ENF sites, particularly those located in central Europe (DE-RuW and CZ-BK1). The interannual variability of RECO is substantial in summer only in northern Europe. However, it can be non-negligible in other seasons in central and southern Europe. There, significant differences between geographically close sites (e.g. FR-Bil and IT-SR2) suggest additional drivers such as soil properties. The interannual variability of NEE is (i) weaker than that of GPP and RECO, likely due to the strong coupling between GPP and RECO, (ii) primarily driven by RECO in winter, and (iii) a complex response of GPP and RECO in the remaining seasons. Note that the pattern of CO₂ flux interannual variability depends on the metric used to assess it. When defined using the coefficient of variation, interannual variability is low in summer and high in winter for GPP and RECO, with increasing variability toward the north. For NEE, variability remains significant throughout the year, particularly from fall to winter (Fig. A1).

Overall, the mean annual cycle of observed CO₂ fluxes in European forests follows a clear spatial pattern driven by climate conditions and land cover class. In contrast, the interannual variability of observed CO₂ fluxes exhibits greater spatial noise across Europe and depends on the metric used. Nevertheless, it remains significant, highlighting the importance of assessing the impact of climate on it throughout the annual cycle.

Figure 4 (a) Mean annual cycle and (b) interannual variability in monthly NEE, GPP and RECO (from left to right, respectively) as measured by eddy covariance over the most extended available period for each site (see Table 1 for details). ENF sites are displayed in red text, while DBF sites are highlighted in green. Sites are ordered from north to south based on their latitude, with black vertical lines indicating the boundaries between northern (NE), central (CE) and southern (SE) Europe.

3.2 Model evaluation in capturing the mean annual cycle and interannual variability of CO₂ fluxes

3.2.1 Mean annual cycle

The model skill in capturing the temporal phasing of the annual cycle and the magnitude in observed CO₂ fluxes is assessed in terms of correlation and mean bias, respectively (see Sect. 2.3 for details). All models accurately capture the observed temporal phasing of GPP and RECO, with correlation values often above 0.8 (Fig. 5a). The model skill is poorer but still correct for NEE, with correlation values remaining above 0.6 for most sites and models. The weaker correlation found for NEE compared to GPP and RECO is not surprising, as accurately estimating the NEE annual cycle requires precise estimation of both the temporal phasing and magnitude of GPP and RECO.

Figure 5 Model skill in capturing (a) the annual cycle and (b) magnitude of observed NEE, GPP and RECO. The model skill in capturing the annual cycle is assessed through the Bravais–Pearson correlation coefficient (R) calculated between the simulated and observed annual cycle of monthly CO₂ fluxes (12 values) for each year of the overlapping period. The multi-year mean R values are shown. The magnitude error is computed as the difference between simulated and observed CO₂ fluxes (model minus observation) for each month of the overlapping period. The mean magnitude error is shown. Results are displayed using coloured dots for single-member models (CarbonSpace, FLUXCOM-X and SMAP-L4C) and boxplots for multi-member FLUXCOM models and the multi-model TRENDY ensemble. The boxes have lines at the lower quartile, median and upper quartile values. The whiskers are lines extending from each end of the boxes to show the extent of the full range of the data, including outliers. The colour attributed to each model is detailed in the legend. ENF sites are displayed in red text, while DBF sites are highlighted in green. Sites are ordered from north to south based on their latitude, with black vertical lines indicating the boundaries between northern (NE), central (CE) and southern (SE) Europe.

Despite reasonable annual cycles, the models struggle in capturing the observed magnitude of RECO and GPP (Fig. 5b). Three groups emerge. The first group includes the FLUXCOM (ERA5, MODIS and X) and TRENDY models, which underestimate both GPP and RECO by about similar amounts, resulting in relatively “weak” positive biases in NEE. The second group corresponds to the SMAP-L4C model, which overestimates RECO by 1 tCO₂ ha⁻¹ month⁻¹ while underestimating GPP by the same amount, leading to a systematic underestimation of CO₂ sequestration by approximately 2 tCO₂ ha⁻¹ month⁻¹. The last group is the CarbonSpace data-driven model, which is the only model that systematically overestimates CO₂ sequestration, by up to 2.5 tCO₂ ha⁻¹ month⁻¹. However, the cause of this overestimation is unclear, as this model does not provide separate GPP and RECO estimates.

Figure 5 highlights key insights into model behaviour. First, the data-driven and SMAP-L4C models generally outperform the TRENDY models in capturing the annual cycle of CO₂ fluxes but do not necessarily provide better estimates of flux magnitude. Second, models that accurately represent the annual cycle can still struggle with magnitude. This is exemplified by the CarbonSpace data-driven model, which ranks among the best for annual cycle representation but severely overestimates CO₂ sequestration at many sites. Third, the machine learning methods used in the FLUXCOM-ERA and FLUXCOM-MODIS ensembles have little impact on both the annual cycle and magnitude of CO₂ fluxes. The input data appear to be more important, with the FLUXCOM models accounting for both vegetation and climate (i.e. ERA5 and X) yielding more reliable results than those accounting for vegetation alone (MODIS). Additionally, FLUXCOM-X improves upon the previous model generation for most sites. Finally, the inter-model spread within the TRENDY ensemble is much smaller for CO₂ flux magnitude than for temporal variability, suggesting that the primary source of uncertainty in DGVMs lies in the temporal phasing of the fluxes rather than their magnitude.

3.2.2 Interannual variability

We qualitatively evaluate how well the models capture the observed interannual variability of CO₂ fluxes in terms of magnitude (monthly bias analysis) and temporal phasing (correlation analysis) over the overlapping period of each model–observation pair.

All models strongly underestimate the magnitude in NEE interannual variability all year long (Fig. 6), particularly during summer, where biases often exceed 2.5 tCO₂ ha⁻¹ month⁻¹. This is a well-known bias of current data-driven and process-based models (e.g. Lin et al., 2023; Nelson et al., 2024). The only exception is the CarbonSpace model, which produces weak positive or negative biases at most sites. Importantly, a biased magnitude of CO₂ flux interannual variability (as measured by the standard deviation) does not preclude the models to capture their temporal co-variability (as measured by correlation) with observed CO₂ fluxes and climate.

Figure 6 Model skill in capturing the observed magnitude in monthly NEE interannual variability. Model skill is assessed by computing biases between modelled and observed NEE interannual variability, as defined by the standard deviation of monthly fluxes. The number of years included varies depending on the model–observation pair, as shown in Fig. 2a. Cells labelled “NA” indicate cases where the overlap criterion is not met. For conciseness, biases are computed using the ensemble mean for the FLUXCOM-ERA5, FLUXCOM-MODIS and TRENDY models. ENF sites are displayed in red text, while DBF sites are highlighted in green. Sites are ordered from north to south based on their latitude, with black vertical lines indicating the boundaries between northern (NE), central (CE) and southern (SE) Europe.

Figure 7 shows the correlations between the modelled and observed interannual variability in monthly NEE. Correlation values are predominantly positive across Europe, though they are often low and not statistically significant (p>0.05). The correlation values tend to be higher in fall for northern Europe sites and all year long in central Europe regardless of the model. The frequent lack of statistical significance in correlations can largely be attributed to the limited number of overlapping years. This is further supported by the fact that, except for Davos, correlation values tend to be higher and more likely to reach the 95 % confidence level in model–observation pairs with a greater number of overlapping years (e.g. DK-Sor, DE-Hai and DE-Tha). Another challenge in capturing the observed temporal phasing of NEE interannual variability is that it requires accurately simulating the interannual variability of both GPP and RECO. The latter is generally better represented by models than the NEE variability itself (compare Fig. 7 with A2–A3). Scale inconsistencies may also contribute to the discrepancies. While flux tower observations reflect local variability, most models represent regional-scale fluxes (and drivers). This is supported by the CarbonSpace model, the only site-scale model used in this study, which produces more satisfactory results, with positive correlation values for almost all sites and all months.

Figure 7 Model skill in capturing the observed temporal phasing in monthly NEE interannual variability. Model skill is assessed using the Bravais–Pearson correlation coefficient (R), calculated for each month where at least 10 years of overlap exist between each model–observation pair. The number of years included varies depending on the model–observation pair, as shown in Fig. 2a. Correlation values marked with “×” are significant at the 95 % confidence level according to the Bravais–Pearson test. Cells labelled “NA” indicate cases where the overlap criterion is not met. For conciseness, the correlation analysis is performed using the ensemble mean for the FLUXCOM-ERA5, FLUXCOM-MODIS and TRENDY models. ENF sites are displayed in red text, while DBF sites are highlighted in green. Sites are ordered from north to south based on their latitude, with black vertical lines indicating the boundaries between northern (NE), central (CE) and southern (SE) Europe.

The spread among members of the FLUXCOM ensembles is low (not shown), regardless of whether they account for climate alone or both climate and vegetation. However, Fig. 7 reveals significant differences between the FLUXCOM products (ERA5, MODIS and X; Fig. 7b–d). These differences are not necessarily due to the type of data used in these data-driven models, as the analysis periods differ between them. In contrast, the correlations between observations and the TRENDY ensemble mean (Fig. 7e) mask substantial variability among individual TRENDY models, with no single model consistently outperforming the others (Fig. A4).

These qualitative results suggest that the interannual variability of simulated NEE is at least partially aligned with that of observed CO₂ fluxes, supporting the use of models to assess the impact of climate on CO₂ flux interannual variability.

3.3 Climate–CO₂ flux relationship

3.3.1 Annual cycle

Figure 8 assesses the co-variability between the annual cycle of CO₂ fluxes and climate variables (T_2m, SM_AVG and VPD) through correlation analysis. For conciseness, correlations are computed using the ensemble mean of CO₂ fluxes for the FLUXCOM (ERA5 and MODIS) and TRENDY models. The inter-member dispersion in FLUXCOM and inter-model dispersion in TRENDY, shown in Fig. A5, are similar to those in Fig. 5 and are therefore not discussed here.

The annual cycle of CO₂ fluxes is closely linked to climate in both observations and models, as evidenced by statistically significant correlations at the 95 % confidence level. The influence of climate on the annual cycle of CO₂ fluxes is generally stronger for GPP and RECO than for NEE, particularly for T_2m (Fig. 8a). Over Europe, GPP, RECO and CO₂ sequestration (i.e. negative NEE) tend to be higher when T_2m and VPD are high and SM_AVG is low. In turn, CO₂ fluxes are amplified in summer and dampened in winter. Among the climate variables, T_2m and VPD exhibit stronger correlations with the CO₂ flux annual cycle than SM_AVG, reflecting their more pronounced seasonality (Fig. 3).

This climate–CO₂ flux relationship is particularly evident in northern and central Europe, where there is strong agreement between observations and models, as well as across different models. The only exception in central Europe is the Alpine site (CH-Dav), where SM_AVG is positively correlated with GPP and RECO and negatively correlated with NEE. This discrepancy arises because CH-Dav is the only site in the study where SM_AVG is higher in summer than in winter (Fig. 3b). In contrast, in southern Europe, the relationship between climate and NEE does not reach the 95 % confidence level in observations. However, the models show a similar relationship to that found in northern and central Europe, albeit with weaker correlations. Whether this disagreement stems from the limited number of observational years available for these sites (7 to 13 years, see Table 1) or from an overestimation of climate impact on the CO₂ flux annual cycle remains an open question.

Figure 8 Co-variability in the annual cycle of CO₂ fluxes and (a) 2 m temperature, (b) soil moisture and (c) vapour pressure deficit across all sites and datasets. Co-variability is assessed through the Bravais–Pearson correlation, computed annually between CO₂ fluxes and ERA5-Land climate data, then averaged across years. The number of years included varies depending on the data, as shown in Fig. 2d. Correlation values marked with “×” are significant at the 95 % confidence level according to the Bravais–Pearson test. Cells labelled “NA” indicate cases where GPP and RECO are unavailable (CarbonSpace) or where CO₂ flux data are corrupted (SMAP-L4C). For conciseness, correlations are computed using the ensemble mean of CO₂ fluxes for the FLUXCOM (ERA5 and MODIS) and TRENDY models. ENF sites are displayed in red text, while DBF sites are highlighted in green. Sites are ordered from north to south based on their latitude, with black vertical lines indicating the boundaries between northern (NE), central (CE) and southern (SE) Europe.

Figure 8 also reveals differences in model behaviour. First, the climate–CO₂ flux relationship appears sensitive to the land cover class only in data-driven models. In the FLUXCOM models (ERA5, MODIS and X), this relationship is weaker for DBFs than for ENFs in central Europe, whereas the CarbonSpace model shows the opposite pattern. This sensitivity to land cover class may reflect model discrepancies, as the observed climate–CO₂ flux relationship does not depict such an ENF–DBF distinction. The reason for this discrepancy may involve differences in climate forcing or model spatial resolution. Second, process-based models tend to overestimate the impact of climate on the annual cycle of CO₂ fluxes across all sites. For example, the SMAP-L4C model amplifies the influence of T_2m and SM_AVG on NEE, particularly in northern Europe (Fig. 8a–b). Meanwhile, in the TRENDY models, VPD emerges as the dominant driver of the NEE annual cycle (Fig. 8c). These results suggest that process-based models may underestimate the role of additional factors, such as soil properties, in shaping the CO₂ flux annual cycle.

3.3.2 Interannual variability

We now investigate the interannual co-variability between CO₂ fluxes and climate using correlation analysis. This is done by analysing the full monthly time series of each dataset in two ways: (i) all months combined after removing the mean annual cycle and (ii) separately for each calendar month. The first approach reveals no robust relationship in the observations and weak correlations in the models (Fig. 9). In contrast, the second approach shows some consistencies in the correlation patterns across datasets, with distinct differences between northern/central Europe and southern Europe (Fig. 10). In the former region, the NEE tends to be positively correlated with T_2m, SM_AVG and VPD in winter and fall. This means that anomalously high T_2m, VPD and SM_AVG favour CO₂ release during the cold seasons. These patterns are similar for GPP and RECO (Figs. A6–A7), but correlations are stronger for RECO, highlighting its significant contribution to NEE interannual variability during the cold seasons. In spring (March–May), the relationship between NEE and $T_{\mathrm{2}\mathrm{m}}/\mathrm{VPD}$ reverses compared to the cold season pattern, while that between NEE and SM_AVG tends to remain positive in central Europe. In turn, anomalously high T_2m and VPD, along with anomalously low SM_AVG (at least in central Europe), tend to be favourable spring conditions for CO₂ sequestration. The increase in spring CO₂ sequestration under anomalously dry soil conditions is likely driven by elevated T_2m and VPD, leading to enhanced evapotranspiration and drier soils. The summer pattern generally shows an opposite sign compared to spring, suggesting stronger CO₂ sequestration under anomalously low T_2M and VPD, and high SM_AVG. In southern Europe, correlation patterns are consistent across datasets in spring and summer only, showing a similar response to that observed in summer in northern/central Europe but extended over a longer period.

Figure 9 Same as Fig. 8 but using monthly anomalies calculated after removing the mean annual cycle.

Figure 10 Interannual co-variability between NEE and (a) 2 m temperature, (b) soil moisture and (c) vapour pressure deficit for each month, all sites and all datasets. The co-variability is assessed through the Bravais–Pearson correlation computed for each month between CO₂ fluxes and ERA5-Land climate data. The number of years included varies depending on the data, as shown in Fig. 2d. Correlation values marked with “×” are significant at the 95 % confidence level according to the Bravais–Pearson test. Cells labelled “NA” indicate cases where the CO₂ flux data are corrupted (SMAP-L4C). For conciseness, correlations are computed using the ensemble mean of CO₂ fluxes for the FLUXCOM (ERA5 and MODIS) and TRENDY models. ENF sites are displayed in red text, while DBF sites are highlighted in green. Sites are ordered from north to south based on their latitude, with black vertical lines indicating the boundaries between northern (NE), central (CE) and southern (SE) Europe.

4 Discussion

This study aims at assessing the impact of climate on the annual cycle and interannual variability of monthly CO₂ fluxes in European DBFs and ENFs through conjointly analysing observations from the FLUXNET network and state-of-the-art data-driven and process-based models.

As a first step, we assess the model abilities to reproduce the observed mean annual cycle and interannual variability of CO₂ fluxes. This evaluation presents two key challenges. First, the temporal coverage of observations in the FLUXNET database is often limited, making it difficult to extract robust signals, particularly for interannual variability. Site-specific characteristics may also cause a disconnect between CO₂ flux variability and regional climate variability (Chu et al., 2017). Second, there is a spatial scale mismatch between site-level observations, representing fluxes from the tower footprint to several square kilometres (Göckede et al., 2008), and most models used in this study, which simulate fluxes at regional to large scales, except for the hectometric-scale CarbonSpace model. Given these constraints, our evaluation should be considered qualitative rather than strictly quantitative.

The models reasonably capture both the annual cycle and interannual variability of observed CO₂ fluxes, though with some magnitude discrepancies. In particular, the models accurately capture the north-to-south gradient, with increased length in the CO₂ sequestration period from northern to southern Europe, as well as the more pronounced annual cycle in DBFs than ENFs. With the exception of the CarbonSpace model, which overestimates CO₂ sequestration across all European sites, most models tend to underestimate GPP and overestimate RECO, resulting in a near systematic underestimation of CO₂ sequestration.

The interannual variability in the models is weaker than in the observations, consistent with previous studies (e.g. Nelson et al., 2024). However, the CarbonSpace data-driven model proved to be the only model tested that does not underestimate the NEE interannual variability. The reasons may involve its high spatial resolution (few hectares) and the use of a Lagrangian particle dispersion model, which allows it to closely align with the flux tower footprints. This results in more precise flux localization, which may improve its response to fine-scale variability. Another reason may also involve the use of an ensemble tree method for regression. This method offers greater flexibility in capturing nonlinear interactions between environmental variables and NEE. Further studies are needed to evaluate these hypotheses. The temporal co-variability between observed and simulated CO₂ fluxes remains correct despite the underestimated magnitude in CO₂ flux interannual variability. This agreement was expected for the data-driven models, as they incorporate FLUXNET observations in their development. However, it was less anticipated for the process-based models, which do not assimilate direct CO₂ flux measurements. Their ability to capture observed interannual variability likely stems from the fact that TRENDY models are driven by observed CO₂ concentrations, land-use changes and climate data, while the SMAP-L4C model benefits from the assimilation of satellite-derived soil moisture observations.

Despite uncertainties, our results highlight that state-of-the-art models are valuable tools for complementing observations, particularly in assessing the impact of climate on the interannual variability of CO₂ fluxes at the European scale. The relationship between CO₂ fluxes and climate is analysed for both the annual cycle and interannual variability using synchronous correlation analyses between each CO₂ flux and individual climate variables, including 2 m temperature (T_2m), vapour pressure deficit (VPD) and soil moisture (SM_AVG).

Regarding the annual cycle, the influence of climate on CO₂ fluxes is stronger in northern and central Europe than in southern Europe, where seasonal climate variations are less pronounced. In southern Europe, both the agreement between observations and models and the consistency among different data-driven and process-based models are weaker. The spread within members (FLUXCOM-ERA5 and FLUXCOM-MODIS) and DGVMs (TRENDY) is also larger. We hypothesize that part of this model uncertainty stems from the limited number of ENF sites in southern Europe, leading to weak constraints for data-driven models and fewer reference points for calibrating DGVMs. However, it is worth noting that the uncertainty associated with the machine learning methods used to develop data-driven models remains low, whereas inter-DGVM spread can be substantial.

Regarding interannual variability, the climate impact on CO₂ fluxes can be summarized in three points. First, climate impacts GPP and RECO more strongly than NEE, regardless of the site and dataset. Because NEE is the difference between GPP and RECO, its interannual variability arises from various combinations of these two components. For instance, reduced CO₂ sequestration can result from a greater increase in RECO compared to GPP, a larger decrease in GPP than in RECO, a decrease in GPP with no change in RECO or an increase in RECO without any change in GPP. Such complexity implies that the influence of climate on NEE is less direct and likely more intricate than its effects on GPP and RECO. Second, the impact of climate on the interannual variability of CO₂ fluxes depends strongly on how it is assessed. It appears weak when monthly anomalies from all months are analysed together but becomes more pronounced when each month is examined separately, revealing seasonal shifts in the sign of the climate–NEE relationship. The timing and direction of these seasonal shifts vary across datasets and regions. In particular, CO₂ sequestration is enhanced by anomalously low T_2m and VPD and anomalously high SM_AVG in summer, whereas the opposite pattern prevails in spring. This finding aligns with previous case studies showing that heatwaves and droughts reduce summer CO₂ sequestration (Ciais et al., 2005; Smith et al., 2020; Thompson et al., 2020; van der Woude et al., 2023). It also highlights that anomalously warm and dry springs may, in some cases, enhance CO₂ sequestration, likely because soil moisture levels remain sufficient during this period, in line with, for example, Delpierre et al. (2009) and Smith et al. (2020). Such a transition from spring to summer is less evident in southern Europe, which instead exhibits consistent patterns from spring to summer. The climate–NEE relationship is much noisier in both space and time in the observations than in the models, and it can vary substantially across different models. This indicates that local flux measurements may not reliably represent regional-scale dynamics, while models may exaggerate the influence of climate on CO₂ flux variability (despite underestimating its magnitude). Further work is needed to disentangle site-specific effects from broader-scale signals, a critical step toward improving the calibration of regional and global models that cannot resolve local heterogeneity.

5 Conclusions

This study makes use of state-of-the-art data-driven and process-based models to complement observations for assessing the impact of climate on the annual cycle and interannual variability of monthly CO₂ fluxes in European DBFs and ENFs. Outputs from different data-driven models (CarbonSpace, FLUXCOM-ERA5, FLUXCOM-MODIS and FLUXCOM-X) and process-based models forced by realistic conditions (assimilation of satellite-derived soil moisture in SMAP-L4C and time-varying CO₂ concentration, land use and climate in TRENDY DGVMs) are analysed conjointly with CO₂ measurements from 19 sites (6 DBFs and 13 ENFs) of the FLUXNET network. Across Europe, a clear north-to-south gradient emerges in the annual cycle of CO₂ fluxes. The length of the CO₂ sequestration season increases southwards in conjunction with more favourable climate conditions for photosynthesis in southern than in northern Europe. This large-scale pattern is perturbed locally by site elevation and other site factors not included in the study (e.g. soil properties and forest age). It is also perturbed by the land cover class, with a more pronounced annual cycle of DBFs than ENFs in central Europe. The interannual variability of CO₂ fluxes does not exhibit such a north-to-south gradient, regardless of the metric used (standard deviation or coefficient of variation). However, it remains strong across all seasons, with spring and summer showing high variability based on standard deviation and autumn and winter showing high variability based on the coefficient of variation.

The models accurately capture the observed features despite magnitude differences. Compared to observations, the CO₂ sequestration is weaker in regional-scale models and stronger in the hectometric-scale CarbonSpace data-driven model. Except for the CarbonSpace model, all models systematically underestimate the interannual variability of CO₂ fluxes, as already reported (Lin et al., 2023; Nelson et al., 2024). Despite biased magnitude, the interannual variability of modelled fluxes correlates well with the observations. This supports the use of models to complement observations, whose limited temporal coverage and site specificities hinder the assessment of climate impacts on CO₂ interannual variability. We show that the influence of climate on CO₂ flux interannual variability is obscured when monthly anomalies are analysed together. This apparent lack of relationship masks distinct seasonal patterns, which are concealed when considering all months together. Winter and fall CO₂ release increases under elevated temperature and VPD in northern and central Europe, while no clear signal emerges in southern Europe. The CO₂ sequestration increases under anomalously hot and dry conditions in spring and cold and wet conditions in summer in northern/central Europe. Anomalously cold and wet conditions also favour CO₂ sequestration in southern Europe from spring to summer. While these seasonal signals appear noisy in the observations, due to limited sample sizes and site-specific variability, they emerge more clearly in the models, albeit with some model-dependent differences.

These results highlight the significant space-time variability in the impact of climate on forest CO₂ fluxes across Europe. This variability underscores the importance of considering regional and seasonal differences when assessing the effects of climate change on CO₂ fluxes. Neglecting these variations could lead to oversimplified conclusions and hinder the development of accurate predictions and effective mitigation strategies.

Appendix A

Table A1 List of TRENDY DGVMs used in this study.

Columns 1 and 2 provide the name of the DGVM and its corresponding reference, respectively.

Figure A1 Same as Fig. 4b but using the coefficient of variation as a metric for interannual variability.

Figure A2 Same as Fig. 7 but for GPP.

Figure A3 Same as Fig. 7 but for RECO.

Figure A4 Same as Fig. 7e but for each TRENDY model considered in this study.

Figure A5 Same as Fig. 5 but for the relationship between the annual cycle of CO₂ fluxes and (a) 2 m temperature, (b) soil moisture and (c) vapour pressure deficit. This relationship is evaluated using the Bravais–Pearson correlation. For details, refer to the legend caption of Fig. 8.

Figure A6 Same as Fig. 10 but for GPP.

Figure A7 Same as Fig. 10 but for RECO.