Systematic underestimation of type-specific ecosystem process variability in the Community Land Model v5 over Europe

Evapotranspiration (ET) and gross primary production (GPP) are critical fluxes contributing to the energy, water, and carbon exchanges between the atmosphere and the land surface. Land surface models such as the Community Land Model v5 (CLM5) quantify these fluxes, estimate the state of carbon budgets and water resources, and contribute to a better understanding of climate change's impact on ecosystems. Past studies have shown the ability of CLM5 to model ET and GPP magnitudes well but emphasized systematic underestimations and lower variability than in the observations.

Here, we evaluated CLM5's predictions of water and energy fluxes using observations from eddy covariance stations from the Integrated Carbon Observation System (ICOS), remote sensing, and reanalysis data sets. We assess simulated ET and GPP from the grid scale (CLM5_grid) and the plant functional type (PFT) scale (CLM5_PFT). CLM5_PFT exhibited a low systematic error in simulating the ET at the ICOS sites (average bias of −4.68 %), indicating that PFT-specific ET closely matches the observations' magnitude. GPP was underestimated by CLM5_PFT, especially in deciduous forests (bias of −43.76 %). The results showed an underestimation of the spatiotemporal variability in the simulated ET and GPP distribution moments across PFTs for both CLM setups compared to reanalysis data and remote-sensing products. These findings provide essential insights for improving land surface models, highlighting the need to enhance the CLM5's ability to capture the spatiotemporal variability in ET and GPP simulations across PFTs.

1 Introduction

Ecosystem processes, such as evapotranspiration (ET) and gross primary production (GPP), play an important role in cycling water, carbon, and energy between ecosystems and the atmosphere. Changes in the magnitude and variability in these fluxes can indicate the ecosystems' inhibited performance due to changing environments (Kühn et al., 2021; Migliavacca et al., 2021). These changes can lead to short-term alterations and long-term trends in water resources and carbon pools in the atmosphere and the land surface. Thus, the accurate quantification of the variability in ecosystem processes is pivotal for developing climate change projections and formulating effective mitigation policies (Friedlingstein et al., 2023; Graf et al., 2023).

Notably, an accurate and functional understanding of land surface processes is essential to identify threatened ecosystems in the present and the future and facilitate carbon budget calculations. Land surface models (LSMs) serve as deterministic and process-based simulators of ecosystems, capturing energy, water, and carbon fluxes, while considering their interactions and the heterogeneity of the land surface (Fisher and Koven, 2020). LSMs can complement point scale observations from in situ research infrastructures by providing spatiotemporally uniform and extensive high-resolution outputs. Their high-resolution process-based simulations contrast the often coarsely resolved remote-sensing data. Hence, LSMs are frequently used tools for investigating and projecting the current understanding of ecosystem processes, such as GPP and ET, on various scales. However, there is uncertainty in the LSM structure, the parameters, the input data, and the initial conditions, which carry over to the simulated variables. Therefore, assessing how well the general simulated ET and GPP variability compares to the observations is crucial. Such evaluations deliver an essential context on LSM biases and form a basis for analyses of more complex ecosystem responses. Recent studies already found discrepancies between LSM simulations of ET and GPP and observations collected in the field and from remote sensing. For instance, these discrepancies are evident in their magnitude and variability (De Pue et al., 2023; Boas et al., 2023; Cheng et al., 2021; Strebel et al., 2024) and their response to drought (Ukkola et al., 2016; Wu et al., 2020; Green et al., 2024). Therefore, assessing the accuracy of LSMs in representing observed GPP and ET fluxes is crucial to test and improve our current understanding of the ecosystem process variability and identify the limitations of state-of-the-art LSMs.

Current land surface models, e.g., the Joint UK Land Environment Simulator (JULES), the Community Land Model 5 (CLM5), or the Community Atmosphere Biosphere Land Exchange Model (CABLE), employ a tiling system within the grid cell to account for functional differences in distinct patches on the land surface. The natural and crop vegetation is grouped into plant functional types (PFTs), the entities for which ecosystem process calculations are resolved (Fisher and Koven, 2020; Bonan et al., 2002; Solomon and Shugart, 1993). Typically, PFTs are defined based on morphological and phenological characteristics of the vegetation (e.g., leaf type and leaf longevity) and climate (Bonan et al., 2002). However, the usefulness of this PFT definition, or at least its current coarsely resolved implementation, is a subject of debate (Caldararu et al., 2016; Van Bodegom et al., 2012). The primary argument against it is that observed plant traits implemented as PFT-related parameters vary to some extent in space and time in response to a changing environment. This spatiotemporal dependence of PFT traits is only marginally represented in LSMs. On top of that, most research assessing LSMs only used a handful of observation sites and did not analyze aggregated values for groups of sites observing the same PFT. Such analyses would provide essential insights; a recent study highlighted the differences between vegetation-type concepts used in observation networks, e.g., the International Geosphere–Biosphere Programme (IGBP) classification, and PFTs used in LSMs and underlined the importance of improving these PFT concepts (Cranko Page et al., 2024).

The phenology of ecosystem processes, i.e., their seasonal cycles and evolution through the year and the growing season length, have shifted in timing due to climate change. A recent study investigated which factors drive the changes in the mean annual dynamics of ecosystem processes in Europe (Rahmati et al., 2023), and many of these discovered feedbacks, for instance, the effect of increased atmospheric dryness on growing season length, are only implemented simplistically in LSMs. Furthermore, robust simulations of LSMs for impact assessments become even more critical as ecosystems experience more disturbances, along with the changing climate. For example, projections show that droughts have recently become more frequent in Europe (Vautard et al., 2023; Rousi et al., 2022) and that these extreme events will become even more frequent and severe in the future (Lehner et al., 2017). While the combined effect of a higher occurrence of compound drought events is currently not fully understood, it is clear from observations that individual drought years, or droughts in general, have already had a profound impact on ecosystem processes in Europe (Graf et al., 2020; Van Der Woude et al., 2023; Poppe Terán et al., 2023). Given that the frequency and severity of extreme events affect GPP and the ET statistical distributions, investigating how the characteristics of the simulated distributions compare with the observed can contextualize findings of modeled ecosystem drought responses in Europe.

One predominantly used LSM is the Community Land Model version 5 (CLM5) (Lawrence et al., 2019, 2018). In the most recent version, CLM5 solves the biogeochemistry (BGC), i.e., the carbon and nitrogen cycles between the atmosphere, vegetation, and soil. CLM5 has been widely employed for quantifying and examining ecosystems at various scales, including global (Xie et al., 2020; Sitch et al., 2015; Lawrence et al., 2019), regional (Cheng et al., 2021; Boas et al., 2023), and site-scale (Strebel et al., 2024; Umair et al., 2020; Song et al., 2020; Fisher et al., 2019) applications. Several studies have highlighted the ability of CLM5 to simulate ecosystem processes close to the observations (Wozniak et al., 2020; Lawrence et al., 2019; Cheng et al., 2021; Zhang et al., 2023; Boas et al., 2023). However, they have also emphasized an underestimated magnitude and variability in the simulations across different timescales and under various conditions.

The present study assesses the CLM5's ability to capture ecosystem processes at a continental scale. To ensure comparability to point scale observations, we conducted high-resolution simulations at 0.0275° (approx. 3 km) resolution over the European Coordinated Regional Climate Downscaling Experiment (CORDEX) domain (Giorgi et al., 2009), resulting in 1544×1592 grid cells. Notably, the output contained variables from the subgrid scale, i.e., from within a 3 km grid cell, for PFTs present in the grid cell. We then compared the CLM5 grid level (CLM5_grid) and PFT level data (CLM5_PFT) to observations from a continental network of sites. The Integrated Carbon Observation System (ICOS) provides the Warm Winter 2020 data (Warm Winter 2020 Team and ICOS Ecosystem Thematic Centre, 2022), which include eddy covariance measurements over a dense network of over 70 sites in Europe. It was named after and curated to support research on the effect of the warm winter of 2020 on terrestrial carbon fluxes. These ICOS data are regarded as the gold standard for calibrating and evaluating process-based models due to their ample spatial coverage as a network encompassing diverse land cover types. Thus, it offers an excellent opportunity to comprehensively assess simulated GPP and ET for specific PFTs from our CLM5 setup over Europe.

Additionally, we include remote-sensing data from the Global Land Surface Satellite (GLASS; Liang et al., 2021) and reanalyses from the European Center for Medium-Range Weather Forecasts Reanalysis 5 – Land (ERA5-Land; hereafter ERA5L; Muñoz Sabater, 2019), as well as from the Global Land Evaporation Amsterdam Model (GLEAM; Martens et al., 2017) in our analyses to identify common patterns of ecosystem process variability between CLM5, in situ observations, reanalysis, and remote-sensing data.

In summary, this study uses ICOS observations as ground truth data. It compares them with grid level and PFT level CLM5 data and terrestrial surface fluxes from reanalyses and remote-sensing derivatives to the following:

1. Compare performance indices (root mean square error and percent bias) between the models and ICOS measurements on a per-site and PFT group basis to assess the systematic error and accuracy of ET and GPP simulations.

2. Investigate how the models represent the observed ET and GPP for different PFTs regarding their sub-annual averaged phenologies, standard deviation, and timing of important phenological events.

3. Evaluate the simulated PFT level ET and GPP statistical distributions and their moments (mean, variance, skewness, and excess kurtosis) to contextualize assessments of factors like droughts, which impact the shape of these distributions.

4. Compare the inter-site differences between ET and GPP time series within PFT groups to estimate how the observed intra-PFT variability is represented in the models.

Thus, these findings offer critical information for comparisons of GPP and ET from the evaluated models. Furthermore, this study also paves the way for a better-informed analysis of the drought response of ET and GPP from the models being assessed over Europe. We expect the following:

1. There is a lower systematic bias, and the simulation is closer to the observations by the PFT scale than the grid scale CLM5 outputs and remote-sensing and reanalysis data.

2. The remotely sensed and modeled data approximate critical events in the phenologies of ET and GPP within the standard deviation of the ICOS measurements for sites of one PFT. However, this ability varies between PFTs.

3. The remotely sensed and modeled ET and GPP data distributions show a lower range among the moments within the PFT groups than the ICOS measurements.

2 Methods and data

2.1 Community Land Model version 5

We use the CLM5 (Lawrence et al., 2018, 2019), which is forced offline with custom input data. The land surface of a region in CLM5 is first disaggregated into grid cells, which are uniformly distributed and simulated individually. These grid cells are tiled into land units (i.e., natural vegetation, crops, lakes, urban areas, and glaciers) with a relative area coverage within the grid cell. Importantly, plants in the naturally vegetated land units compete for water in a single soil column. The vegetation is grouped into PFTs (Lawrence and Chase, 2007), which are distinguished through leaf habit (evergreen or deciduous), morphology (needle- and broadleaves, grass, and shrubs), and the bioclimate of the grid cell location (boreal, temperate, and tropical). While competition for soil moisture includes interactions among different PFTs, this is closer to natural conditions than separated soil columns and encourages evaluations on the PFT scale. Here, we use CLM5–BGC, which calculates vertical carbon and nitrogen pools and fluxes between the vegetation, soil, and atmosphere. In the following subsections, we briefly describe the essential processes in CLM5 that are particularly relevant to this study, as well as the input data and leading features of the European CLM5 setup.

2.1.1 Gross primary production and evapotranspiration

The stomatal conductance of plants (g_s) couples the water exchange with carbon uptake between vegetation and the atmosphere. In the CLM5, g_s is calculated by the Medlyn stomatal conductance model (Medlyn et al., 2011) as follows:

$$\begin{array}{}\text{(1)}& g_{\mathrm{s}}=g_{\mathrm{0}}+\mathrm{1.6}\left(\mathrm{1}+{\displaystyle \frac{g_{\mathrm{1}}}{\sqrt{D}}}\right){\displaystyle \frac{A}{c_{\mathrm{s}}}},\end{array}$$

where g₀ is the Medlyn intercept and defaults to 100 mol m⁻² s⁻¹, and g₁ is the Medlyn slope, a PFT-specific parameter. D is the vapor pressure deficit indicating atmospheric water demand, and c_s is the CO₂ partial pressure at the leaf surface relative to the total atmospheric pressure. A is the carbon assimilated through photosynthesis, as follows:

$$\begin{array}{}\text{(2)}& A={\displaystyle \frac{c_{\mathrm{s}}-c_{\mathrm{i}}}{\mathrm{1.6}{r}_{\mathrm{s}}}}.\end{array}$$

The calculation of A is adapted from Bonan et al. (2011). It is based on the Farquhar model (Farquhar et al., 1980) and limited by the photosynthetic capacity given by the LUNA model (Ali et al., 2016). It requires knowledge of the gradient of the CO₂ concentration from the outside to the inside of the leaf and neglects the CO₂ storage at the leaf surface. c_s and c_i are the leaf surface and internal partial CO₂ pressures, and r_s is the stomatal resistance, which is the inverse of g_s. Furthermore, c_s and c_i are calculated.

$$\begin{array}{}\text{(3)}& {\displaystyle }{c}_{\mathrm{s}}=c_{\mathrm{a}}-\mathrm{1.4}{r}_{\mathrm{b}}A\text{(4)}& {\displaystyle }{c}_{\mathrm{i}}=c_{\mathrm{a}}-\left(\mathrm{1.4}{r}_{\mathrm{b}}+\mathrm{1.6}{r}_{\mathrm{s}}\right)A\end{array}$$

The factor of 1.4 refers to the diffusivity ratio between CO₂ and H₂O gases in the leaf boundary, and 1.6 is the same ratio in the stomata. The equations for A, g_s, c_i, and c_s are computed iteratively until c_i converges, using a hybrid algorithm with the secant and Brent methods (Lawrence et al., 2018). The photosynthesis is scaled to the canopy GPP by considering the effect of sunlit-to-shaded-area ratios of the total leaf area.

The water input from the atmosphere to the land surface can be snow accumulating on the ground, streamflow, lake water, moisture intercepted by the vegetation canopy, or water infiltrating the ground. The water in the ground percolates through 20 soil layers and is stored, directly evaporated, or taken up by plant roots relative to their transpiration demand. Hydraulic stress in a plant is calculated in a hydraulic framework using Darcy's law for transient porous media flow (Bonan et al., 2014).

The transpiration flux T is calculated with the resulting r_s from above.

$$\begin{array}{}\text{(5)}& T={\displaystyle \frac{e_{\mathrm{s}}-e_{\mathrm{i}}}{r_{\mathrm{s}}}}\end{array}$$

e_s is the H₂O vapor pressure at the leaf surface, and e_i is the saturation H₂O vapor pressure resulting from the leaf temperature. If T cannot meet the atmospheric water demand because of a soil moisture shortage, CLM5–BGC introduces water stress and attenuates g_s based on that transpiration deficit factor. Through decreased g_s, water stress also regulates photosynthesis, A.

Total evapotranspiration is then determined by summing the transpiration and evaporation from vegetation interception, surface water, the ground, and, potentially, snow.

2.1.2 Setup of the European CLM5

The European Coordinated Regional Climate Downscaling Experiment (CORDEX; Giorgi et al., 2009) domain delimited the extent of this study, matching with the extent of regional atmospheric models. With a resolution of 3 km (0.0275°), our grid contains 1544×1592 grid cells, including the ocean. We used standalone CLM5 with the activated BGC module and stub models for ice, sea, and waves.

The simulations were forced by the Consortium for Small-Scale Modeling (COSMO) Reanalysis 6 (Bollmeyer et al., 2015; Wahl et al., 2017), a 6 km resolution data set providing meteorological variables over the European CORDEX domain from 1995 to 2019. The main advantage of using this reanalysis is the high resolution and a better representation of seasonal precipitation intensities compared to a coarser resolved global reanalysis (Bollmeyer et al., 2015). Using this forcing in high-resolution LSM simulations should lead to a more accurate simulation of sub-surface and surface hydrological fluxes, especially in regions with a relatively heterogeneous land surface (Wahl et al., 2017; Prein et al., 2016).

The static surface information was initialized for the year 2000 and was determined using input data from a standard repository (Lawrence et al., 2018). These data include land use information (Hurtt et al., 2020), PFT distribution maps (Lawrence and Chase, 2007), soil texture (IGBP, 2000), and slope and elevation information (Earth Resources Observation and Science (EROS) Center, 2017).

The CLM5–BGC needs initial conditions for the carbon pools. For that reason, a spin-up workflow is necessary to bring the carbon pools and fluxes of carbon to a steady state before starting with production simulations. The spin-up method consists of two steps. First, an accelerated decomposition simulation step, where carbon pools are artificially minimized. Second, a conventional simulation step, growing the carbon pools to the desired equilibrium state. During both spin-up steps, the atmospheric forcing from 1995 to 2012 was cycled (i.e., a cycling period of 18 years). The progress towards a steady state is monitored by assessing the difference in total carbon fixed in the ecosystem between a selected year within the last 18-year cycling period and the same year in the previous cycling period. C_tot,y is the total ecosystem carbon (including vegetation and soil) in the year y, and $C_{\mathrm{tot},y-t}$ is the complete ecosystem carbon in the year y−t. A grid cell's carbon pools are in carbon equilibrium if the following is fulfilled:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{\Delta }{C}_{\mathrm{tot}}}{t}}<\mathrm{1}\mathrm{g}\mathrm{C}{\mathrm{m}}^{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}.\end{array}$$

The following conditions define the final steady state on the continental scale.

1. A total of 97 % of the grid cells (and the total area) is in equilibrium.

2. The change in continental ecosystem carbon across the continent is lower than 2 Tg C yr⁻¹ for the three preceding cycle periods.

The soil organic matter carbon pools in high northern latitudes were the slowest to reach equilibrium, which was reached after just about 1500 simulation years.

After the spin-up, we conducted a 24-year (1995 until 2018) transient simulation starting with the initial conditions established by the spin-up. We output the simulated variables from two model levels for the analyses.

1. CLM5_PFT. This is the model's native resolution of vegetation-related states and fluxes calculation. Using output at this level (not the default configuration) allows for multiple time series per grid cell, with each corresponding to a single PFT. This enables a selection of modeled data as needed. For instance, when comparing model data to ecosystem level measurements, CLM5_PFT relates to the simulated time series of the corresponding PFT, resulting in an adequate assessment of model functions. When comparing to in situ observations, we will refer to CLM5_PFT when we subset the ICOS site location and the agreeing PFT from the CLM5 data.

2. CLM5_grid. The grid cell level output aggregates the PFT and the other tiles (i.e., croplands, urban areas, and lakes) that compose the grid cell area. Consequently, these data do not relate to a single functional type. Instead, it informs us about the average state and fluxes in the grid cell area. In this study, CLM5_grid designates CLM5 data extracted from the grid cell closest to the station's location.

2.2 Evaluation data

2.2.1 Station data

As ground truth data in the comparisons, we used the ICOS research infrastructure, which has a station observation network spanning 14 European countries (ICOS RI, 2021). Each station has at least one eddy covariance measurement tower and incorporates a processing workflow following a standardized protocol. We use the curated data, the Warm Winter 2020 data set (Warm Winter 2020 Team and ICOS Ecosystem Thematic Centre, 2022), which consists of homogenized variable time series following the ONEFlux data pipeline (Pastorello et al., 2020). The ICOS Warm Winter 2020 data have measurements of 73 stations totaling over 800 station years (available years are station-dependent) corresponding to multiple land cover types (see Fig. 1 for a map with the station locations and Table S1 for more information on the available years per station). Note that the land cover type indicated by the ICOS site metadata and represented in the measurements refers to the predominant PFT in the footprint of the eddy covariance station. We omitted the stations over wetland and mixed forest land cover types to ensure a coherent analysis because no PFT counterpart is implemented in CLM5_PFT. Also, shrub PFTs were not included in our analyses because there were insufficient shrubland sites in the ICOS data to support a robust evaluation. The analyses also excluded stations whose land cover type was not included in metadata sites (e.g., DEIMS-SDR at https://deims.org, last access: 20 February 2024), leaving a total of 42 stations for our analyses. Because the land cover types from the selected sites correspond well with PFTs in CLM5, we will also refer to them as PFTs.

The processing workflow of the Warm Winter 2020 data extracts daily time series for GPP, partitioned from the net ecosystem exchange (NEE), using the nighttime method and a dependence on a variable friction velocity threshold (in g C d⁻¹; GPP_NT_VUT_REF). We retained negative GPP values in these data, which stem from the uncertainty in the NEE measurements and partitioning method, to avoid introducing bias into the GPP distributions (Reichstein et al., 2012; Pastorello et al., 2020). For the ET evaluation, we also extracted the gap-filled latent heat flux (W m⁻²; LE_F_MDS). Importantly, we verified our results by checking for inconsistencies in the analysis of ICOS NEE (NEE_VUT_REF), ecosystem respiration (RECO_NT_VUT_REF), and energy-balance-corrected latent heat flux (LE_CORR).

The conversion of latent heat (W m⁻²) into ET (mm d⁻¹) is achieved by multiplying with the factor 0.035, assuming a constant enthalpy of vaporization decoupled from temperature because variable enthalpy has a negligible effect on the overall outcome of the conversion.

Last, we use the leaf area index (LAI) from the ICOS archive final quality data set (ETC L2 Archive). LAI is measured but only sparsely available, starting from 2017, and thus only has 2 years intersecting with our study period (2017 and 2018). Furthermore, the data within this intersection period are only available for a smaller number (in relation to the eddy covariance (EC) data above) at evergreen needleleaf forest and cropland sites. Therefore, we do not include the analysis in the main text but include these results only in the Supplement for the context of the main analyses of ET and GPP.

2.2.2 Remote-sensing and reanalysis data

To assess the CLM5 performance in the context of additional complementary data products, we include remotely sensed GPP data from the Global Land Surface Satellite (GLASS; Liang et al., 2021). The GLASS GPP product uses the Moderate Resolution Imaging Spectroradiometer (MODIS) and Advanced Very High Resolution Radiometer (AVHRR) sensors and the revised light use efficiency (LUE) model (Zheng et al., 2020) in eight-times-daily resolution in time and 0.05° resolution in space. We also compare the CLM5 outputs with GLASS ET data, which apply a multi-model ensemble (e.g., MODIS-ET and remote-sensing Penman–Monteith ET) to remote-sensing information to estimate eight-times-daily latent heat on a 0.05° grid. We convert latent heat to ET, as described in Sect. 2.2.1. Similarly, MODIS-derived GLASS LAI data (Ma and Liang, 2022) are used in this study to provide context to the ET and GPP analyses (same 0.05° grid and eight-times-daily resolution).

Last, we use LAI and ET reanalysis data for evaluation, which fuse observations and models. They are the European Center for Medium-Range Weather Forecasts Reanalysis 5 – Land product (ERA5L; Muñoz Sabater, 2019), which has a spatial resolution of 0.1° and hourly temporal resolution, and the Global Land Evaporation Amsterdam Model (only ET; GLEAM version 3.5a; Martens et al., 2017), which has a spatial resolution of 0.25° and daily temporal resolution.

2.3 Data processing

First, the remote-sensing and reanalysis data are bilinearly remapped to the 3 km European CORDEX grid and interpolated to eight-times-daily means for 1995–2018. The ICOS observation time series are interpolated to eight-times-daily means for each station whose data availability overlaps with our study period. Then, we extracted the CLM5_grid, GLASS, ERA5L, and GLEAM data from the grid cell closest to the location of each selected ICOS station. Furthermore, we select the time series in CLM5_PFT that coincides with that grid cell and the station's predominant PFT. Importantly, we focus only on the four predominant PFTs represented in the entire ICOS station network, namely evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), and croplands (CRO), as outlined in Table 1. Finally, the periods for which station data are absent or of bad quality (determined by the corresponding measurement or gap-filling quality flag in the ICOS data) are discarded from the simulations to ensure we are comparing the same set of conditions.

Table 1 The predominant plant functional types (PFTs) in the Integrated Carbon Observation System (ICOS) Warm Winter 2020 observation data set that correspond with the International Geosphere–Biosphere Programme (IGBP) land cover classifications, the number of corresponding sites, and the accordant PFTs in the European Community Land Model v5 (CLM5) setup.

The ICOS observations were also interpolated to eight-times-daily means, encompassing a timescale with significant variability in ecosystem processes (De Pue et al., 2023), to match the coarsest time resolution of other data sets (i.e., GLASS remote sensing) and thus to facilitate comparison of processes at the same scale. For a consistent comparison, the analyses only account for time steps for which valid values are present for all data sources. We evaluate the data for each variable over each station and groups of stations with the same PFT.

2.4 Analyses

2.4.1 Yearly evolution and statistical distributions

We calculate ET and GPP PFT-specific phenology (mean sub-annual dynamics), resulting in day-of-year (DOY) plots. This is done by averaging the same eight-times-daily time step across years for each site and calculating the mean and standard deviation of site-specific DOY belonging to one PFT.

Furthermore, we determined the statistical distributions as probability density functions resulting from the Gaussian kernel density estimate (Scott, 1992). Subsequently, the distribution moments (mean, variance, skewness, and excess kurtosis) are calculated. The distributions and their moments are based on all eight-times-daily values corresponding to one PFT for each data source. The uncertainties in the distribution moments are calculated based on Harding et al. (2014).

2.4.2 Shift in phenological events

The three analyzed phenological events of ET and GPP – the start of the growing season, the peak, and the end of the growing season – are determined for each PFT group and data source as the average DOY of the event among the stations and available years within that PFT group for each variable. The eight-times-daily time series of each variable was first smoothed with a 1-dimensional Gaussian filter to rule out potential errors due to small-scale variability and dampen the effect of potential outliers. More specifically, the peak timing is the mean DOY of the overall maxima of the smoothed averaged yearly evolution across stations for each PFT and data source. The start and the end of the growing season were determined by the mean DOY of the two infliction points (Li et al., 2023; Lian et al., 2020; Whitcraft et al., 2015) of the smoothed yearly averaged evolution across stations for each PFT and data source. The shift in these events is simply the difference in the determined mean PFT-specific DOY between the models and the observations. As a measure of uncertainty in the mean PFT-specific DOY, we also calculate the standard deviation of the DOY of the events across stations in each PFT group.

2.4.3 Performance metrics

The percent bias (PBIAS) measures systematic model error and is calculated as follows:

$$\begin{array}{}\text{(7)}& \text{PBIAS}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{X}_{\mathrm{S},i}-X_{\mathrm{O},i}}{{\sum }_{i=\mathrm{1}}^n{X}_{\mathrm{O},i}}}\times \mathrm{100},\end{array}$$

where n is the number of time steps, X_S,i is the simulated value of the variable X at the time i, and X_O,i is the observed value of the variable X at the time i. If the PBIAS for variable X is positive then the model overestimates; if it is negative, it underestimates the observed variable X. In our analysis, X_i is the interpolated eight-times-daily mean.

Furthermore, we estimated the root mean square error (RMSE) to indicate model accuracy and the root mean square difference (RMSD) to indicate similarity. RMSE and RMSD are calculated the same. However, the term “error” assumes the truthfulness of the reference data. Hence, we use the RMSD when comparing data only between models.

$$\begin{array}{}\text{(8)}& \text{RMSE}=\text{RMSD}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(X_{\mathrm{S},i}-X_{\mathrm{O},i}\right)}^{\mathrm{2}}}{n}}}\end{array}$$

A RMSE close to zero indicates that the model approximates the observations nicely. Similarly, a low RMSD reveals a high similarity between the two analyzed series. We calculate these metrics on a per-station basis and a set of stations belonging to the same PFT.

2.4.4 Modified Taylor diagrams

The Taylor diagram (Taylor, 2001) depicts multiple model performance indices in a single diagram by making use of the relationship of the calculation terms of the standard deviation, correlation, and RMSE. Their relationship can be summarized in the following equation of error propagation:

$$\begin{array}{}\text{(9)}& {\text{RMSE}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{O}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{S}}^{\mathrm{2}}-\mathrm{2}{\mathit{\sigma }}_{\mathrm{O}}{\mathit{\sigma }}_{\mathrm{S}}r,\end{array}$$

where σ_O is the standard deviation of the observations, σ_S is the standard deviation of the simulation, and r is the Pearson correlation coefficient. The multi-variate diagram can be constructed due to the geometric relationship between these statistical indices through the law of cosines. Thereby, plotting the calculated Pearson correlation against the standard deviation of the models and the observation on a trigonometric polar plane, the RMSE manifests as the polar Euclidean distance from the reference observations. We calculate the standard deviation and the Pearson correlation on the PFT-grouped stacked time series and plot these for e for each data source on one Taylor diagram per PFT. We modify the default Taylor diagram by scaling each marker's size by the absolute PBIAS for the corresponding source and PFT.

Figure 1 The share of represented plant functional types (by color, with evergreen needleleaf forest (ENF; green), deciduous broadleaf forest (DBF; orange), grasslands (GRA; purple), and croplands (CRO; pink)) in panel (a) in the ICOS station network used in subsequent analyses and panel (b) in the corresponding grid cells in our European CLM5 setup. In panel (c) is a map showing the locations of the ICOS stations, with the marker type indicating their PFT and the color of the marker indicating their hydro-climate (adapted from Jafari et al., 2018) and based on the mean annual precipitation from the Consortium for Small-Scale Modeling (COSMO) – Reanalysis 6. Our 3 km European CLM5 simulation domain corresponds to the entire map box in panel (c).

3 Results

3.1 Land surface representation

Before evaluating the GPP and ET variables from CLM5 and how they are compared with observations, we first assess if the PFT composition of the entire ICOS station network is comparable to the PFT composition in the respective cells selected in CLM5_grid. This is important as GPP and ET magnitudes, variability, seasonality, drought responses, and trends strongly depend on the present vegetation type. In Fig. 1, we observe that ENF, the PFT of almost half of the present ICOS stations, represents only around a quarter of the corresponding CLM5_grid area. DBF also covers a smaller share of the area in those grid cells than in the ICOS station network. On the other hand, GRA and CRO are overrepresented in CLM5_grid compared to the share of respective ICOS stations. Consequently, when comparing with the ICOS observations, the selected data from CLM5_grid data are, on average, over-representing the functionality of GRA and CRO and underrepresenting ENF and DBF, which hampers the evaluation of CLM5_grid with in situ ET and GPP. Hence, we also included the respective CLM5_PFT GPP and ET in the subsequent analysis, enabling an accurate assessment of the functionality and relationships between PFT in the model. Additionally, we assess the similarities and differences between the two model scales, CLM_grid and CLM_PFT, and their approximation to the observations.

3.2 General model performance

This section presents model performance indices correlation, RMSE, and PBIAS, comparing each model's ET and GPP estimates with measurements from the ICOS sites. We compared the RMSE and PBIAS on a per-site basis (Tables S2 and S3), which yielded good results for most sites. The focus of this study, though, is the performance of PFT aggregations through combining data from sites that belong to the same PFT.

Figure 2 Modified Taylor diagrams with observations from the Integrated Carbon Observation System (ICOS) of evapotranspiration as reference (black markers), showing model performances between the years 1996–2018 (years varying by station; see Table S1. Data sources by color include the Community Land Model v5 (CLM5), CLM5_grid in red; CLM5_PFT in yellow; Global Land Surface Satellite (GLASS) in green; European Center for Medium-Range Weather Forecasts Reanalysis 5 – Land (ERA5L) in brown; and Global Land Evaporation Amsterdam Model (GLEAM) in purple). Each diagram shows these plots for one plant functional type. (a) Evergreen needleleaf forest (ENF; circles). (b) Deciduous broadleaf forest (DBF; triangles). (c) Grasslands (GRA; squares). (d) Croplands (CRO; crosses). The azimuth angle indicates the Pearson correlation with the ICOS data, the radial distance is the standard deviation, and the semicircles centered at the reference standard deviation show the root mean square error (RMSE). The size of each marker indicates the percent bias (PBIAS).

Figure 2 shows modified Taylor plots visualizing the performance indices of the model ET against observations for each PFT. For more specific information, Table S4 lists the number of ET eight-times-daily time steps that went into calculating these indices and their values. For ENF, all of the models indicate a correlation of around 0.8 with the ICOS observations, and CLM5_grid, CLM5_PFT, and GLEAM have a similar variability to ICOS. CLM5_PFT has a higher absolute RMSE and a smaller absolute PBIAS than CLM5_grid for ET across PFTs, except in CRO. Notably, the systematic bias in CLM5 is generally negative, with the same exception. On the other hand, ERA5L, GLASS, and GLEAM exhibit a generally positive systematic bias for ET. ERA5L and GLASS show more significant deviations from the ICOS ET observations at ENF and DBF than CLM5_PFT and CLM_grid but have smaller RMSE values at GRA and CRO. GLEAM has generally low RMSEs and performs best among the models simulating ET at ENF and CRO. The most considerable systematic ET biases are found for ERA5L at CRO and DBF sites, followed by GLASS for the same PFTs. The low absolute PBIAS of CLM5_PFT across all PFTs and the lower correlation than the other model data at GRA and CRO points to potentially missing or simplistic representations of ecohydrological processes or management. Besides, all models approximate the ICOS ET observations fairly well, with correlations mostly over 0.8 but with partly high systematic biases by ERA5L at DBF and CRO sites.

Figure 3 Modified Taylor diagrams with observations from the Integrated Carbon Observation System (ICOS) of the gross primary production as a reference (black markers) and showing model performances between the years 1996–2018 (years varying by station; see Table S1. For colors, labels, and acronyms, refer to Fig. 2).

In Fig. 3, we show modified Taylor diagrams with the GPP performance indices of the models against the ICOS observations for each of the selected PFTs. For more specific information, Table S5 lists the number of GPP eight-times-daily time steps that went into calculating these indices and their values. CLM5_PFT performed better than CLM5_grid in approximating the ICOS GPP observations at DBF sites, showing a higher correlation and lower RMSE and GRA sites. Conversely, CLM5_grid is closer to the observations for ENF and CRO PFTs. The GLASS data show the lowest GPP RMSEs and highest correlation values concerning ICOS measurements across all PFTs. All models approximated the ICOS GPP best (lowest RMSE) at ENF, and the worst performance was at CRO sites. Furthermore, all models exhibit a negative, systematic bias in simulating the observed GPP across all PFTs. Especially at DBF and GRA PFTs, CLM5_grid, CLM5_PFT, and GLASS show large systematic underestimations of the measurements. CLM5_PFT has a notably small PBIAS related to the ICOS data for ENF and CRO sites. Especially at CRO sites, all models showcase comparatively low correlation values (<0.7). While the correlation is high (>0.75) for all models at DBF and GRA sites, especially for CLM5_PFT and GLASS at DBF sites (0.93 and 0.92), the high PBIAS hints that modeled data do not incorporate important processes or management practices that cause to the high carbon uptake at DBF sites over the long term. Because of the slowly evolving carbon states in the terrestrial ecosystems, the initial conditions of the carbon pools (e.g., soil organic matter and carbon in plant organs in the vegetation) could be a cause for the difference in the magnitude of the GPP.

Figure 4 In the left column are the yearly evapotranspiration (ET) evolutions averaged across stations belonging to one plant functional type (rows show evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), and croplands (CRO)) and across the years (available years vary per station; see Table S1). We differentiate the data source by color (Integrated Carbon Observation System (ICOS) observations in black; Community Land Model v5 (CLM5), CLM5_grid in red; CLM5_PFT in yellow; Global Land Surface Satellite (GLASS) in green; European Center for Medium-Range Weather Forecasting Reanalysis 5 – Land (ERA5L) in brown; and Global Land Evaporation Amsterdam Model (GLEAM) in purple). The corresponding standard deviations across the sites and across the years are plotted in the right column to measure the spread around this mean.

Figure 5 Mean shifts in ET phenological events (the start of the growing season, peak, and the end of the growing season) between the Integrated Carbon Observation System (ICOS) observations (solid black line) and the models (by color, namely the Community Land Model v5 (CLM5), CLM5_grid in red; CLM5_PFT in yellow; Global Land Surface Satellite (GLASS) in green; European Center for Medium-Range Weather Forecasts Reanalysis 5 Land (ERA5L) in brown; and Global Land Evaporation Amsterdam Model (GLEAM) in purple) among sites belonging to one plant functional type, namely evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), and croplands (CRO). On the x axis is the day of the year of the event. Error bars in grey correspond to the standard deviation of the day of the event in the models across sites of one plant functional type, and the error bars in black correspond to the standard deviation across the respective observations.

3.3 PFT phenology and its variability

3.3.1 ET

This section describes the results of the investigation on the mean and the standard deviation of the yearly evolution of ET across PFTs and data sources (Fig. 4a, c, e, g). We will analyze the ET mean and standard deviation for each PFT sequentially. On average, the annual evolution of ET for CLM5_grid and CLM5_PFT compares well to the ICOS measurements, as already hinted by the good correlation values in the previous section. They also capture the observed seasonal transitions between low winter ET and high summer ET well. Except for CRO sites, CLM5_grid and CLM5_PFT ET are slightly lower than the ICOS observations throughout the year but especially in summer (mean PBIAS of −13.08 % and −18.70 %, respectively; see Table S4). ERA5L and GLASS overestimate ET at sites of all PFTs, most predominantly in the ENF and DBF sites and during summer (mean PBIAS of +28.64 % and +18.25 %, respectively). The magnitude of variation across sites within each PFT (Fig. 4b, d, f, h) is captured well, generally showing smaller variations at DBF and CRO sites and larger variations at ENF and GRA. Some specific aspects of this variation across sites are captured best by CLM5_PFT. The bimodality of the intra-station variation at GRA sites is seen across the year (Fig. 4f), and the peak variability across stations at CRO sites is seen in the second half of the year (Fig. 4h). This exhibits the ability of CLM5_PFT to differentiate ET between stations and the PFTs better than CLM5_grid and the other models. The GLASS ET variability across stations compares remarkably well to the observed across ENF at DBF sites (Fig. 4b and d).

Figure 5 reveals the shift in the timings of key phenological events based on ET (growing season start, summer peak, and growing season end) between each model and the ICOS observations. Generally, for ENF and DBF sites (Fig. 5a, b), all models show the earlier occurrence, and at CRO (Fig. 5d), they show a later occurrence of these phenological events than the measurements. CLM5_PFT has the mean timing of the events within the standard deviation of the ICOS timing across all PFTs. However, it shows a substantial variability, larger than observed in the event timings across GRA sites. Similarly, GLASS and CLM5_grid show close approximations to the observed timings but simulate all these events significantly earlier at DBF sites and significantly later at CRO sites, with little variation in the timings across sites. The ERA5L and GLEAM data exhibit a much earlier growing season start (24 and 20 d earlier) and summer peak (16 and 12 d earlier, respectively) than observed by ICOS at GRA sites.

Figure 6 In the left column are the yearly gross primary production (GPP) evolutions averaged across stations belonging to one plant functional type (rows show the evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), and croplands (CRO)) and across the years (available years vary per station; see Table S1). We differentiate the data source by color (Integrated Carbon Observation System (ICOS) observations in black; Community Land Model v5 (CLM5), CLM5_grid in red; CLM5_PFT in yellow; and Global Land Surface Satellite (GLASS) in green). The corresponding standard deviations across the sites and across the years are plotted in the right column to measure the spread around this mean.

Figure 7 Mean shifts in GPP phenological events (the start of the growing season, peak, and the end of the growing season) between the Integrated Carbon Observation System (ICOS) observations (solid black line) and the models (by color, namely the Community Land Model v5 (CLM5), CLM5_grid in red; CLM5_PFT in yellow; and Global Land Surface Satellite (GLASS) in green) among sites belonging to one plant functional type, namely evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), and croplands (CRO). On the x axis is the day of the year of the event. Error bars in grey correspond to the standard deviation of the day of the event in the models across sites of one plant functional type, and the error bars in black correspond to the standard deviation of the respective observations.

3.3.2 GPP

The GPP values of all PFTs show a summer peak and a low period in winter (Fig. 6). The negative values present in the ICOS measurements are caused by the processing of the measurements by ICOS and are, therefore, not represented by CLM5 or GLASS. Again, a general underestimation of observed GPP is shown across all PFTs (Fig. 6a, c, e, g), particularly during the summer months from all models. CLM5_PFT shows a larger GPP than CLM5_grid and, therefore, has a lower systematic bias (mean PBIAS across PFTs of −19.61 % and −27.65 %; see Table S5). GLASS GPP is closer to the ICOS GPP at ENF, DBF, and GRA and has the lowest mean PBIAS across PFTs of −16.67. The most substantial underestimations are at DBF during summer (Fig. 6c), where CLM5_grid and CLM5_PFT have a PBIAS of −38.88 % and −43.76 %, and GLASS −24.52 %. The GPP variability across sites is, similar to ET, the lowest at DBF sites. Notably, GLASS remote-sensing GPP underestimates the variability among sites of one PFT substantially throughout the year at GRA and CRO sites (Fig. 6f, h). The observed variability dynamics across the year, e.g., the bimodality at GRA sites (Fig. 6f) that was also visible for ET, is captured best by CLM5_PFT. However, not all models capture the behavior of CRO GPP inter-site variability (Fig. 6h). This supports the suspicion of the influence of management and missing processes in CRO in the models, possibly concerning the timings of planting, fertilizing, and harvesting the crops as the cause of these mismatches. The overall negative systematic bias in the models points at potentially missing sensitivities to or lower levels of, e.g., atmospheric CO₂ and vapour pressure deficit (VPD) that have been recently found to increase the water use efficiency and carbon assimilation (Poppe Terán et al., 2023; Friedlingstein et al., 2023).

Shifts in phenological events between the observations and the models are already noticeable in Fig. 6 but are quantified and visualized in detail in Fig. 7. CLM5_PFT and CLM5_grid predominantly simulate the timing of these events within the standard deviation across ICOS stations for each PFT. In the GLASS GPP data, the events are shifted further from the measurements, most notably at DBF sites (16 d earlier growing season start and 11 d earlier summer peak) and at CRO sites (22 d belated peak and 45 d belated end of the growing season). Generally, in both CLM5 scales, the shifts to the ICOS observations were the largest in CRO. Similar to the ET event timings, CLM5_PFT shows the largest variation in these timings among the models, especially at GRA sites, and also considerable differences in the timing of the growing season end of ENF sites. These findings confirm the ability of CLM5_PFT to approximate PFT-specific variation in ecosystem processes, but the contrasting results of the model performance indices will be further reviewed in Sect. 4.

Figure 8 The probability density curves for all evapotranspiration (ET) values from stations belonging to the selected plant functional types, namely evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), and croplands (CRO). The data source differs by color (Integrated Carbon Observation System (ICOS) observations in black; Community Land Model v5 (CLM5), CLM5_grid in red; CLM5_PFT in yellow; Global Land Surface Satellite (GLASS) in green; European Center for Medium-Range Weather Forecasts Reanalysis 5 Land (ERA5L) in brown; and Global Land Evaporation Amsterdam Model (GLEAM) in purple).

3.4 Statistical distributions

3.4.1 ET

In this section, we describe the results of the statistical distributions of ET in the model and the observations for each PFT. Then, we give more details on the moments of these distributions and how the models compare to the observations. Generally, the models approximate well the shape of the distributions (Fig. 8), with a pronounced peak in the occurrence of positive ET values close to 0 that represent low winter values across all PFTs and, moreover, the slowly decreasing frequency of values towards the high ET summer values, which is more variable among the models. The variability in the summer peak magnitude (see Sect. 3.3) among stations of the same PFT causes the ICOS and CLM5 ET distributions to have only a slightly pronounced second mode at the high summer ET values. On the other hand, the ERA5L and GLEAM ET distributions show a very pronounced second mode at the higher ET values for each PFT. This hints at the lower variability in the summer peak magnitude among these stations, which misrepresents the observed high variation in ICOS.

The moments of these distributions give more insights into their specific characteristics. Furthermore, differences in moments between the observations and the models can yield important information on potential misrepresentations (Fig. 9). For example, a differing mean between ICOS and a model points to a general shift in the distribution, specifically its center of mass. Therefore, we confirm a shift in ET distributions of ERA5L, GLASS, and GLEAM towards higher values for all PFTs in reference to ICOS. CLM5_grid and CLM5_PFT have lower means, except for CLM5_PFT at CRO. The second moment, the variance, informs us about the variability in values. Notably, GLEAM data underestimate, and GLASS data overestimate, the observed variability in ET at all PFTs. CLM5_PFT has a broad range of variability across PFTs, which corresponds well with ICOS observations, while CLM5_grid and the other models show a very similar level of variability independent of the PFT. All models agree with the observed positive sign of the skewness (indicating a longer right tail of the distribution) for all PFTs. And while all the models simulate a platykurtic (negative excess kurtosis and pronounced relative tail) characteristic of the distributions across PFTs, ICOS shows a leptokurtic (positive excess kurtosis, less pronounced tails, and more pronounced peak) behavior at ENF and CRO sites. Furthermore, the variation in the reach model's skewnesses and kurtoses (y axis ranges for each color in Fig. 9c and d) across the PFTs is considerably lower than the observed ranges (corresponding x axis ranges). Altogether, these findings showcase the ability of CLM5_PFT to model intra- and inter-PFT ET variance better than the other considered models, on the one hand, but also showcase the shortcomings of all the considered models in representing the variation in the extreme ends of the ET distributions across all PFTs on the other hand.

Figure 9 The mean (a), variance (b), skewness (c), and excess kurtosis (d) of the evapotranspiration (ET) distributions (visualized in Fig. 8) from the models (y axis; colors, namely Community Land Model v5 (CLM5), CLM5 grid in red; CLM5PFT in yellow; Global Land Surface Satellite (GLASS) in green; European Center for Medium-Range Weather Forecasts Reanalysis 5 Land (ERA5L) in brown; and Global Land Evaporation Amsterdam Model (GLEAM) in purple) as opposed to the corresponding values from observations (x axis) aggregated for each plant functional type (marker type), namely evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), and croplands (CRO). The error bars are the standard errors in the respective moment, depending on the sample size.

Figure 10 The probability density curves for all gross primary production (GPP) values from stations belonging to the selected plant functional types are shown, including evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), and croplands (CRO). The data source differs by color (Integrated Carbon Observation System (ICOS) observations in black; Community Land Model v5 (CLM5), CLM5_grid in red; CLM5_PFT in yellow; and Global Land Surface Satellite (GLASS) in green).

3.4.2 GPP

We continue to delineate the results of the same analyses for the GPP distributions and their moments (Fig. 10). The frequency peaks at the low GPP values, which correspond to the base winter GPP, are overestimated by all models at ENF, DBF, and GRA. This could partly be explained by negative GPP values in the ICOS data, which the models do not represent. By definition, there is no negative GPP. However, these negative values are given through the uncertainty range of the NEE partitioning method and are retained in the analysis to preserve the partitioning distribution (Reichstein et al., 2012; Pastorello et al., 2020). This is probably related to underestimating the observed winter GPP in ENF and GRA sites seen in Fig. 6a and e. Another striking finding is the missing occurrence of the highest observed GPP values in the models at all PFTs, but most noticeable at DBF sites, where the upper-half range of GPP values (>12 g C d⁻¹) is not represented in any model. The overrepresented mid-range GPP values and the partly pronounced second modes in the mid-range GPP values across PFTs are possibly caused by the low summer peaks and low variability across sites (see Fig. 6).

The models show lower GPP means than the ICOS measurements for all PFTs in Fig. 11a. Similarly, for all models across all PFTs, the underestimated GPP variance indicates a lower spread of the PFT distributions than in ICOS. While models agree on the positive skewness of the GPP distribution (skewed to the left), the largest skewness at CRO sites is not well represented by all the models. Finally, similar to the findings with ET kurtosis, the models fail to distinguish the distinct leptokurtic characteristics (tails that are less heavy) of the GPP distribution of CRO sites compared to the other PFTs, as seen in the observations. Across PFTs and for all models, the ranges spanned by the intra-PFT distribution moments are smaller than those observed. Most strikingly, the GPP variance range across PFTs, which is, among the models, the largest for CLM5_PFT (between 8 and 12 g C d⁻¹), is much smaller than for ICOS (11 to 27 g C d⁻¹). This suggests the models do not simulate GPP differently enough between the PFT groupings. Thus, model development and parameter optimization studies that aim to improve these representations should focus on enhancing the variability in DBF.

Figure 11 The mean (a), variance (b), skewness (c), and excess kurtosis (d) of the gross primary production (GPP) distributions (visualized in Fig. 10) from the models (y axis; colors, namely Community Land Model v5 (CLM5), CLM5 grid in red; CLM5PFT in yellow; and Global Land Surface Satellite (GLASS) in green), as opposed to the corresponding values from observations (x axis) aggregated for each plant functional type (marker type), namely evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), and croplands (CRO). The error bars are the standard errors in the respective moment, depending on the sample size.

3.5 The inter-site similarity of PFT groups

To support the interpretations of our findings, we quantify the similarity of ET and GPP across sites of the same PFT and compare the differences between the models and the observations. In this section, we analyze the mean RMSD of each PFT per ET and GPP data sources. A low RMSD indicates that the stations corresponding to one PFT are similar, while a high RMSD hints at a greater diversity within the PFT. By comparing the mean RMSD per PFT for ET and GPP across data sources, we can evaluate how much diversity is captured in the data of a particular PFT in the observations and models. The standard deviation of the RMSD for each PFT gives information on the spread of the inter-site RMSDs within the PFT group around that mean.

Figure 12 The bars indicate the mean of the root mean square difference (RSMD) of evapotranspiration calculated for sites with the same plant functional type. The error bars are their standard deviation. Low values indicate high similarity between the sites, and high values show high dissimilarity. The color of the bars differentiates the data source (Integrated Carbon Observation System (ICOS) in black; Community Land Model v5 (CLM5), CLM5_grid in red; CLM5_PFT in yellow; Global Land Surface Satellite (GLASS) in green; European Center for Medium-Range Weather Forecasts Reanalysis 5 Land (ERA5L) in brown; and Global Land Evaporation Amsterdam Model (GLEAM) in purple).

Figure 12 shows that CLM5_grid and GLEAM have lower ET time series differences between the corresponding sites for all PFTs than ICOS. CLM_PFT has a lower mean RMSD than CLM5_grid among ENF and DBF sites. Both CLM5_PFT and CLM5_grid underestimate the observed diversity of ET at ENF and DBF sites. Interestingly, the variation in ERA5L and GLASS ET time series for ENF is higher than observed, and they also show the most significant variation in RMSD. Meanwhile, the DBF mean RMSD of all models is lower than that of ICOS. CLM5_PFT shows a higher diversity of ET between GRA sites and CRO sites than CLM5_grid. The CLM5_PFT surpasses the observed mean RMSD for the GRA PFT, highlighting the potential to simulate GRA sites variably. All other models underestimate it slightly (CLM5_grid, ERA5L) or in a more pronounced manner (GLASS and GLEAM). Particularly at CRO sites, the ET RMSD of CLM5_PFT is substantially higher than the other models and at a similar level to ICOS observations. In contrast, all other models show significantly lower mean RMSDs there. Generally, a higher ET RMSD mean in a PFT group comes with a higher spread (higher standard deviation) for all data sources. The RSMD in ET between stations is lower for CLM5_grid and GLEAM than for ICOS for all PFTs.

Figure 13 The bars indicate the mean of the root mean square difference (RSMD) of gross primary production calculated for sites with the same plant functional type. The error bars are their standard deviation. Low values indicate high similarity between the sites, and high values show high dissimilarity. The color of the bars differentiates the data source (Integrated Carbon Observation System in black; Community Land Model v5 (CLM5), CLM5grid in red; CLM5PFT in yellow; and Global Land Surface Satellite (GLASS) in green).

Figure 13 shows that for GPP, the models generally have a lower mean RSMD than ICOS across stations for all PFTs, except for CLM5_grid at DBF. CLM5_PFT has a more diversely simulated ET across ENF, GRA, and CRO sites than CLM5_grid. Interestingly, the observed magnitude of the RMSD is lowest for DBF and highest for CRO and has a more extensive range across PFTs than the models. For example, the RMSDs of ICOS data differ by approximately 1.3 g C d⁻¹ between GRA and CRO, while CLM_grid, CLM5_PFT, and GLASS indicate similar RMSDs for those PFTs. Especially CLM5_grid shows a constant within-PFT variability in around 1.5 g C d⁻¹ independent of the PFT. Higher mean GPP RMSD values also come with a higher standard deviation. These results hint at a complex relationship of variability representation within the PFTs. The higher RMSE values of CLM5_PFT in the general model performance analysis (Sect. 3.2) suggest that the variation across sites of one PFT seen here does not directly translate to better model performance. Apart from the magnitude of the variability, its accurate and proportionate timing is pivotal for enhanced model performance.

4 Discussion

Our results show that CLM5_grid and CLM5_PFT approximate the ET observations from ICOS better than GLASS remote-sensing and ERA5L reanalysis but worse than GLEAM reanalysis. Moreover, especially for CLM5_PFT, the systematic error in simulating ET is lower than all other evaluated data sets. For GPP, we found that CLM5_grid and CLM5_PFT performed worse than GLASS data, indicated by a larger PBIAS and larger RMSE. Surprisingly, CLM5_PFT generally had a higher RMSE than CLM_grid but, at the same time, a lower PBIAS. Averaged ET and GPP phenologies were relatively well simulated but exhibited underestimations across all PFTs, especially in DBF, compared to ICOS measurements. CLM5_PFT better captured the PFT-specific mean and standard deviation of the ET and GPP annual dynamics than CLM5_grid and reanalyses and remote-sensing data. The GPP and ET distributions analysis showed underestimations of their observed variability for all models, namely CLM5_grid, CLM5_PFT, GLASS, ERA5L, and GLEAM. Last, we found that for most PFTs, the modeled and remotely sensed data were too similar between stations of the same PFT group compared to the ICOS observations.

4.1 Uncertainty

4.1.1 Observations

Notably, the EC measurements carry uncertainties that might affect the results of this study, especially related to the systematic errors in the simulations. For instance, EC measurements neglect the energy from large eddies. To check for possible inconsistencies, we evaluated the energy-balance-corrected ET (ET_corr) from the ICOS sites (Pastorello et al., 2020). This methodology assumes a constant Bowen ratio to close the energy imbalance. Simulated ET underestimates ET_corr to a greater degree than the non-corrected ICOS ET (Figs. S1, S2), suggesting a higher systematic error than in the analysis of non-corrected ET. Besides that, we discovered the same patterns with the corrected ET, concluding that the energy balance error did not introduce significant bias to our results and the interpretations. Furthermore, GPP is not directly measured but partitioned from NEE. The NEE partitioning method has an underlying uncertainty stemming from potentially unfulfilled assumptions that propagate to the GPP and ER variables in the ICOS data. So, we also ensured that our results remained consistent by evaluating the non-partitioned NEE and the ER variables (Figs. S3, S4, S5, S6). We discovered a substantial underestimation and missing variability in NEE and ER across PFTs in CLM5, confirming the systematic underestimation in our analysis of GPP. While we believe that our analyses have followed meticulous approaches to ensure robust results by applying the ICOS quality flags and comparing these additional variables, many studies still emphasized the biases arising from a shifting footprint with varying wind direction and wind speed and the energy balance correction method assuming a constant Bowen ratio (Jung et al., 2020; Eshonkulov et al., 2019; Chu et al., 2021). Therefore, we encourage developing and using novel and more accurate energy-balance-closure methods (Zhang et al., 2024). Furthermore, dropping bad-quality gap-filled data from the ET and GPP time series might introduce a bias that underrepresents periods of low-friction velocity and atmospheric inversion conditions. Last, based on the geographical distribution of the ICOS station network, the results might misrepresent southern and eastern Europe and semi-arid and arid hydro-climates (Fig. 1; also read Ohnemus et al., 2024). Those factors might have influenced the diversity of ET and GPP values and the ranges of their distributions.

4.1.2 Forcing

Importantly, discrepancies between the COSMO Reanalysis used to force the European CLM5, and the station observations might introduce deviation into our analyses that could hamper interpretations of our results regarding the model functionality. While the high-resolution forcing data already include information from observations through data assimilation, particular locations and conditions might be less well represented than others, and a resulting bias in the meteorological variables would propagate to the simulation of ET and GPP. However, data assimilation approaches minimize the systematic error in the atmospheric model to the observations. Furthermore, the probability and potential influence of including a bias from the forcing of a single location is lowered by considering multiple sites in the performance and statistics of the PFTs. Nevertheless, we assessed the meteorological variables from the COSMO Reanalysis 6 (temperature, shortwave incoming radiation, precipitation, and relative humidity) with the ICOS station data to scrutinize potential errors arising from the forcing. We used the same approach for the GPP and ET evaluation (Figs. S7–S14). We discovered that the forcing variables' average yearly dynamics and distributions represent the ICOS observations well. More minor yet notable misrepresentations include underestimations of shortwave downward radiation and precipitation in summer and relative humidity over GRA and CRO sites throughout the year compared to the measurements. This could explain some of our analyses' ET and GPP underestimations by CLM5. Notably, the mean and variance across the PFTs and their ranking are represented reasonably well for all forcing variables compared to our results with GPP and ET. Furthermore, the skewness and excess kurtosis of the forcing temperature and shortwave downward radiation compare well to the ones from ICOS, indicating well-matching distributions between the COSMO Reanalysis 6 and the observations. However, in particular, the higher-degree moments of the distribution are not well simulated for precipitation and relative humidity. These characteristics of the distributions affect the CLM5 simulations of GPP and ET and might have influenced our results. Further considerations, including ensemble simulations with perturbed forcings, are required to capture the uncertainty introduced into CLM5 fully, but this is beyond the scope of this study.

4.1.3 Static information and initial conditions

The static surface information, including the soil texture, elevation, aspect, land unit, and PFT distributions, affect the simulation of ET and GPP in CLM5. The soil texture composition will define how water is stored and conducted in the soil, contributing to the evaporation from the soil, an essential ET component. Furthermore, the soil texture will influence root water uptake if vegetation is present in the soil column, indirectly impacting plants' transpiration, another critical ET component. Furthermore, ET is regulated by the available energy, which is determined by how the canopy, the elevation, and the aspect of that location influence the incoming radiation. Especially the diversity between these input variables across the locations of the ICOS stations might have played an essential role in the simulation of the PFT-specific ET and GPP distributions.

Last, particularly for CLM5_grid, GLASS, GLEAM, and ERA5L, the distribution of PFTs across the domain and in the grid cells corresponding to the ICOS stations define the equations and parameters that will be used for the calculation of ET and GPP. Consequently, if the grid cells corresponding to ICOS stations are dominated by PFTs that do not comply with the stations' footprints, the simulations of specific PFTs in the model are negatively affected. Importantly, this does not apply to CLM5_PFT because we could select the data that belong to the adequate PFT. Therefore, interpretations of our results relating directly to vegetation functions implemented in CLM5 are here primarily focusing on the CLM5_PFT data.

The initial conditions of the carbon cycle, most notably the size of the soil and vegetation carbon pools, are another source of uncertainty. Essentially, our spin-up and production simulations were restricted to the years where the high-resolution forcing was available (1995–2018). The spin-up simulations, therefore, recycle atmospheric forcings for a substantial period, which we also use in the production simulations. Hence, the production simulations adopted the equilibrium state (incoming carbon equals outgoing carbon) required to conclude the spin-up. However, in natural conditions, there was no carbon equilibrium in the simulated years. Instead, the carbon cycle experiences dynamic changes, such as long-term trends resulting from changing environmental conditions. Many European ecosystems exhibited a net carbon uptake, thus acting as a carbon sink (Pilli et al., 2017; Winkler et al., 2023), and were measured in ICOS accordingly. The negative long-term mean NEE indicates carbon sources, evident across all PFTs in the EC observations (Fig. S4a). On the other hand, the simulations show a NEE close to zero for all PFTs, directly showing the effect of the equilibrium state of the land surface in the model. The results of DBF, which is the most significant carbon sink in the ICOS data and simultaneously shows the largest GPP underestimations by CLM5, underline a potentially important role of the carbon equilibrium in our results. Future work will conduct a more comprehensive spin-up under conditions closer to a real-world carbon equilibrium (the 1950s or earlier) and a transition run before the production simulations to capture the dynamic trends of the land surface processes. Possibly, the bias in the EC measurements towards conditions with low-friction velocity and atmospheric inversion might also cause overestimations of GPP and the resulting carbon sink in ICOS.

4.2 PFT-specific evaluation

While CLM5_PFT showed a smaller systematic error than CLM_grid for most PFT compared to the observations (lower absolute PBIAS), the ability to approximate the observation time series is worse (higher RMSE). A shifting sign in the bias of the CLM5_PFT data explains these counterintuitive results. The presence of both positive and negative bias (in time and across stations) cancels out and yields an overall low PBIAS. In summary, we find in the evaluation that the ET time series of CLM_PFT are not closer to observations than CLM5_grid for any PFT, but CLM5_PFT generally approximates the ET sum over time better than CLM5_grid for ENF, DBF, and CRO. However, it is also clear that, on average, the phenology of CLM5_PFT is closer to the observed than CLM5_grid, for instance, for both ET and GPP at DBF and GRA sites. Furthermore, the timings of the phenological events in CLM5_PFT are most often closer to the observed than in CLM5_grid. Importantly, critical PFT-specific characteristics, like the timing of DBF's steep spring GPP increase, are only captured by CLM5_PFT and the inter-site variability in ET and GPP throughout a standard year. This discrepancy between the evaluation metrics and the vegetation phenology suggests that CLM5_PFT could better capture the PFT-specific variability that ICOS observes. However, this variability is modeled in a way that did not contribute to a low RMSE that is, for instance, shifted in time or space, so the averaged PFT-specific comparisons (the phenology and the distribution moments) compare better with ICOS than CLM5_grid. Further evidence of this explanation is that CLM5_PFT generally captures more variability (higher ET and GPP standard deviation across sites throughout the year for ENF, GRA, and CRO and higher variance for each PFT). This ability to capture more variability than the other models, closer to the observed variability, can improve the represented variability in CLM5_PFT if the suitable variation can be modeled at the right time and location. This spatiotemporal discrepancy of simulated and observed GPP and ET variability could potentially be solved with optimized PFT parameters (Baker et al., 2022; Birch et al., 2021; Cheng et al., 2021; Dagon et al., 2020; Deng et al., 2021; Fisher et al., 2019).

Several past studies also indicated the underestimation of ET and GPP in CLM5 compared to observations (Boas et al., 2023; Strebel et al., 2024; Cheng et al., 2021; Birch et al., 2021), which we confirm in this study. Parameter improvements could also alleviate these general underestimations of GPP and ET across PFTs, especially during summer (Dagon et al., 2020). However, optimal parameters might vary from site to site (Lin et al., 2015), even if they have the same PFT. Thus, CLM5 and, more generally, LSMs that implement plant traits as parameters on the PFT level cannot capture this intrinsic PFT variability resulting from these traits. Although optimized parameters might still reduce the bias on the continental level, a more comprehensive approach to the spatiotemporal variability in plant traits might improve regional simulations drastically (Anderegg et al., 2022; Van Bodegom et al., 2014; Kattge et al., 2011).

Given the hydraulic role of vegetation leaves in controlling transpiration, there is a tight relationship between ET, GPP, and LAI. In CLM5–BGC, the assimilated carbon by GPP gets further partitioned to respiration and the carbon storage in the plant organs, i.e., leaves, roots, and stems. Furthermore, the leaf carbon then controls the development and state of the vegetation leaves and, thus, the LAI. On the other hand, LAI controls GPP by determining the upscaling factors from leaf photosynthesis to the canopy, thereby driving canopy conductance. Unfortunately, no large-scale LAI in situ measurements and no CLM5_PFT simulated LAI are available, and comparisons between CLM5_grid LAI and reanalysis or remote-sensing LAI suffer from known biases in the latter and yield no further context for our evaluation based on ground truth information. We adhered to an LAI evaluation of CLM5 with sparse but systematic ICOS measurements, ERA5L reanalysis, and GLASS based on MODIS (Fig. S15). Notably, the ICOS LAI measurements are only available for 2 years of the study period (2017 to 2018) and are limited to ENF and CRO sites. Additionally, LAI measurements' expensive and time-intensive nature restricts the time resolution to a few yearly measurement points. As a result, the data points for comparison are few, and the uncertainties are larger (noticeable larger error bars in Fig. S15). Another caveat is the potential mismatch of the land surface representation between the EC tower footprint (ET and GPP measurements) and the area covered by the LAI measurement campaigns. However, some key findings from this analysis are still robust. For example, all models overestimate LAI at ENF and CRO sites (Fig. 14a), contrasting the results of GPP and ET. The variance in ENF sites is much more significant in GLASS and ERA5L than in CLM5_grid, which is closest to the observations. The higher-order moments are more uncertain because of the small number of data points. The contrasting results, especially between the LAI and GPP PFT level averages, suggest that processes and parameters connecting the assimilated carbon to the leaf area, depending on the environmental conditions, must be revisited. However, we make an even stronger case for systematic, long-term, and high-resolution LAI in situ measurements (for example, using drones; Bates et al., 2021), which would support a more robust and diverse evaluation of the simulations of this essential variable.

4.3 Inter-site similarity of PFT groups

For all models (CLM5_grid, CLM5_PFT, ERA5L, GLASS, and GLEAM), the distributions of ET and GPP across PFTs are very similar, which is not the case for the observations. This is especially true for their variances (i.e., their spread around the mean) but also notable for the means, skewnesses, and kurtoses. We expected CLM5_PFT to show more significant variability than CLM5_grid and the other grid scale models because the aggregated mixed-PFT data of the grid cell would homogenize the variables and cancel out some of the variability. While CLM5_PFT shows a more extensive range of variation in ET and GPP across PFTs than CLM5_grid, ERA5L, GLASS, and GLEAM, it still vastly underestimates the observed range of variance by ICOS, especially for GPP (Figs. 9, 11).

The mean RMSD across sites of the same PFT indicates that ET across sites can be as different in CLM_PFT for GRA and CRO as in the observations (Fig. 12). However, the ET differences across sites with the same PFT were underestimated at ENF and DBF. GPP differences across sites with the same PFTs were underestimated for all PFTs (Fig. 13). This suggests that the missed variance could mainly stem from missed PFT internal inter-site differences or unresolved differences in site-specific abiotic conditions (e.g., soil depth and texture). Possibly, this could not be improved through optimization of PFT-specific parameters, as these sites would still share the same set of parameters. An enhanced concept of functional types in vegetation, focusing on the spatiotemporal variability in observed plant traits, could better facilitate improvements that raise the simulated ET and GPP variance in space and time.

4.4 Data requirements

As outlined above, beyond parameter optimizations, a comprehensive implementation of functional ecosystem diversity could significantly improve the LSM simulation outputs regarding multiple aspects of their distributions. This could introduce a state-of-the-art understanding of vegetation function into LSMs, which is essential to evaluate different theories of plant trait evolution and their effect on current and future energy, water, and carbon cycles.

In that light, we encourage sites to co-locate research infrastructures (Futter et al., 2023), like ICOS and the integrated European Long-Term Ecosystem Research Infrastructure (eLTER RI). Thereby, sites cover additional observation spheres like biodiversity (e.g., functional diversity of plants) and socio-ecology (through forest and crop management and driving land use change) and establish a strong base for studies to increase the understanding of the whole system (Ohnemus et al., 2024; Mirtl et al., 2018, 2021; Baatz et al., 2018). Furthermore, this would promote large-scale observations needed to introduce more trait variability into LSMs. Last, combining LSMs and these holistic observations by data assimilation, going beyond decoupled modeling efforts (Bloom et al., 2020) and resulting in an ecosystem reanalysis (Baatz et al., 2021), would provide essential, explicit and accurate data on the carbon cycle, which are currently unavailable.

4.5 Distribution moments and droughts

Investigating the influence of drought on the analyses, or generally the ability of the models to simulate drought and the vegetation response, is complex due to the differences in drought response functionality. For instance, plant water stress might occur due to different magnitudes of water deficit in the soil on different aggregation timescales and with a variable lag to the water deficit. A future study will investigate the PFT-scale drought responses from the model and how the drought propagates through the ecohydrological sphere and compare it to observations. However, drought frequency, duration, and severity affect the shapes of the distribution of the precipitation and, eventually, the ecosystem processes. Thus, we briefly discuss possible insights into their drought responses.

Importantly, the skewness and excess kurtosis moments, which inform us about the characteristics of the distribution tails (relativity between the tails and the general tailed aspect, respectively) of precipitation (Guo, 2022), as well as vegetation states and function (Kanavi et al., 2020; Liu et al., 2022; Cooley et al., 2022), are influenced by dry conditions, depending on their frequency, duration, and severity. We found a low variability in the skewness and excess kurtosis of the precipitation used to force our CLM5 simulations (Fig. S10c and d), specifically a significantly lower skewness and excess kurtosis at ENF and DBF sites. A lower positive skewness than the observations means that the distribution is less skewed towards lower values, and a lower positive excess kurtosis than the observations indicates generally larger tails. A possible interpretation of these differences in the distribution moments is that the atmospheric forcings show more frequent, longer, and more severe extreme precipitation events, while the ICOS measurements are more concentrated around their mean. While the propagation of these extreme events could be complex and non-linear, we generally found the same results (lower skewness and smaller absolute excess kurtosis) for the simulated distributions of ET and GPP for almost all models and PFTs (Figs. 9c–d and 11c–d), suggesting a more direct relationship. However, because of the possible non-linearity and the influence of other factors, the detailed relationship between these findings and the ability of CLM5 to simulate ecosystem drought responses must be examined in future studies. In any case, the missing accuracy in representing higher distribution moments in the atmospheric forcings and in land surface models must be considered in studies using these to investigate drought.

5 Conclusions

We evaluated the simulated evapotranspiration (ET) and gross primary production (GPP) from a 3 km resolved Community Land Model v5 (CLM5) set up over the European CORDEX domain. We differentiated the model outputs between the grid scale (CLM5_grid) and the plant-functional-type scale (CLM5_PFT) and compared them with ICOS station data as ground truth data. Furthermore, we compared with ET and GPP from remote-sensing-derived data from the Global Land Surface Satellite (GLASS) and reanalysis products such as the European Centre for Medium-Range Weather Forecast Reanalysis 5 – Land (ERA5L) and the Global Land Evaporation Amsterdam Model (GLEAM). CLM5_grid and CLM5_PFT exhibited promising skills in approximating the observations and often performed better than ERA5L, GLASS, and GLEAM. CLM5_PFT showed a lower systematic bias (lower percent bias) but approximated the ICOS observations in a generally worse manner (larger root mean square error) than CLM5_grid (Figs. 2 and 3 and Tables S4 and S5). ET and GPP were systematically underestimated across all PFTs throughout the year for both model scales. Especially during summer at DBF sites, GPP was substantially lower for CLM5_PFT and CLM5_grid than for ICOS observations (Figs. 4, 6).

Essentially, CLM5_PFT and, to a greater degree, CLM5_grid, ERA5L, GLEAM, and GLASS showed a lower spatiotemporal variability in ET and GPP than the measurements exhibited by a lower range of all the modeled ET and GPP distribution moments across PFTs than in ICOS. This smaller range and a lower root mean square difference between sites of one PFT group suggests that CLM5_grid and, more surprisingly, CLM5_PFT simulate GPP and ET more similarly across PFTs than the ICOS measurements.

Future studies should investigate whether optimizing parameters in CLM5_PFT with observation data increases the diversity of ET and GPP values or whether this is a structurally induced bias. This work provides essential insights for studies that aim to find optimized parameters and meaningful context for analyses of more specific ET and GPP dynamics using the evaluated data.