2.1 CMIP6 model data

We compare vegetation productivity responses in models from the CMIP6 esm-hist experiment (Eyring et al., 2016). This experiment covers the recent historical period from 1850 to 2014. Only data from 2000–2014 are included in our analysis in order to permit a fair comparison with the available satellite observations. To investigate vegetation productivity responses to intraseasonal rainfall variability, we require daily or sub-daily model output data for precipitation and GPP. These data are available for the esm-hist experiment for only five models: ACCESS-ESM1-5, BCC-CSM2-MR, CNRM-ESM2-1, NorESM2-LM, and UKESM1-0-LL. We use the precipitation fields provided at daily resolution and take the daily mean of the GPP data, which is made available at 3 h resolution. These models also all provide daily surface soil moisture (mrsos) data: this is the mass of water in the upper 10 cm of soil. We compute daily near-surface VPD from near-surface air temperature (tas) and relative humidity (hurs). This is done by first computing the saturation vapour pressure e_s according to Tetens' formula (Bolton, 1980):

where tas is expressed in degrees Celsius (°C). The vapour pressure deficit is then

Near-surface relative humidity is not available at daily resolution for BCC-CSM2-MR, so this model is excluded from the VPD-focused sections of our analysis. All model data are regridded to 1°×1° horizontal resolution using land-area-weighted averaging.

2.2 Observational datasets

Intraseasonal rainfall events are identified using the Integrated Multi-satellitE Retrievals for GPM (IMERG) V06 daily product (Huffman et al., 2019). Surface soil moisture is assessed using both the ESA CCI Soil Moisture combined active–passive microwave product v06.1 (Dorigo et al., 2017; Gruber et al., 2019) and the Global Land Evaporation Amsterdam Model (GLEAM) v3.6a (Miralles et al., 2011; Martens et al., 2017). GLEAM models the distribution of soil water content based on Multi-Source Weighted-Ensemble Precipitation (MSWEP) v2.8 (Beck et al., 2017) and corrects the soil moisture of the top model layer by assimilating the ESA CCI Soil Moisture combined product. The two surface soil moisture datasets are therefore not independent, but both are tested here to ascertain whether the differences in the products affect the assessment of the CMIP6 models' surface soil moisture. ESA CCI Soil Moisture has the advantage of being less reliant on model algorithms, while the surface soil layer in GLEAM is 10 cm deep and may therefore be more representative of the CMIP mrsos variable, which is defined over the same depth. The microwave observations used in ESA CCI Soil Moisture measure varying soil depths but are typically taken to quantify soil moisture in the top 2–5 cm (Ulaby et al., 1982). We compute near-surface VPD from ERA5 reanalysis data for 2 m air temperature and 2 m dew-point temperature (Hersbach et al., 2020).

Daily GPP data are obtained from FLUXCOM RS+METEO (Jung et al., 2020). FLUXCOM RS+METEO uses machine learning to upscale estimates of terrestrial carbon fluxes from eddy covariance flux towers to create a global gridded GPP product based on remote sensing and meteorological forcing data. We test versions of FLUXCOM RS+METEO using two different meteorological forcing datasets: ERA5 and CRU JRA v1.1 (Harris et al., 2014; Kobayashi et al., 2015; Harris, 2019). For each of these forcing datasets, we use the ensemble median GPP over three machine learning methods and two flux partitioning methods, as detailed by Tramontana et al. (2016) and Jung et al. (2019). Note that while the meteorological forcing data are updated daily, the remote sensing driving data are based on mean seasonal cycles (Jung et al., 2020).

Whilst daily data are preferable for investigating intraseasonal variability, information can also be obtained from other sub-monthly GPP datasets. Incorporating these datasets into the analysis allows a better understanding of the uncertainty in our results that arises from the need to derive GPP from direct observations. We include the 8-daily MODIS Terra GPP product (Running et al., 2015), as quality-controlled and regridded by Kern (2021). We also analyse the 8-daily VPM GPP product (Zhang et al., 2017a), which is intended to provide an alternative to MODIS Terra GPP by using the improved light use efficiency theory of the vegetation photosynthesis model (VPM). The VODCA2GPP dataset (Wild et al., 2022), which estimates GPP based on vegetation optical depth (VOD) retrieved from passive microwave observations, is also produced at 8 d resolution.

Satellite observations of solar-induced chlorophyll fluorescence (SIF) provide an additional way of estimating GPP at the global scale. SIF is approximately linearly correlated with GPP (Frankenberg et al., 2011), although the relationship varies with season and vegetation type (Chen et al., 2021). We therefore also compare the modelled GPP responses with SIF data that are spatially downscaled from GOME-2 observations (Duveiller et al., 2020). This product includes data based on two different methods of retrieval from GOME-2: Joiner et al. (2013), henceforth labelled JJ, and Köhler et al. (2015), labelled PK. These datasets have a temporal resolution of 8 d using a 16 d rolling window.

All observational products are available for the complete period 2000–2014, except for IMERG, which is available from June 2000 onwards, and the GOME-2 SIF data, which are available from 2007 onwards. All products are regridded by land-area-weighted averaging to a horizontal resolution of 1°×1°.

2.3 Compositing around intraseasonal precipitation events

Intraseasonal precipitation events are identified using the method of Harris et al. (2022). Long-term linear temporal trends in the data are removed before processing. For each 1°×1° grid box, we apply a 25 d low-pass Lanczos filter to the daily precipitation anomaly, where the anomaly is computed relative to the climatology for a 7 d rolling window. This removes variability at frequencies higher than 25 d from the anomaly time series. Intraseasonal precipitation events are then defined as local maxima of the filtered time series that lie above 1 standard deviation from the mean. This method identifies wetter-than-usual intraseasonal periods with a characteristic timescale of approximately 25 d. The elevated precipitation associated with an intraseasonal wet event may be a single large rainfall event or it may comprise several shorter bursts of rainfall. Using such a filtering approach ensures that the driving precipitation composites have a consistent timescale so that they provide a fair basis on which to compare the timescales of the GPP responses between the models and the observations. The dates of these events are identified separately for each CMIP6 model and for IMERG. Due to gaps in the IMERG dataset at high latitudes, this procedure results in far fewer observed events than modelled events outside the latitude band 60° S–60° N. We therefore restrict all analyses of both modelled and observational data to 60° S–60° N in order to ensure a fair comparison between the two. We then composite daily standardised anomalies (relative to 7 d rolling climatology) of precipitation, SSM, VPD, and GPP around the dates of the intraseasonal precipitation events. For the observational GPP products with 8 d resolution, the standardised anomalies are computed relative to a 31 d rolling climatology instead of 7 d to ensure sufficient data for the climatological averaging. It should be noted that the composited variables do not undergo Lanczos filtering; the filter is only used to determine the dates of the precipitation events around which to composite.

The results are aggregated over regions containing multiple grid boxes to increase the number of precipitation events contributing to each composite. Only grid boxes with valid ESA CCI Soil Moisture observations are included to maximise our observational knowledge of the surface soil moisture perturbations that are driving the GPP responses. Although this eliminates tropical forest regions from the analysis, it is found to make only a small difference to the resulting GPP composites, since vegetation in these regions does not respond strongly to intraseasonal wet events compared to vegetation in water-limited regions (Harris et al., 2022). We also remove frozen grid boxes from the composites by discarding events in months when the median of maximum 2 m air temperature is below 0 °C.

3.1 Global evaluation

Global (60° S–60° N) composites of precipitation, surface soil moisture, and GPP around intraseasonal wet events are shown in Fig. 1. The composites of precipitation are similar between the CMIP6 models and observations (Fig. 1a), indicating that our method successfully creates similar standardised anomalies of precipitation driving the land surface responses in each case. The surface soil moisture responses following these wet events are also similar between models (Fig. 1b), with consistent maximum standardised anomalies 2–3 d after the wet event peak and comparable longer-term anomalies of elevated soil moisture out to 60 d. The observed surface soil moisture composites from the GLEAM and ESA CCI products are also a reasonable match to the CMIP6 models. The models are in better agreement with the observed composite obtained from GLEAM, which may be due to GLEAM representing the same depth of soil surface layer as the model output.

Figure 1Composites of (a) precipitation, (b) surface soil moisture, and (c, d) GPP around intraseasonal wet events. All land pixels between 60° S and 60° N are included, except in months when the ground is frozen or in regions without valid ESA CCI Soil Moisture observations. The surface soil moisture and GPP composites have been scaled by the maximum amplitude of their corresponding precipitation composite.

Download

However, the responses of GPP to the intraseasonal wet events show large inter-model differences (Fig. 1c). All models exhibit a positive GPP anomaly in the days following the peak of wet events, but the amplitude and timing of this anomaly vary greatly. CNRM-ESM and UKESM produce a larger standardised anomaly in GPP than the other models, while ACCESS-ESM and NorESM show the smallest responses. NorESM also has the slowest response in GPP, with the peak GPP anomaly occurring 22 d after the peak precipitation anomaly compared to 4 d for BCC-CSM, the fastest-responding model. The models also show contrasting behaviour leading up to the peak of the wet events: some have a much larger negative GPP anomaly than others, associated with a reduction in downwelling shortwave radiation during the wet events. However, in this study we focus on the post-event characteristics of the GPP response.

To investigate which of the modelled productivity responses are realistic, the wet event GPP composites are shown in Fig. 1d for a variety of datasets based on observations (as described in Sect. 2). While composites obtained from the same family of observational products (e.g. the two SIF retrieval methods) are similar to one another, there is a large spread in the responses between the types of observational products, with anomaly magnitudes and timings covering most of the behaviour seen in the models. The exception to this is the faster response of BCC-CSM, for which we do not find any observational evidence.

Whilst a horizontal resolution of 1° is relevant for evaluating ESMs, which represent an average coarse-scale behaviour of the Earth system, it is possible that the intraseasonal precipitation events and the vegetation responses do not spatially co-occur within this 1° grid box in observations. Therefore, to test whether the 1° GPP responses are attributable to the precipitation anomalies, we also composite the GPP response around intraseasonal wet events using the IMERG and VPM datasets at 0.25 and 0.1° horizontal resolution. The global composite GPP response is not found to be sensitive to spatial resolution (Fig. S1 in the Supplement). This indicates that the coarser 1° response, which represents the vegetation response at a scale relevant to ESMs, is representative of the processes occurring at much smaller spatial scales in the real world.

It is therefore clear that there is great uncertainty in the response of GPP to intraseasonal rainfall events in CMIP6 Earth system models but that the observations also fail to constrain this response. This makes it challenging to assess which of the model responses are most realistic from a global perspective. In order to further investigate whether the models are providing appropriate GPP responses, we therefore move to a regional evaluation.

3.2 Regional evaluation

To provide a regional assessment of GPP responses to intraseasonal rainfall events, we create composites for each IPCC AR6 land region (Iturbide et al., 2020) using the same method as in Sect. 3.1. This gives 46 regions, which are designed to represent areas with consistent climate features (Iturbide et al., 2020). Only 40 regions are included in this analysis: GIC (Greenland and Iceland), RAR (Russian Arctic), and the two Antarctic regions lie entirely outside the 60° S–60° N latitude band in which we identify intraseasonal precipitation events. We also remove the desert regions SAH (Sahara) and ARP (Arabian Peninsula) because these show outlying relationships between modelled and observed GPP anomalies, which is likely to be due to the very low GPP in these regions. All regions included in the analysis have a minimum of 500 precipitation events identified in each model and in observations. Full details of the sample sizes of events are included in Fig. S2.

In order to summarise the GPP responses in each region, we focus on three key properties of the composites, which are illustrated in Fig. 2a. The peak amplitude in GPP is defined as the maximum standardised anomaly occurring after the peak of the rainfall event composite (i.e. after day 0). The lag is the number of days after the rainfall peak that this peak amplitude occurs. The post-event amplitude is the mean standardised anomaly over days 40–60 after the rainfall peak, a period when there is no longer an anomaly in precipitation. The amplitude metrics characterise the contribution of intraseasonal wet events to GPP variability, while the lag measures how quickly GPP responds to a wet event.

Figure 2(a) Illustration of the three properties of the GPP response composite that will be analysed. The example composite used for the illustration is the GPP response from UKESM in the northern Australia region. (b–d) Comparison of the peak amplitude, post-event amplitude, and lag of the GPP response to intraseasonal wet events between the CMIP6 models and the FLUXCOM-CRUJRAv1 observational product. Each scatter point represents an AR6 region. Dashed lines show linear best fits and panel legends detail the linear correlation coefficients. Asterisks after correlation coefficients denote significance at the 95 % level. The grey lines indicate the identity line in each panel, i.e. where agreement between modelled and observed values would occur. The peak and post-event amplitudes in GPP have been scaled by the peak and post-event amplitudes in surface soil moisture, respectively.

Download

Figure 2 compares the modelled peak amplitude, lag, and post-event amplitude of GPP in each AR6 region to values from one of the observation-based products. The observational product used for this comparison is FLUXCOM RS+METEO driven by CRU JRA v1 reanalysis (henceforth referred to as FLUXCOM-CRUJRAv1). This product is chosen as a starting point because FLUXCOM RS+METEO provides the only available daily GPP data, but other observational products will be analysed later in this section. The peak and post-event amplitudes of GPP are scaled according to the observed peak and post-event amplitudes in surface soil moisture, respectively (i.e. the modelled peak amplitude in GPP is multiplied by the ratio between the observed and modelled peak amplitudes in surface soil moisture). Surface soil moisture is chosen as the scaling variable because the GPP amplitude is strongly linearly related to the surface soil moisture amplitude (see Fig. 4 and accompanying analysis later in the section). ESA CCI Soil Moisture is used as the benchmark observation. This scaling accounts for the possibility that individual models may under- or overestimate the regional surface soil moisture perturbation at the regional scale following wet events. Whilst a larger surface soil moisture perturbation is expected to lead to a larger-amplitude response in GPP, the scaled data allow us to assess the sensitivity of the GPP response to surface soil moisture by adjusting to identical surface soil moisture perturbations for each model.

All the analysed CMIP6 models show positive correlation for the peak and post-event amplitudes of regional GPP responses with FLUXCOM-CRUJRAv1, significant at the 95 % level (Fig. 2b and c, respectively). Importantly, this indicates that the models are generally able to correctly represent which regions develop larger anomalies in GPP following intraseasonal rainfall events. However, it is clear that the actual modelled values of response amplitude are very different between models. For example, in regions that show higher peak amplitudes in FLUXCOM-CRUJRAv1, UKESM is able to represent these higher amplitudes, whereas the amplitudes modelled by ACCESS-ESM and NorESM are much lower.

The lag of modelled regional GPP responses (Fig. 2d) is significantly correlated with observed lags for ACCESS-ESM, BCC-CSM, CNRM-ESM, and UKESM. No significant linear relationship is found between observed and modelled lags for NorESM, indicating that the timing of GPP responses to rainfall may not be realistic in this model. In many regions, the GPP lag modelled by NorESM is much longer (20–60 d) than the observed lag in any region. This suggests that the NorESM GPP response is generally too slow, which is consistent with this model showing the slowest response in the global composites in Fig. 1. Even for the models where a significant correlation in lag is found, the correlation coefficients are much lower than for the amplitude relationships, showing that the models match the observations better for the magnitude of the GPP response than the timing.

Given the uncertainty in observed GPP responses established in Sect. 3.1, it is important to assess the effect of changing the observational GPP product on the results of the assessment carried out in Fig. 2. The Taylor diagrams presented in Fig. 3 show the correlation between regional responses for each possible model–observation pair, in addition to the inter-regional standard deviation normalised by the observed standard deviation $\left(\frac{{\mathit{\sigma }}_{\text{model}}}{{\mathit{\sigma }}_{\text{obs}}}\right)$ and the root mean squared error normalised by the observed standard deviation. For example, in Fig. 3a, the orange cross summarises the linear relationship between UKESM and FLUXCOM-CRUJRAv1 that is shown in full in Fig. 2b, which has a correlation of 0.89 and a slightly higher spread in peak GPP amplitudes between regions for the model than the observations $\left(\text{i.e.}\frac{{\mathit{\sigma }}_{\text{model}}}{{\mathit{\sigma }}_{\text{obs}}}\approx \mathrm{1.2}\right)$.

Figure 3Taylor diagrams of regional relationships between modelled and observed GPP responses. Markers show the correlation and ratio of standard deviations $\frac{{\mathit{\sigma }}_{\text{model}}}{{\mathit{\sigma }}_{\text{obs}}}$ between the model and each observation-based GPP product for the (a) peak amplitude, (b) post-event amplitude, and (c) lag of the GPP response to intraseasonal wet events. The colour of the marker indicates which model is being evaluated, and the shape of the marker indicates the observation-based product used. The labelled grey contours indicate the root mean squared error normalised by observed standard deviation. GPP peak and post-event amplitudes have been scaled by the modelled peak and post-event amplitudes in surface soil moisture, respectively.

Download

When focusing on the peak amplitude of the GPP response, the correlation between the model and observations shows some variation based on the choice of observational product, but the value typically lies between 0.8 and 0.9. The exception is NorESM, which is consistently a worse match for the observed regional distribution of GPP response amplitudes (r≈0.6), regardless of the observational product selected. The extent to which the inter-regional standard deviation in peak amplitudes matches the observed standard deviation is much more dependent on the choice of model than on the choice of observational product. This reflects the fact that some models produce much higher maximum regional peak amplitudes than others (as seen in Fig. 2). NorESM produces the lowest inter-regional variation in GPP response amplitudes, related to its inability to represent the larger responses; σ_model is far too low compared to the observed standard deviation, regardless of product. ACCESS-ESM also underestimates the inter-regional variation in response amplitudes compared to the observational products, while UKESM overestimates it. CNRM-ESM and BCC-CSM both reasonably capture the variation, falling within the spread of the observational products. Combined with their high correlation coefficients (r>0.8), these two models therefore provide the most realistic representation of the differences between regions in the peak amplitude of GPP responses to intraseasonal rainfall events.

The correlation between models and observations is more dependent on the choice of observational product for the post-event amplitude (Fig. 3b – note the change in the radial axis scale from Fig. 3a) than for the peak amplitude. This indicates inter-product disagreement on which regions experience elevated GPP in the months following wet events. CNRM-ESM outperforms the other models here in terms of correlation and also produces a standard deviation consistent with several of the observational datasets. Similarly to the results for the peak amplitude, NorESM and ACCESS-ESM tend to underestimate the inter-regional variation in response, while UKESM overestimates it. In other words, intraseasonal rainfall events in NorESM and ACCESS-ESM play too minor a role in GPP variability compared to in observations.

Whereas all pairs of models and observations show correlations with r>0.5 for both the peak and post-event amplitudes, the modelled lags in the regional GPP responses show much less agreement with the observed lags (Fig. 3c). Many of the relationships do not show a correlation significant at the 95 % level, with some even having a negative correlation. For each model, there are also large spreads in both the correlation and $\frac{{\mathit{\sigma }}_{\text{model}}}{{\mathit{\sigma }}_{\text{obs}}}$ dependent on the observations used for the assessment. The timing of GPP responses to intraseasonal rainfall events is much less well constrained by observations than the amplitude.

We now consider the possible causes of such large differences between models. Different inter-model GPP responses may arise due to differing land cover maps (e.g. if forests show a weaker response to rainfall events than other land cover, then a model with more extensive forest cover will produce a weaker globally composited GPP response). Alternatively, differences in the GPP response could be due to variations in the models' representation of processes linking water availability and vegetation productivity so that changes in GPP differ between models even when identical land cover is being considered.

To assess whether the first possibility is solely responsible for the differing GPP responses seen in Fig. 1c, we compare regional land cover fractions between the models (Appendix A). We find that inter-model differences in the amplitude and lag of GPP responses can be seen even in regions where the models have similar land cover. Therefore, land cover differences cannot by themselves explain the discrepancies in the GPP response, so we further investigate the differences in water–carbon coupling processes. Two important processes that have known uncertainty in their representation in ESMs are the control of soil moisture stress on GPP and the response of stomatal conductance to changes in vapour pressure deficit (VPD). We therefore investigate the regional relationships between the GPP response and the perturbations in surface soil moisture (SSM) and near-surface VPD following the wet events.

Figure 4Regional GPP responses to intraseasonal wet events compared to the driving perturbations in surface soil moisture (SSM: a, c, e) and vapour pressure deficit (VPD: b, d, f). Each scatter point represents an IPCC AR6 region. Dashed lines show linear best fits. The legend for each panel indicates the linear correlation coefficient between the driving perturbation and the GPP response for each model. Asterisks after correlation coefficients denote significance at the 95 % level. The uppermost legend, showing which model is denoted by each colour, applies to all panels. The observational products (OBS) used are FLUXCOM-CRUJRAv1 (GPP), ESA CCI (surface soil moisture), and ERA5 (VPD).

Download

Figure 4 explicitly compares these relationships in each model and in observations. Again, FLUXCOM-CRUJRAv1 is used as an illustrative observational GPP product, but the other products will be included later in the section to test the robustness of the relationships. ESA CCI Soil Moisture is used for the SSM observations, but the results are qualitatively unchanged if it is substituted with GLEAM (not shown). When defining the peak amplitude of the VPD response, the minimum rather than the maximum standardised anomaly is taken, since VPD typically decreases during and after intraseasonal wet events. Note that daily VPD data were not available for BCC-CSM, so this model is not included in the VPD analysis. SSM and VPD perturbations are not independent from one another; the regional amplitudes of SSM and VPD responses are strongly negatively correlated (not shown). However, we test the relationship of GPP with both variables since they influence the GPP through different modelled processes.

All models represent the observed positive correlation between the SSM peak amplitude and the GPP peak amplitude: regions with a larger standardised anomaly in SSM following rainfall events show a larger standardised anomaly in GPP due to the greater change in water availability. The models also capture the observed negative correlation between the peak amplitudes of the VPD and GPP responses. Regions with a smaller drop in VPD standardised anomaly after wet events are experiencing less evaporative demand compared to climatological conditions, so stomata can remain open to allow photosynthesis without experiencing as much water loss to the atmosphere. These correlations are significant (p<0.05) in all the observational products tested (Fig. S4a and b). However, even though all models represent these basic relationships, there are large differences between the correlation coefficients, meaning that GPP responses to rainfall are more tightly coupled to the associated perturbations in SSM and VPD in some models than others. For example, in CNRM-ESM2-1, the correlation between the peak amplitude of the SSM perturbations and the peak amplitude of the GPP response is r=0.94, meaning that the magnitude of the SSM perturbations can explain 89 % of the variance in the GPP peak amplitude between regions (r²=0.89), whereas in UKESM1-0-LL it only explains 52 % of the variance (r=0.72, r²=0.52).

There are also differences between models in the slopes of the linear relationships. This demonstrates that the models' GPP peak amplitudes have varying sensitivities to changes in the peak amplitude of the SSM and VPD perturbations. For example, for a given increase in the peak standardised anomaly of SSM, NorESM exhibits a much smaller increase in the associated GPP peak anomaly than the other models. We test the significance of these differences in sensitivity among all the models and all the observational GPP products. This is done using estimated marginal means of linear trends and correcting for multiple testing using the Tukey method. Figure 5 shows the resulting 95 % confidence intervals for the sensitivity of the GPP responses to the driving perturbations in SSM and VPD. This analysis shows that for the relationship between the peak amplitudes of VPD and GPP (Fig. 5b), the sensitivity of NorESM is significantly different (p<0.1) to the sensitivity of all other models and all observational GPP products except SIF-GOME2-JJ. These differences are also all significant at p<0.05 other than NorESM-VODCA2GPP. Therefore, we conclude that the amplitude of GPP responses to wet events in NorESM is not sufficiently sensitive to the changes in vapour pressure deficit.

Figure 5Sensitivity of regional GPP responses to intraseasonal wet events compared to the driving perturbations in surface soil moisture (a, c, e) and VPD (b, d, f). Error bars denote 95 % confidence intervals. The sensitivity of the observational GPP products is measured with respect to ESA CCI for surface soil moisture and ERA5 for VPD.

Download

The peak amplitude of GPP responses is also less sensitive to SSM perturbations in NorESM than in the other models (Fig. 5a), showing a significantly (p<0.05) different sensitivity to BCC-CSM, CNRM-ESM, and UKESM. However, of the seven observational products, the SSM sensitivity in NorESM is only significantly different to FLUXCOM-ERA5, so there is insufficient evidence here to state that it is unrealistic compared to the range of observations. In contrast, the peak amplitude of GPP responses in BCC-CSM shows excessive sensitivity to SSM perturbations, with significant (p<0.05) differences to MODIS Terra, VPM, SIF-GOME2-JJ, and VODCA2GPP (and p<0.1 for the difference with FLUXCOM-CRUJRAv1). It additionally shows a significantly higher sensitivity than two of the other models (NorESM and ACCESS-ESM).

Similar relationships between SSM and VPD as well as GPP are also found for the post-event amplitudes, but the correlation is generally weaker than for the peak amplitudes. The peak amplitude of GPP is more strongly constrained by the concurrent SSM and VPD anomalies than the post-event amplitude. In the case of NorESM, the correlation is sufficiently weak that no significant relationship is found between SSM and GPP amplitudes after the event when testing at the 95 % level; vegetation productivity is much less strongly linked to surface soil water conditions than in the other models and in the observations (again, the correlations are significant at p<0.05 in all the observational products tested). The weaker relationships between post-event amplitudes can be partly explained by the fact that over the 40–60 d timescale, changes in root zone soil moisture explain more of the variation in the GPP amplitudes than changes only in the surface layer (whereas the surface layer is generally more important for the peak amplitude). Comparing the post-event amplitude of soil moisture in the top 1 m of soil with the post-event amplitude of GPP gives higher correlation coefficients than seen in Fig. 4c for all models (other than BCC-CSM, for which the 1 m layer data are not available), as shown in Fig. S5. In the months after the peak of the wet event, the rainfall infiltrates the soil and provides additional moisture availability to vegetation at deeper soil levels. Figures 4c and 5c show that NorESM is less sensitive to the SSM post-event amplitude than the other models, although its sensitivity is significantly (p<0.05) different only to UKESM, CNRM-ESM, and four out of the seven observational products. The modelled sensitivity to VPD post-event amplitude is not inconsistent with observations for any model.

The lag in GPP response is compared to the peak amplitude rather than the lag of SSM and VPD because there is very little variation in SSM and VPD lags between regions due to these variables responding very quickly (1–3 d) after rainfall. In the models, the lag of the GPP response tends to be shorter in regions with larger perturbations in SSM and VPD. All models exhibit this relationship except NorESM. However, the evidence for the relationship in the observations is much weaker (see Fig. S4). Only two of the seven observational products tested show a significant (p<0.05) negative correlation between SSM peak amplitude and GPP lag (FLUXCOM-ERA5 and VODCA2GPP), while VPM shows a significant positive correlation. Four out of the seven products show a significant positive correlation between VPD peak amplitude and GPP lag (FLUXCOM-ERA5, FLUXCOM-CRUJRAv1, SIF-GOME2-PK, and VODCA2GPP). Even in the observational products that do show these significant relationships, the correlation coefficients are much lower than those obtained from the models. The highest correlation coefficients from observations are from VODCA2GPP at $r=-\mathrm{0.50}$ for the correlation between SSM peak amplitude and GPP lag and r=0.56 for VPD peak amplitude and GPP lag. These correlations are considerably weaker than those found in the models (other than NorESM), which have $\left|r\right|$ between 0.58 and 0.73. It is therefore possible that the timing of the GPP responses in the models is too tightly constrained by the amplitude of the regional perturbations in SSM and VPD.