Unrecognised water limitation is a main source of uncertainty for models of terrestrial photosynthesis

Quantification of environmental controls on ecosystem photosynthesis is essential to understand the impacts of climate change and extreme events on the carbon cycle and the provisioning of ecosystem services. Machine learning models have become popular for simulating ecosystem terrestrial photosynthesis because of their predictive skill, but often do not consider temporal dependencies in the data, even though process understanding suggests that these should exist. Here, we investigate how models that account for temporal structure impact the prediction of ecosystem photosynthesis. Using time-series measurements of ecosystem fluxes paired with measurements of meteorological variables from a network of globally distributed sites (N=104) and remotely sensed vegetation indices, we train three different models to predict ecosystem gross primary production (GPP): a mechanistic, theory-based photosynthesis model, a memoryless multilayer perceptron (MLP) and a recurrent neural network (Long Short-Term Memory, LSTM). Through comparisons of patterns in model error, we assess the ability of these models to predict GPP across a wide diversity of ecosystems and climates, and to account for temporal dependencies, with a focus on effects by low rooting zone moisture and freezing air temperatures. While both deep learning models outperform the mechanistic model, we find their overall performance is similar, with an R² of 0.79 spatial out-of-sample predictions for both models. Overall, model skill is consistently good across moist sites with strong seasonality. During periods affected by temporal patterns such as drought and frost, the LSTM shows lower model error than the MLP as well as an LSTM with shuffled input, showing that there is an advantage from learned temporal dependencies. Generalisation patterns reveal that the LSTM tends to be more successful than the (time-agnostic) MLP in simulating GPP in dry environments. However, a large variability in model skill across relatively dry sites remained. This was not resolved by the inclusion of additional earth observation data, although this improved overall model performance. Insufficient information on the exposure and response to water stress and related effects on GPP appear to be dominant sources of error for modelling ecosystem fluxes across the globe. With the increasing frequency of hydroclimatic extreme events, effects of water limitation are expected to become more prevalent, which calls for models that better represent its impact on ecosystem function.

1 Introduction

Photosynthesis plays a major role in the global carbon cycle and drives important ecosystem functions (Beer et al., 2010). Ecosystem-level gross CO₂ uptake through photosynthesis is referred to as gross primary production (GPP) and varies in response to the environment. Understanding its variations across space and time as well as its dependencies on environmental conditions is key for predicting changes and feedbacks in the terrestrial biosphere (Booth et al., 2012).

GPP variations are driven by multiple, simultaneously varying environmental factors and the physiological and structural responses of plants to these conditions. Solar radiation supplies the energy for photosynthesis and acts as a dominant driver of GPP, depending on light absorption (Monteith, 1972). Temperature, light and water availability trigger phenological changes and regulate seasonal cycles of active leaf surface area and therefore seasonal changes in light absorption. Air temperature affects leaf temperatures, which in turn govern enzymatic rates and photosynthesis (Berry and Bjorkman, 1980; Kattge and Knorr, 2007; Kumarathunge et al., 2019; Bernacchi et al., 2003). Low moisture availability across the rooting zone, in combination with a high vapour pressure of air at the leaf surface, determines the effects of water stress and can lead to GPP reductions (Stocker et al., 2018; Novick et al., 2016).

Continuous GPP estimates can be obtained from eddy covariance measurements of ecosystem gas exchange (Baldocchi, 2020) and capture surface-atmosphere exchange fluxes, integrated over a radius on the order of a kilometre around the site of measurement (Chu et al., 2021). These measurements, paired with observations of meteorological variables and soil conditions, are made available through different networks and initiatives (e.g., AmeriFlux, 2025; ICOS, 2025; OzFlux, 2025). The combination of data from multiple regional networks has led to large datasets with standardized processing of eddy covariance measurements from a large number of sites (Pastorello et al., 2020; Hufkens and Stocker, 2025; Abramowitz et al., 2024). The availability of large datasets of GPP along with their environmental covariates, and paired with remotely sensed variables, has made machine learning (ML) a widely used approach for predicting spatio-temporal variations of ecosystem-atmosphere exchange fluxes (Kang et al., 2025; Gaber et al., 2024; Tramontana et al., 2016; Montero et al., 2024; Papale et al., 2015; Yang et al., 2007; Jung et al., 2011; Joiner and Yoshida, 2020; Zheng et al., 2020).

Process understanding and empirical patterns of GPP dynamics suggest that there should be temporal dependencies in data of GPP and its predictors. Temporal dependencies arise as a result of several processes. Low soil moisture can reduce GPP (Stocker et al., 2018) and reflects the history of precipitation, radiation, and leaf phenology over the preceding weeks to months. Plant hydraulic processes induce a temporal hysteresis effect over the course of diurnal cycles (Tuzet et al., 2003). Physiological changes are caused, e.g., by the seasonal acclimation of the photosynthetic apparatus to varying levels of radiation inputs and temperature (Kumarathunge et al., 2019; Luo and Keenan, 2020; Liu et al., 2024b; Berry and Bjorkman, 1980). Ecosystems in cold climates have been found to delay springtime photosynthesis resumption early in the season through photoprotective processes, despite high levels of solar radiation (Luo et al., 2023, henceforth referred to as “cold acclimation”). Stress by extreme environmental conditions can cause delayed and long-lasting effects, such as impaired transpiration and reduced CO₂ assimilation (Barber and Andersson, 1992; Reichstein et al., 2013; McDowell et al., 2022; Bastos et al., 2020; Yu et al., 2022).

Several published machine learning models for GPP treat values of GPP time series as independent and identically distributed observations and therefore do not account for temporal dependencies (Nelson et al., 2024; Kang et al., 2025; Tramontana et al., 2016; Gaber et al., 2024). This limitation may be relieved by temporal aggregation to daily-monthly time scales and by pairing data with additional, remotely sensed observations that capture phenological changes and variations in the amount of active, light-intercepting foliage area (Baldocchi, 2018). However, additional physiological changes that affect the efficiency of light utilization for CO₂ assimilation at the leaf level are more challenging to capture by remotely sensed reflectance data (Ryu et al., 2019; Stocker et al., 2018). As a consequence, substantial unexplained GPP variation is expected to remain at the seasonal and diurnal time scales.

A potential solution for this problem is the use of time-aware machine learning algorithms that can learn non-stationary relationships and temporal dependencies. Such algorithms have been introduced for modelling GPP and related fluxes (Nakagawa et al., 2023; Besnard et al., 2019; Kraft et al., 2025; Montero et al., 2024). Montero et al. (2024) compared the performance of three recurrent architectures for GPP modelling and evaluated them on GPP extremes. Besnard et al. (2019) evaluated a Long Short-Term Memory (LSTM) network to study net ecosystem CO₂ exchange. Kraft et al. (2025) assessed sequential models for global upscaling of evapotranspiration. Nakagawa et al. (2023) introduced a temporal fusion transformer for global upscaling of GPP. In these previous studies, the impact of using such an architecture for modelling known temporal effects was either not evaluated or inconclusive.

In contrast to ML models, mechanistic GPP models, such as the P-model (Stocker et al., 2020), embody process understanding and provide a theory-based prediction that accounts for these known temporal dependencies. The foundation in plant physiology may also help these models to generalise more robustly, as the underlying relations remain valid when extrapolating to new conditions not seen in the training data. The price to pay is that these models lack the flexibility to pick up any patterns that were not anticipated during their design, whereas the high representation power of ML models gives them the ability to uncover and respect such patterns.

While the P-model only uses a single greenness index derived from remote sensing, previous work on flux modelling has shown that additional remotely sensed variables can be informative (Nelson et al., 2024; Kraft et al., 2025). The signal from thermal remote sensing may reflect changes in photosynthesis that are driven by physiological responses and stomatal regulation, affecting transpiration and therefore surface energy partitioning and surface heating. Therefore, land surface temperature (LST) may be useful information for GPP prediction. Common mechanistic and light use efficiency models don't consider this information as additional forcing. Furthermore, the full information of surface reflectance in all available individual bands may contain additional information about GPP changes as pigments deployed under stress conditions or the leaf water content can affect surface reflectance beyond what commonly used single greenness-based indices reflect (Ceccato et al., 2001; Gamon et al., 2016).

In this study, we evaluate the use of an LSTM (Hochreiter and Schmidhuber, 1997) as a predictor of GPP. LSTMs have been shown to be successful at tasks where memory effects across a range of temporal scales are involved, such as sea surface temperature prediction (Zhang et al., 2017), rainfall-runoff modelling (Kratzert et al., 2018) and canopy greenness modelling (Liu et al., 2024a). To contrast the recurrent and purely data-driven design, we compare against a standard, non-recurrent multilayer perceptron (MLP), an LSTM with permuted input, as well as to the process-based P-model. The P-model serves as a benchmark with known treatment of temporal effects. We compare these models based on model performance and generalisation capabilities. To investigate the ability of the LSTM to account for temporal dependencies, we assess seasonal patterns of cold acclimation and water limitation effects in dry conditions. Additionally, we analyse spatial patterns of model generalisability (spatial out-of-sample performance) with respect to different environmental factors.

In addition to the standard set of predictors, we include LST and multiple bands of surface reflectance. To aid the models in simulating temporal dependencies, we also provide additional features related to soil water availability. In view of the known influence of root zone moisture on GPP (Stocker et al., 2018) and the inability of the MLP to account for the precipitation and radiation history (and thus implicitly for the evolution of root zone moisture), we test if its performance improves when (simulated) soil moisture is provided as a complementary predictor. Providing observation-derived soil moisture or a general index of water availability as an additional predictor is a common approach taken also for other memoryless GPP models (Nelson et al., 2024; Kang et al., 2025; Tramontana et al., 2016; Gaber et al., 2024). While the LSTM is expected to be able to learn the effects of soil moisture limitation, we also test whether it benefits from an Earth-observation derived estimate of root zone water storage capacity (Stocker et al., 2023) as additional time-invariant context.

2 Materials and methods

2.1 Data

We sourced daily GPP data from a collection of eddy covariance flux sites gathered from the PLUMBER2 framework (Ukkola et al., 2022), which includes sites from OzFlux (Isaac et al., 2017), FLUXNET2015 (Pastorello et al., 2020), LaThuile (2025) as well as AmeriFlux (2025), ICOS Warm Winter 2020 (Warm Winter 2020 Team and ICOS Ecosystem Thematic Centre, 2022) and ICOS Drought 2018 (Drought 2018 Team and ICOS Ecosystem Thematic Centre, 2020). Site selection was performed through several steps. Sites located in cropland or wetland ecosystems were excluded. For each site, we only included full calendar years of data. Some years of data were excluded due to evident inconsistencies found by visual inspection. We selected sites with at least five consecutive years of high-quality, gap-free data. GPP data were included if at least 50 % of all half-hourly measurements were of good quality (either measured or gap-filled with high confidence).

We used GPP estimates generated by the nighttime partitioning method (GPP_NT_VUT_REF; Reichstein et al., 2005). Half-hourly GPP estimates were aggregated to obtain daily GPP values. In addition to GPP, meteorological variables were obtained, which were measured directly at the flux sites. We used the following meteorological variables: air temperature (TA_F_MDS), daytime air temperature (TA_DAY_F_MDS), shortwave incoming radiation (SW_IN_F_MDS), longwave incoming radiation (LW_IN_F_MDS), daytime vapour pressure deficit (VPD_DAY_F_MDS), air pressure (PA_F), precipitation (P_F) and wind speed (WS_F). Observations that were either missing or had insufficient quality (< 50 % measured or gap-filled half-hourly measurements with high confidence) were gap-filled with linear interpolation for air temperature, respectively k-nearest neighbour imputation per site. This only affected around 10 % of all predictor values.

Along with local site-level measurements, we used remotely sensed estimates of the fraction of absorbed photosynthetically active radiation (fAPAR), extracted from the MODIS FPAR MCD15A2H Collection 6.1 product (Myneni et al., 2021). fAPAR captures variations in phenology and represents the amount of solar radiation absorbed by the canopy and usable for photosynthesis. fAPAR data were extracted for the pixel (500 × 500 m² area) that contains the flux measurement site and for the eight pixels immediately surrounding it. The nine values where combined through a weighted average, using as weights the inverse of their variance as per the data product. The fAPAR sequences were gap-filled based on the mean seasonal cycle, then smoothed and interpolated to the time resolution of the flux data with a LOESS spline.

We retrieved additional pre-processed remote sensing variables from the MODIS instruments through the FluxnetEO dataset (Walther et al., 2022). We included the first seven reflectance bands (RED, NIR, BLUE, GREEN, SWIR1, SWIR2, SWIR3) as well as daytime and nighttime land surface temperature (LST) at a viewing zenith angle of 0° (LST_TERRA_Day_VZA0, LST_TERRA_Night_VZA0).

The collection of flux sites was further filtered to include only sites with data available from the FluxnetEO dataset. This process resulted in a collection of 104 sites (Fig. 1) with a total of 1020 site-years of data. Detailed site information is given in Appendix A. The selected sites cover a wide range of environmental factors, in particular aridity (Figs. 1 and D1). To measure aridity, we calculated the moisture index (MI) for each site as total $P/\mathrm{PET}$, where P is the precipitation measured at the site and PET is the potential evapotranspiration following Priestley-Taylor (Priestley and Taylor, 1972), as implemented in the SPLASH ecosystem water balance model (Davis et al., 2017). The root zone water holding capacity for each site was sourced from Stocker et al. (2023).

Figure 1 Site locations and their moisture indices. For each category, the quantised moisture index values and number of sites are given in brackets.

2.2 Models

We have implemented two deep learning models to evaluate different architectures for the prediction task. To account for temporal dependencies, we used a Long Short-Term Memory (LSTM) network (Hochreiter and Schmidhuber, 1997). The network included LSTM cells with layer normalization, known to stabilize hidden dynamics and reduce training time (Ba et al., 2016). The number of layers and hidden dimension of the network were tuned as hyperparameters on validation data. The LSTM layers were followed by a variable number of linear layers with GELU activations (Hendrycks and Gimpel, 2023) that each halve the dimension until reaching a size of 16 neurons.

To isolate the impact of recurrence, we implemented a second neural network model without any memory mechanism, namely a standard multilayer perceptron (MLP). Its architecture is identical to keeping only the linear layers of the LSTM. The hidden dimension of the first layer was tuned as a hyperparameter. The architectural similarity means that the difference between predictions of the two networks is a good indication for the influence of information from past time steps. Hyperparameters were tuned separately to ensure both variants reach their best performance. We also evaluated a version of the LSTM where we permute the temporal ordering (LSTM_perm). This model has the same architecture and capacity but it cannot rely on memory as the temporal patterns are eliminated from the input data. This version serves as an additional evaluation that disentangles the influence of the memory of the LSTM from other model properties such as the overall capacity.

The third model was the P-model, a mechanistic, theory-based representation of ecosystem-level photosynthesis acclimation and GPP (Stocker et al., 2020; Wang et al., 2017; Prentice et al., 2014). It builds on the widely used Farquhar-von Caemmmerer-Berry (FvCB) model for leaf-level C₃ photosynthesis (Farquhar et al., 1980). The FvCB model is combined with an optimal balancing of the costs of carbon assimilation and transpiration (Prentice et al., 2014). Furthermore, the P-model implements the coordination hypothesis, which states that photosynthesis is balanced at the intersection of light and Rubisco-limited assimilation rates during average daytime conditions (Maire et al., 2012). Based on these relations, the P-model predicts photosynthesis acclimation parameters to describe the processes that determine the light use efficiency (LUE). GPP is then modelled as the product of LUE and absorbed photosynthetically active radiation (APAR), which in turn is taken to be the product of the photosynthetic photon flux density (PPFD) and fAPAR. The forcings for the P-model correspond to the input data of the LSTM and the MLP. We used the FULL model setup as described in Stocker et al. (2020), which includes an empirical soil moisture stress function and temperature dependency of the intrinsic quantum yield. These two components were calibrated to the data with an optimization of four parameters through minimising the root mean squared error with the generalised simulated annealing algorithm, as implemented in the GenSA R package (Xiang et al., 2013). The P-model is implemented in the R package rsofun (Stocker et al., 2024).

2.3 Experimental setup

We assessed the three models for their ability to handle temporal dependencies and their ability to generalise to new sites with different environmental conditions.

2.3.1 Global model

We first evaluate each model in a spatial cross-validation setup, which measures performance at sites that were not seen during model training (where “training” of the P-model means calibration). For the spatial cross-validation, we assigned each site to one of five folds, stratified based on the per-site mean air temperature and the moisture index to achieve a similar distribution of climate types in all folds. In turn, four folds served as training data to fit the model weights and tune the hyperparameters, then GPP predictions were produced for the test sites in the fifth, held-out fold.

2.3.2 Site-specific model

To separate changes in environmental conditions (which can be covered by the global model) from potential variations of the functional relationships between different sites (which cannot be represented by a single set of model parameters), we also fit separate per-site models and evaluate them with a temporal cross-validation. In that setup, the temporal sequences of predictors and GPP at every individual site are split into years, setting the start of the year to the coldest month for sites in temperate, continental and polar climates, and to the wettest month for tropical and arid sites (following the Köppen-Geiger climate code). Cross-validation then proceeds by holding out every year in turn, and fitting the model on the remaining years.

2.3.3 Model training

In both setups, the models are trained on chunks of 128 d, whereas testing was performed on the full sequences (i.e., all data of a site in the spatial cross-validation setup, respectively individual years in the temporal cross-validation setup). The training chunks were created per site, with a random start date of the initial chunk within the first 96 d and a sliding window with regular overlap of 32 d. The features were standardized using the mean and standard deviation of the training folds.

Models were trained by minimising the mean squared error with the Adam optimizer (Kingma and Ba, 2017). Training was performed for a maximum 50 iterations and stopped after 10 iterations without improvement in the validation loss. To limit overfitting, an L2 penalty was applied on the parameter updates, and dropout (Hinton et al., 2012) was used after each LSTM layer except for the last layer. The learning rate was adaptively reduced when the loss no longer improved for several iterations. The batch size, the weight of the L2 penalty, the dropout rate, the initial learning rate, the patience before reducing the learning rate as well as the reduction factor were all tuned in an inner cross-validation loop within each data fold. In the spatial cross-validation setup, a 3-fold inner cross-validation was performed within each of the 5 data folds. The sites were distributed based on mean air temperature and moisture index to ensure equal representation of climate types across folds. The hyperparameters were tuned using random sampling of 20 configurations. In the temporal cross-validation setup, hyperparameters were tuned for each held-out test year based on the remaining years of data for each site. A cross-validation with 5 randomly chosen held-out validation years was performed. In this setup, the hyperparameters were tuned using random sampling of 40 configurations. In both setups, the configuration with the lowest average validation RMSE across the held-out folds was selected. The set of options for each hyperparameter is listed in Table 1. After hyperparameter tuning, models were trained on the full data folds with the chosen hyperparameter configurations, with 20 % of the sites used as validation data for early stopping of the training.

2.3.4 Evaluation

GPP predictions were evaluated using the squared Pearson's correlation coefficient (R²) and root mean squared error (RMSE). In addition to assessing the daily predictions, we aggregated predictions and observations to different scales. We calculated the mean seasonal cycle by averaging over all years observed at a site to obtain a mean value per day of the year. Predictions and observations were also temporally aggregated to site-level means. Moreover we calculated daily anomalies, defined as deviations between the daily values from the mean seasonal cycle; as well as yearly anomalies, defined as deviations between a site's annual mean values and its global, multi-year mean.

The site-specific models were evaluated for test years that start at the wettest or coldest month. For days before the first day of the first such month, no predictions where made. When comparing the site-specific model and the global model at the site level, we therefore filter the predictions of the global model to comprise the exact same test days as the site-specific model.

When evaluating the P-model, the same spatial cross-validation was used as for the machine learning models, with the model parameters calibrated separately for each fold (Stocker et al., 2024).

We investigated the ability of the models to capture two different, well-known temporal effects. To test how well soil moisture effects are reproduced, we looked at the (absolute) percentage error of the model predictions as the potential cumulative water deficit (PCWD) increases. PCWD was calculated from the flux data as the cumulative difference between potential evapotranspiration (PET) and precipitation (Stocker, 2021). PET was estimated based on Priestley-Taylor (Priestley and Taylor, 1972), as implemented in the SPLASH ecosystem water balance model (Davis et al., 2017). Values were pooled from all sites and all test days per PCWD interval, with intervals chosen such that they have at least 100 data points. We evaluated the temporal effect from soil moisture from a different angle through an analysis of model bias during drought events identified from the full dataset. Drought events were collected through the identification of sites and days where light use efficiency was reduced for at least 3 d (Stocker et al., 2018). All events were aligned at the onset. The bias was calculated as the difference between model predictions and observations and aggregated per day for the period between 10 d before the start of each event and up until 100 d after the start of each event. The bias values were normalised by subtracting each value by the median value during the window covering between 1 and 10 d before each event.

Second, we evaluated the models' ability to reproduce cold acclimation effects. We selected four sites that have been found to have a reduced light use efficiency and thus a delayed increase in GPP at the start of the growing season (DE-Hai, US-Ha1, US-MMS, US-PFa). Prediction errors at these sites were contrasted with those at four sites that did not exhibit any GPP delay (BE-Vie, FI-Hyy, NL-Loo, RU-Fyo), based on the findings of Luo et al. (2023). For the two groups of sites, we aggregated and compared predictions per day of the year with different models.

For the global LSTM model, we compared site-level performance across various environmental conditions: moisture index ($P/$PET), Köppen-Geiger climate zone (Beck et al., 2018), and IGBP vegetation type (International Geosphere-Biosphere Programme). To further investigate the generalisation of the models across sites, we compared site-level performance between the global model and the site-specific model, by computing the relative difference $\mathrm{\Delta }{R}^{\mathrm{2}}=R_{\mathrm{global}}^{\mathrm{2}}-R_{\mathrm{site}}^{\mathrm{2}}$ and the ratio $\mathrm{rRMSE}={\mathrm{RMSE}}_{\mathrm{global}}/{\mathrm{RMSE}}_{\mathrm{site}}$.

We evaluated two different sets of predictor variables. The standard set included the site-level meteorological measurements as well as fAPAR. The expanded set included LST variables and the MODIS reflectance bands. The model evaluations with this expanded set are denoted by a + after the model name. We compared this setup to the standard predictor set in terms of overall performance. Separately to the standard setups, we fed additional predictors to the feature set of the deep learning models. For the MLP, we added soil moisture. Due to the limited quality of measured soil moisture at many flux sites, we used modelled soil moisture from the SPLASH water balance model (Davis et al., 2017). For the LSTM, we added the root zone water storage capacity, extracted from the global map of Stocker et al. (2023).

Finally, we evaluated a Temporal Convolutional Network (Bai et al., 2018) and an LSTM with attention layer (Vaswani et al., 2017) as two alternative sequence models (Appendix C). We compared their overall performance as well as their event response performance to the LSTM.

The details of the data, models and experimental setup are summarised in Tables B1, B2 and B3.

Table 1 Hyperparameter search space for the LSTM and MLP models.

3 Results

3.1 Overall performance

Overall, we found that both machine learning models in the global setting (i.e., a single, fixed model trained on multiple sites) predict GPP more accurately than the process-based P-model (Table 2). Furthermore, the model variants with the additional remote sensing variables (+ version) performed better than the models with only the standard predictor set. For daily predictions, the R² calculated from pooled data of all sites was 0.79 for both the LSTM+ and the MLP+, compared to 0.64 for the P-model. The RMSE was 22.5 % lower for the LSTM+ than for the P-model. Both the LSTM+ and MLP+ modelled the seasonal cycle well, with R² values of 0.88. The P-model achieved an R² of 0.78 for modelling the mean seasonal cycle. The prediction of anomalies was more challenging for all three evaluated models. The LSTM+ achieved an R² of 0.28 for daily anomalies and R² of 0.1 for annual anomalies. While the differences in both R² and RMSE were minimal between the LSTM+ and MLP+, the P-model was outperformed in both aspects.

Model performance varied substantially between sites (Fig. 2). For 94 (out of 104) sites the LSTM+ reached a higher R² value than the P-model. The LSTM+ also outperformed the MLP+ at slightly more than half of the sites. While the LSTM+ clearly outperformed the LSTM in terms of overall metrics, the performance didn't improve at every site. At 68 out of 104 sites, the LSTM+ outperformed the LSTM. In further analyses, we focus on the + versions of each model, as the overall performance evaluation indicates that they are better predictors of GPP.

Table 2 Performance metrics (R² and RMSE) of the global model at different aggregation levels. Metrics are calculated from pooled data of all sites. “Daily” evaluates daily predictions and observations. “Seasonal” aggregates by day of year per site. “Spatial” uses site means. “Daily anom.” is the deviation from each site’s mean seasonal cycle. “Annual anom.” is the deviation of the annual mean from the multi-year mean per site. The best values per row are printed in bold.

Figure 2 Comparison of the LSTM+ against the P-model (a), against the MLP+ (b), and against the LSTM (c). R² of predicted versus observed values of daily GPP per site from the global cross-validation is shown for the LSTM+ along the y-axis and the P-model/MLP/LSTM along the x-axis. The dotted line indicates equal performance.

More salient differences between the models were observed when inspecting the predicted mean seasonal cycles within different climate zones (Fig. 3). The deep learning models were better at predicting the timing of early spring GPP increase in several climates (Köppen-Geiger codes Dfb, Dfc, Cfa). They also outperformed the P-model in desert and semi-arid climates.

Figure 3 Mean seasonal cycle of GPP and model predictions by climate zone and hemisphere. Predictions from the global cross-validation for the LSTM+, LSTM_perm+, MLP+ and P-model are compared against GPP observations. Climate zone boundaries are from Beck et al. (2018).

3.2 Temporal patterns in model error

From the overall performance metrics, it appeared that the LSTM+ is not clearly better than the MLP+, despite its ability to learn temporal patterns. However, the improved performance of the LSTM+ compared to the LSTM_perm+ indicates that temporal dependencies were learned. Targeted evaluations of temporal patterns give a clearer insight into the differences between the MLP+, the LSTM_perm+ and the LSTM+.

The LSTM+ showed different error characteristics than the P-model at high values of PCWD. While relative errors increased with higher PCWD for the P-model, they stayed relatively constant for the LSTM+ (Fig. 4a) For all pooled data, the four models showed similar error distributions at lower levels of PCWD, again the errors of the P-model increased at higher PCWD values, whereas they did not for the LSTM+, LSTM_perm+ and the MLP+. Separating the analysis of model errors versus PCWD by vegetation type revealed differences at higher PCWD values (Fig. 4b). For evergreen forests, relative errors increased for the P-model, the MLP+ and the LSTM_perm+ from a PCWD of 800 mm, but also decreased again. The relative errors increased less for the LSTM+ than for the other models. For non-evergreen forests, both the LSTM+ and the MLP+ showed lower relative error than the P-model above a PCWD of 1000 mm.

Figure 4 Error distribution across different amounts of cumulative water deficit (PCWD, quantized into 100 mm bins). Bins with less than 100 samples are not shown. Predictions are by the global models. Lines denote the median values, shaded regions lie between the lower (25 %) and upper (75 %) quartiles. (a) all data pooled. (b) Evergreen sites, including Evergreen Needleleaf Forests and Evergreen Broadleaf Forests. (c) Non-evergreen sites, including all other vegetation types (Deciduous Forests, Shrublands, Savannas and Grasslands).

The comparison between sites with and without cold acclimation (delayed GPP) revealed clear differences w.r.t. the predicted seasonal cycles (Fig. 5) of the different models and the seasonal cycle of model bias. For sites without a delay in springtime GPP increase, the LSTM+ performed best during spring and summer. For sites with cold acclimation, the deep learning models capture the delay better than the P-model. The LSTM+ predicts the evolution of GPP best during springtime, although that edge is mostly after the start of spring, whereas some bias remains at the onset of GPP increase as well as during summer.

Figure 5 Mean seasonal cycle of GPP and model predictions of the LSTM+, LSTM_perm+, MLP+ and P-model for four sites with delayed GPP (a) and four sites without delayed GPP (b), as well as mean bias per day of the year for sites with delayed GPP (c) and non-delayed GPP (d). The bias is calculated as the difference between model predictions and observations.

Figure 6 The model bias per day before and after the start of drought events. The bias is calculated as the difference between model predictions and observations. The bias is normalised with respect to the median value during the 1 to 10 d before each event onset. Shaded areas mark the area between the 33rd and 66th percentiles.

The event analysis (Fig.  6) showed that the deep learning models predicted the drought response equally well during the first 20 d after an event. The LSTM+ predicted the response best from 20 d after the start of a drought event.

Evaluations of the alternative sequence models (Appendix C) showed that while the TCN+ performed better at 60 out of 109 sites (Fig. C1), no difference was found between models in terms of their performance during drought events (Fig. C2), supporting the further evaluation of the LSTM+ for temporal patterns.

3.3 Spatial patterns in model performance

The observed error patterns suggest a qualitatively different behaviour of machine learning models, especially the LSTM+, during conditions where temporal effects are known to occur. To investigate this further at the site level, we plot model performance per site against relevant site characteristics (Fig. 7).

Figure 7 Performance as R² (a) and normalised RMSE (b) of the LSTM+ (global model) per site against different site characteristics. Section (a) and (b) each show comparisons of model performance against different site characteristics: moisture index, vegetation type and climate zone. The plots on the left show model performance per site against the moisture index. The plots on the right show model performance per site for different categories of vegetation type and climate zone, as well as box plots summarising the performance in each category. GRA: grassland, SAV: savanna, EBF: evergreen broadleaf forest, WSA: woody savanna, MF: mixed forest, CSH: closed shrubland, ENF: evergreen needleleaf forest, DBF: deciduous broadleaf forest, OSH: open shrubland.

The LSTM+ performed best for relatively moist sites. Across sites with moisture index $P/$PET ≥,0.75 the median R² is 0.82, whereas it was only 0.65 for more arid sites (MI < 0.75). The (normalised) RMSE follows a similar pattern, with a median value of 0.62 for sites with MI < 0.75, compared to 0.44 for moist sites (Fig. D2).

We found only a single site with poor predictions and a high moisture index, conditions at that site are likely not adequately represented since it is one of the few very moist evergreen sites in the dataset. All other sites with low performance (R²<0.5) have a moisture index below 0.75. The opposite is not true: the R² varies greatly among the drier sites, not all sites with low moisture index exhibit poor performance.

All sites with low R² (< 0.5) belonged to desert and semi-arid (BWh, BWk, BSh, BSk), Mediterranean (Csa, Csb) and hot and warm temperate (Cfa, Cfb) climates. In terms of normalised RMSE, predictions are especially poor in cold semi-arid climates (BSk). Also sites in polar climates (ET) performed relatively poorly, followed by other desert and semi-arid climates.

Among vegetation types, the LSTM+ works best for sites with deciduous broadleaf forest vegetation, with an average R² of 0.84. Sites with inaccurate predictions (R²<0.5) are spread across most vegetation types, with the exception of closed shrublands and deciduous broadleaf forest. The vegetation types with the lowest overall performance are evergreen broadleaved forests and open shrublands.

3.4 Generalisation across space

The ability of the global LSTM+ to generalise across space varies across different values of the moisture index, vegetation types, and climate zones (Fig. 8). For most sites, we observe a positive rRMSE, meaning that the global model is less accurate than the site-specific ones. On the other hand, ΔR² exhibits a less clear dependence on the moisture index, vegetation types, and climate zones.

The LSTM+ appears to show worse generalisability (more negative ΔR² and rRMSE further above 1) at relatively arid sites (MI < 0.75), compared to more moist sites. For relatively arid sites, site-specific models mostly have a lower RMSE than the global one (5 out of 40 sites, or 12.5 %), whereas the advantage is less pronounced for moist sites (21 out of 64, or 33 %). Out of the 13 sites with a considerable drop in performance of the global model compared to the site-specific model (poor generalisability of the model to those sites, $\mathrm{\Delta }{R}^{\mathrm{2}}<-\mathrm{0.1}$), 11 have a MI of less than 0.75. Out of those 11 sites, 4 are classified as a hot summer Mediterranean climate (Csa) and 6 as a semi-arid climate (BSh, BSk). For the Mediterranean, desert, and semi-arid climates, generalisation capabilities vary substantially, with the global model outperforming the site-specific one in some sites but not in others.

In terms of generalisation, there are no clear overall differences between the LSTM+ and the MLP+ (Fig. 9). However, a pattern emerges w.r.t. aridity: sites for which the global LSTM+ model performed better than the MLP+ have a mean moisture index of 0.85, while sites for which the opposite was found (global MLP+ model better than global LSTM+ model) have a mean MI of 1.11.

Figure 8 Patterns in the LSTM+'s ability to generalise, measured by ΔR² (a) and rRMSE (b). Positive ΔR² values and rRMSE values  < 1 mean the global model performs better than the site-specific one. Section (a) and (b) each show comparisons of model performance against different site characteristics: moisture index, vegetation type and climate zone. The plots on the left show generalisation ability per site against the moisture index. The plots on the right show generalisation ability per site for different categories of vegetation type and climate zone, as well as box plots summarising the performance in each category. GRA: grassland, SAV: savanna, EBF: evergreen broadleaf forest, WSA: woody savanna, MF: mixed forest, CSH: closed shrubland, ENF: evergreen needleleaf forest, DBF: deciduous broadleaf forest, OSH: open shrubland.

Figure 9 Per-site comparison of ΔR² LSTM+ vs. MLP+. Colours encode the moisture index. The number of sites on either side of the diagonal and their mean moisture index are displayed in the corners. At the top we show the significance of the difference in MI, as per the two-sided Mann-Whitney U-test.

3.5 Performance with additional features

Based on the observed temporal error patterns in response to water deficit, we tested whether the memoryless MLP+ would benefit from soil moisture as an added predictor (Fig. 10). Including soil moisture information from the SPLASH water balance model led to a minimal difference in overall performance. There is a trend that adding soil moisture improves the MLP+ prediction at drier sites (mean MI 0.93), but leads to a small performance loss at moist sites (mean MI 1.10).

We also tested the LSTM+ with the estimated root zone water holding capacity as additional (time-invariant) predictor (Fig. 10). While this led to minimal overall differences, the performance improved at both very moist and very dry sites.

Figure 10 Impact of complementary input features (soil moisture for the MLP+, water holding capacity for the LSTM+). (a) to (c) show the R² and (d) to (f) show the RMSE. (a) and (d) compare the MLP+ with soil moisture to the MLP+. (b) and (e) compare the LSTM+ to the MLP+ with soil moisture. (c) and (f) compare the LSTM+ with water holding capacity to the LSTM+. Each dot represents the performance for one site. The number of sites on either side of the diagonal and their mean moisture index are displayed in the corners. At the top we show the significance of the observed differences, as per the two-sided Mann-Whitney U-test.

4 Discussion

4.1 Key insights

We have compared different model types for predicting GPP from a shared set of predictors. To summarise, the following key results were obtained:

• The neural networks (LSTM and MLP) are better GPP predictors than the mechanistic model (P-model) (Table 2, Figs. 2, 3).

• The LSTM improves GPP prediction compared to the LSTM_perm (and MLP) under temporal stress patterns by leveraging learned memory (Figs. 4, 5, 6).

• When averaged across conditions, the LSTM is not clearly better than the MLP (Table 2). However, the LSTM shows advantages under drought conditions (Figs. 4, 6, 9). This can be mitigated to some extent by providing a simple soil moisture index (Fig. 10).

• With added earth observation variables, there is a clear improvement in model performance (Table 2, Fig. 2). However, there is still a large variability in model performance in dry conditions (Figs. 7, 8).

In the following, we discuss these insights in more detail.

4.2 Neural networks are skilled GPP simulators

An initial conclusion from our experiments is that neural network models have higher predictive skill than the theory-based P-model, across all levels of aggregation (Table 2). The main advantage of neural models is their capacity to represent complex functional dependencies, including effects that may not have been anticipated when deriving a model from plant physiological theory. Importantly, the neural networks predict GPP more accurately at unseen test sites. In other words, learning does not overfit the specific data streams at the training sites but discovers transferable patterns that are valid across space, and thus implicitly across environmental gradients. We attribute this robustness to the diversity of sites in our dataset, and to careful (fully automatic and data-driven) hyperparameter tuning.

In contrast to the deep learning models (LSTM+ and MLP+), the P-model implements rigid functional dependencies derived from a simplified depiction of the underlying processes (e.g., the big-leaf representation of canopy light absorption and a schematic, empirical treatment of water stress effects). The model has only few parameters to calibrate to the data (in our case four). On the one hand, so few degrees of freedom deprive the model of the ability to adapt to small but persistent effects present in the data; making it less accurate. On the other hand, they prevent it from going too far astray in the face of unexpected inputs.

Considering overall evaluations across all sites and dates, the two neural models perform equally well. However, the advantage of the LSTM+ over the other models was clear under certain conditions that we expected from the outset to underlie temporal structure in the data – water stress (Fig. 4) and frost/cold acclimation (Fig. 5). Under these conditions, the LSTM+ outperforms both the MLP+ and also the mechanistic model. This indicates a lack of current mechanistic understanding of processes affecting GPP under these conditions.

It should also be noted that the collection of sites used here represents only a limited subset of all relevant environments on Earth. None of the sites are located in a tropical ever-wet climate and certain conditions and combinations of vegetation types and species, environments, and plant growth conditions (soil, subsurface hydrology) may not be covered by our spatial cross-validation setup. The limited data availability is exacerbated by the fact that seasonal variations in GPP and environmental conditions tend to be very small in the tropics. Hence, it remains unclear whether the learned models extrapolate to conditions not well represented in our data collection (Ludwig et al., 2023; Meyer and Pebesma, 2022). While it is technically possible to predict GPP wherever the input predictors are available, such upscaling should be done with caution, and the limited reliability under particular conditions should be considered. Our results suggest that predictions are least reliable in regions with pronounced seasonal or perennial water limitation.

Analysing patterns in prediction error of the different models and model performance of out-of-sample predictions has revealed several key insights for (data-driven) modelling of terrestrial photosynthesis and its limitations. In the following, we discuss them in more detail.

4.3 Modelling cumulative and lagged effects on GPP benefits from a recurrent model

The prediction errors for the MLP+, LSTM_perm+ and P-model tended to grow with increasing water deficit, but errors of the LSTM+ remained smaller for moderate levels of water deficits compared to the other models (Fig. 4). Event analysis also showed that the recurrent architecture helps predict the response in GPP to drought events (Fig. 6). The benefit was observed primarily from 20 d after the start of a drought event. The similarity in prediction errors among models during the start of drought events and at lower levels of PCWD appears to indicate that the recurrent architecture only benefits the prediction of GPP at higher levels of water stress. On the other hand, as the LSTM+ was trained with sequences of 128 d, it could adapt to water limited conditions that build up over a period of up to 4 months. Due to the limited sequence length, longer periods of water stress were presumably not learned, which may contribute to rising errors towards the high end of cumulative water deficits and the relatively poor prediction of annual anomalies which may be driven by variable hydroclimatic conditions across years.

In this context, we point out that several factors likely degrade the prediction of annual anomalies: Inter-annual variability of ecosystem fluxes likely reflects effects by specific site histories which are not reflected in the predictor variables (Abramowitz et al., 2024), from inconsistencies in measurements of fluxes and meteorological covariates across years (e.g., sensor replacements), or by lagged effects of climatic extreme events (Zscheischler et al., 2014). By their nature, such effects are difficult to learn from example data spanning at most a few decades.

Adding soil moisture as a predictor to compensate for the MLP+'s lack of memory did not clearly boost overall performance. This indicates that (simulated) soil moisture does not fully account for the effects of gradually changing water stress, calling into question a widely used practice (Nelson et al., 2024; Kang et al., 2025; Tramontana et al., 2016; Gaber et al., 2024) and likely relates to a general challenge in accurately modelling water stress effects, which we discuss in more detail below. Soil moisture information did, however, improve GPP prediction at relatively arid sites at the cost of a slight drop at moist sites – nudging the behaviour of the MLP+ towards that of the LSTM+. We speculate that this trade-off could hint at a dependence of the functional relationships on aridity. The dependence of GPP on a soil moisture optimum that shifts in response to the growing season soil moisture (Peng et al., 2024) could also contribute to an advantage of the LSTM+ compared to the MLP+ with the current value of soil moisture.

Conversely, the LSTM+ had a (small) disadvantage in moist regions compared to the MLP+ (Fig. 9); while there was no obvious relation between the moisture index and the preference for global or site-specific modelling. Taken together, it seems that the LSTM+ more consistently generalises across different aridity levels than the non-recurrent model. This could be an indication that the functional relationships it uncovers hold over a wider range of aridity regimes. These patterns are unlikely to stem from uneven spatial representation, data sampling artifacts or differences in training data length. These were accounted for by ensuring even representation and sequence length of moist and arid sites (Fig. D1) as well as even distribution of moisture index values when splitting sites during training. Therefore, we attribute performance gaps to the ability of the models to account for the ecological complexity under arid conditions.

The LSTM+ also better captured delayed GPP increase in spring due to the cold acclimation effect (Fig. 5). Luo et al. (2023) found that a reduced efficiency of photosynthetic light utilisation during springtime was a consequence of a combination of low minimum temperatures and high radiation during the weeks and months leading up to and during the start of the growing season. A recurrent deep learning model offers a basis for more accurately modelling GPP under such conditions than non-recurrent architectures. Unresolved challenges remain, though, in the form of a remaining marked bias in the early part of spring also for the LSTM+.

Overall, we found that the LSTM+ benefits GPP prediction during periods that are affected by temporal patterns. For the temporal patterns evaluated, the LSTM_perm+ and MLP+ behaved similarly, while the LSTM+ showed a performance benefit. Furthermore, the LSTM_perm+ had the same model capacity as the LSTM+. Combined, this indicates that the improved prediction of LSTM+ is due to a learned memory, compared to the LSTM_perm+ and MLP+ which both only see contemporaneous predictor variables.

4.4 Unknown effects of water stress are a dominating source of model error

We found that GPP can relatively reliably be predicted across relatively moist, winter-cold sites. For sites with a moisture index of 0.75 and above, the mean R² for spatial out-of-sample predictions was 0.80. This indicates that – at least for the abiotic and biotic conditions represented in our dataset – GPP can be reliably simulated. Generalised models that spatially upscale yield relatively reliable results under such conditions with the R² of spatial out-of-sample tests falling between 0.4 and 0.94.

However, under more arid conditions (MI < 0.75) we found very variable performance of a generalised model. Different factors may cause poor generalisability across sites. Poor data quality with systematic differences of measurement errors across sites (Abramowitz et al., 2024), differences in functional relationships between GPP and its predictors across different species and vegetation types, or insufficient information in predictor variables all may underlie the variable performance of the global model across sites. Our results suggest that variable model performance is not clearly related to vegetation types (Fig. 7). A tendency of poorer model performance in evergreen vegetation is likely related to limited information in remotely sensed greenness which is provided as a predictor (fAPAR). A clearer relationship of model performance was found across climate zones and across the gradient of the site’s average aridity (Fig. 7). Together with our finding of a clear relationship between the model prediction error and potential cumulative water deficits (Fig. 4) and length of drought events (Fig. 6), this suggests that poor model generalisability is linked to variable exposure and response to water stress across sites. Apparently, the history of precipitation and radiation is not sufficient to accurately model vegetation water stress exposure and responses, and effects on GPP.

Two factors are likely to undermine generalisability. First, responses to declining water potentials in the rooting zone are highly variable across species and linked to plant hydraulic traits and water use strategies. Even within broad classes of vegetation types, hydraulic relations of different plant species exhibit a wide variety (Choat et al., 2012; Joshi et al., 2022; Anderegg et al., 2018; Xu et al., 2016; Whitley et al., 2017; Konings and Gentine, 2017). Particularly in dry-adapted ecosystems (e.g., savannas and shrublands), contributions of different species to ecosystem-level integrated fluxes may also change over the season as a result of species-specific responses of leaf area to dryness (Xu et al., 2016; Whitley et al., 2017). Without related information provided to models, this complexity and the resulting variability of GPP responses to dryness cannot accurately be modelled across ecosystems with different species compositions.

Second, the exposure to water stress is highly variable across sites as a result of the surrounding topography and subsurface hydrology. Giardina et al. (2023) found strong variations of the functional relationship between evapotranspiration and cumulative water deficits, suggesting strongly variable rooting zone water storage capacities and plant access to groundwater across sites (Fan et al., 2017). Subsurface hydrology, groundwater influence, and belowground moisture convergence also appear to lead to large differences in ecosystem water balances at relatively dry sites (Hahm et al., 2019; McCormick et al., 2021). Several flux measurement sites have been identified as having greater mean annual precipitation than evapotranspiration, suggesting subsurface moisture convergence and the influence of groundwater (Abramowitz et al., 2024). Due to the close link between evapotranspiration and GPP, these relations affect vegetation activity in general, including GPP. The surrounding topography, water holding capacity of the soil and weathered bedrock, groundwater table depth and rooting depth are either insufficiently known or specified from the predictors used for modelling GPP. Thus, related effects on GPP and the associated variability of water stress exposure cannot accurately be modelled across space. Variable water stress exposure and response affect GPP in relatively moist and energy-limited sites to a lesser degree than sites with frequent water limitation. Hence, a clear relation of model generalisability across aridity (Fig. 8) emerges.

Another factor that contributes to these generalisation challenges is remotely sensed fAPAR, which is a predictor of GPP and tends to be less accurate in arid regions, where satellite sensors struggle to capture the large spatial heterogeneity at sub-pixel scales and may be affected by light absorption by non-photosynthetically active tissue (Kannenberg et al., 2024). In particular, noise due to sparsity, vegetation senescence and soil background lead to a frequent overestimation of fAPAR in arid and semiarid regions (Smith et al., 2019). This is exacerbated by the relative lack of ground observations in drylands, needed for calibration.

We expected that site-specific responses would be more effectively modelled by site-specific models compared to a generalised, global model. However, this was not unanimously the case. While the mean prediction error (RMSE) was generally lower for site-specific models, these models often predicted a smaller fraction of variation in the data than the global models. This was most clearly found for sites for which relatively short time series were available for model training. Hence, site-specific responses appear to be learnable, given sufficient data. Our interpretation is that R² measures the ability to explain the variance in GPP, a task that becomes easier as the model sees more data. Put differently, a site-specific model can more accurately memorize the mean seasonal cycle of one particular site; but may not learn as well to deduce daily variations from observed changes in light and meteorology, due to its restricted sample. Indeed, sites where the global model was better (positive ΔR²) invariably had relatively short observation periods (Fig. 8), which increases the need to learn parts of the functional relation from other sites.

While the data-driven models showed advantages in certain situations with temporal effects, it is clear that challenges remain. Even with memory and accessible precipitation history, exposure and responses to water stress is not sufficiently modelled. This applies to models of GPP in general, including mechanistic models. As outlined above, it is likely that models need additional information to simulate GPP reliably at arid sites. Our experiments indicated that addition of commonly used satellite-derived variables boosts overall model performance without resolving the large variability in model performance across dry sites. This highlights the need for the inclusion of further information that is relevant for accurate drought response simulations. We expect that a similar reasoning and potential for model improvement applies also for mechanistic land surface models. The integration of additional satellite data to inform models for reliable simulations of water stress exposure in land surface models remains a largely unsolved task. Typically, land surface models account for water limitation by predicting leaf area index and stomatal conductance internally. Therefore, unknown exposure to water stress will likely remain a large source of uncertainty in land surface models. Similarly, mechanistic models like the P-model use fAPAR, but no other remote sensing data due to challenges in linking them in process representations.

5 Conclusions

We have demonstrated that an LSTM – a popular type of recurrent deep neural network – is a powerful model to predict ecosystem GPP from local meteorological observations and remotely sensed fAPAR. Based on a spatial and temporal out-of-sample evaluation, we find that the model has significantly higher predictive skill than the theory-based P-model and outperforms a non-recurrent deep learning model under conditions of low root-zone moisture availability and very low temperatures in preceding weeks. The LSTM reliably simulates GPP dynamics across a range of environmental conditions and vegetation types (no agricultural vegetation tested here) at relatively moist sites (MI > 0.75).

Through a detailed analysis of error patterns, we find that a recurrent model more accurately captures the GPP response to longer-term, cumulative impact. In particular, the LSTM adapts better to arid environments affected by water stress, a condition that builds up over time. Yet, we find that there is still a large variability in model skill across relatively arid sites, even if it outperforms both the mechanistic P-model and a memoryless neural network. The variability remains even with the inclusion of additional earth observation data, although this improves general model performance. This suggests that the model lacks information on variations in exposure and response to water stress and related effects on GPP. The inclusion of additional remotely-sensed and temporally varying information (e.g., sun-induced fluorescence (Li et al., 2018), vegetation optical depth (Konings and Gentine, 2017)) or static information about the topography, average groundwater table depth (Fan et al., 2013), and subsurface structure (Pelletier et al., 2016) as predictors for the deep learning models bears the potential for reducing errors and yielding more reliable GPP simulations in dry environments. As ecosystems are becoming more exposed to water limitation due to climate change (Denissen et al., 2022; Fu et al., 2024), it remains an important research topic to improve the predictability of ecosystem fluxes in the context of water stress.

Appendix A: Site information

Table A1 Overview of sites included in the study. “MI” is the moisture index. “Clim.” is the Köppen-Geiger climate zone. “Veg.” is the IGBP vegetation type. “Evergeen” is the classification of sites into categories used in Fig. 4. “Delayed GPP” is the classification of sites into categories used in Fig. 5.

Appendix B: Methods details

Beck et al. (2018)

Table B1 Overview of the data sources.

Table B2 Overview of the implemented models.

Table B3 Cross-validation and training setups.

Appendix C: Sequence models

We evaluated two additional sequence models using the global model setup.

• Temporal Convolutional Network. We implemented a Temporal Convolutional Network (TCN) following Bai et al. (2018). Similar to the LSTM, we used GELU activations (Hendrycks and Gimpel, 2023) and LayerNorm (Ba et al., 2016).

• LSTM + attention. We also implemented a second version of the LSTM that uses a multi-head (8) attention mechanism (Vaswani et al., 2017) with two linear layers after the LSTM layer.

Figure C1 Comparison of the LSTM against the TCN (a) and against the LSTM + attention (b). R² of predicted versus observed values of daily GPP per site from the global cross-validation is shown for the LSTM along the y-axis and the TCN/LSTM + attention along the x-axis. The dotted line indicates equal performance.

Figure C2 The model bias per day before and after the start of drought events for three sequence models. The bias is calculated as the difference between model predictions and observations. The bias is normalised with respect to the median value during the 1 to 10 d before each event onset. Shaded areas mark the area between the 33rd and 66th percentiles.

Appendix D: Additional figures

Figure D1 Number of sites as well as sequence lengths per site for different levels of the moisture index. The number of sites are counted per bin of width 0.2 of the moisture index.

Figure D2 The distribution of performance metrics for sites with a moisture index below 0.75, and sites with a moisture index above 0.75. The individual values as well as box plots are shown.