Half-hourly eddy covariance (EC) flux data were used in this study to estimate the Amax of global grasslands. Data were retrieved from the FLUXNET 2015 Tier 1 dataset. ¹ Intotal, 15 EC flux sites were analyzed in this study. These sites were selected based on the data availability (i.e., sites with more than 60% of Amax_a data over the available years). Meteorological data measured alongside the flux data were also used, including soil water content (SWC), vapor pressure deficit (VPD), air temperature (Ta), and solar shortwave radiation (Rad). Detailed information for the sites used in this study is listed in Table 1.

Leaf area index data (LAI) were extracted from the moderate resolution imaging spectroradiometer (MODIS) Leaf Area Index/FPAR Products version 6.1 (MOD15A2H v061). This product was an 8-day composite dataset with a spatial resolution of 500 m. We filtered out those data points where significant clouds or snow were present, or where the data quality was poor according to the quality flags (FparLai_QC and FparExtra_QC) provided by this dataset. The average LAI of pixels within a 500 m radius around the flux tower was used to represent the LAI of the site, because Chu et al., 2021 suggested that target areas with a smaller radial distance from the flux tower tend to be less biased in representing the land surface characteristics as covered by the flux footprint climatology. Then, the LAI time series were smoothed using the Savitzky–Golay (SG) filter, which has widely been used to remove the noises in MODIS LAI data (Zhu et al., 2016; Huang et al., 2021).

Normalized difference vegetation index (NDVI) data were downloaded from the MODIS C6 (MOD13C1) which had a resolution of 16-day and 500 m. The average NDVI of pixels that are within a 500 m radius around the flux tower was used to represent the NDVI of the site. To match the temporal resolution of LAI, we interpolated the NDVI time series to the LAI dates using linear interpolation.

We used the daytime flux partitioning method to derive the parameters of photosynthesis and respiration. The daytime flux partitioning method uses a light response curve to fit the flux of CO₂. We followed the method proposed by Lasslop et al. (2010), with a small modification of separately estimating day-time and night-time reference respiration rates to account for light inhibition of leaf respiration (Keenan et al., 2019):

where α (μmol C J⁻¹) is the canopy-scale quantum yields; β (μmol C m⁻² s⁻¹) is canopy light-saturated CO₂ assimilation rate; R_g (W m⁻²) is incoming shortwave radiation; and F_r (μmol C m⁻² s⁻¹) is the ecosystem respiration term. Parameter β is estimated as an exponentially decreasing function of VPD, to account for the effect of VPD on apparent photosynthesis:

where β₀ and k are fitted parameters and VPD₀ is set to 10 hPa, as in Lasslop et al., (2010). β₀ is the canopy light-saturated CO₂ assimilation rate after removing the influence of VPD, which is equivalent to Amax.

Ecosystem respiration F_r is described by the Lloyd and Taylor model as

where R_ref (μmol C m⁻² s⁻¹) is the reference respiration rate at the reference temperature (T_ref = 15°C) and E₀ (°C) is the temperature sensitivity. T_air is the air temperature and the parameter T₀ (°C) is set to a constant of −46.02°C following Lloyd and Taylor (1994).

To get the seasonal variations of the parameters in the above model, we fitted the equations above to a dynamic window of 2–14 days of Fc depending on the availability of flux measurements, and assumed that every day in the same time window has the same daily parameters. We retrieved the time series of all parameters in the above model (i.e., α, β₀, k, R_ref, and E₀) for each site by implementing Equations (1)–(3) using the REddyProc R package. ² The estimated Amax time series were smoothed with a nonparametric method for estimating regression surfaces (LOESS).

The light-saturated photosynthetic rate per unit leaf area (Amax_a) was calculated as the ratio of Amax to LAI. The unit of Amax is μmol m⁻² ground area s⁻¹ and the unit of Amax_a is μmol m⁻² leaf area s⁻¹. Then, the reason why the ratio of Amax to LAI is Amax_a is demonstrated in Equation (4):

We used linear regression coefficients to quantify the environmental control on the seasonal variations in Amax, Amax_a, and LAI. Environmental variables included air temperature (Ta), solar radiation (Rad), SWC, VPD, and day length (Dlen). Day length was defined as the number of hours with Rad >4 W m⁻² each day. The effects of other environmental factors were also analyzed, but the other environmental factors did not show significant effects on Amax, Amax_a, or LAI. For Amax_a, we also added NDVI into the independent variables, because environmental variables had little explanatory power for Amax_a.

Data used for the seasonal analysis were integrated for every adjacent and non-overlapped eight-day window. Because the time series were multiyear, we divided the 8-day time series each year by the annual mean of that year to remove the inter-annual variations in the variables (Baldocchi et al., 2021). To make the regression coefficients of different variables comparable, we applied z-score standardization to the normalized time series before we conducted the linear regression.

The log–log scaling slope analysis was used to quantify the relative contribution of Amax_a and LAI to variations in Amax because Amax_a and LAI are related to Amax in a multiplicative manner (i.e., Amax = Amax_a × LAI) (Renton and Poorter, 2011). The relative contributions of Amax_a and LAI can be calculated as the regression coefficient of log(Amax) in the following equations:

where a_Amaxa, a_LAI, b_Amaxa, and b_LAI are coefficients of the linear regression between Amax and Amax_a and LAI, in which b_Amaxa and b_LAI represent the relative contribution of Amax_a and LAI to the variations of Amax. b_Amaxa + b_LAI will always equal 1. If b_Amaxa is close to 1 and b_LAI is close to 0, then Amax_a is largely responsible for variations in Amax, whereas if b_Amaxa is close to 0 and b_LAI is close to 1, then LAI is largely responsible for variations in Amax. If b_Amaxa and b_LAI are both close to 0.5, then Amax_a and LAI are similarly responsible. More information about the derivation and effectiveness of this method can be found in Renton and Poorter (2011).

The seasonal dynamics of Amax, Amax_a, and LAI and their multiyear averages for each site are shown in Figure 1 and Table 2, respectively. Amax had strong seasonal variations, mainly showing unimodal curves over the season. Amax at North Hemisphere sites usually exhibited a cyclical character each year, but some grasslands in Australia had variable temporal patterns of Amax, which were different for each year, such as at AU-Emr, AU-Stp, and AU-TTE. The seasonal pattern of Amax_a varied substantially among sites. CH-Cha and IT-MBo did not exhibit any clear seasonal trend in Amax_a, while Amax_a at other sites had strong seasonal variations, and often had a similar trend with Amax over the growing season. LAI also varied a lot over the season.

The climate influenced the seasonal patterns of Amax, Amax_a, and LAI. In grasslands located in no dry season climate according to Köppen-Geiger climate classification (e.g., AT-Neu, CH-Cha, CH-Fru, IT-MBo, US-AR1, US-ARc, and US-Goo), Amax, Amax_a, and LAI were higher in higher temperature and radiation season. In contrast, in arid or dry summer sites, the seasonal pattern of Amax, Amax_a, and LAI was driven more by water availability (Figure 2). For example, at the US-SRG, US-Var, and US-Wkg sites, Amax was significantly greater during higher water availability periods, but lower in high temperature and radiation periods because of the associated low water availability.

In addition to climate, agricultural management practices also impacted the seasonal variation of Amax, Amax_a, and LAI. At sites under intensive management practices (i.e., AT-Neu, CH-Cha, CH-Fru, IT-MBo, US-Goo), the unimodal curve of Amax and LAI changed by grass cutting and grazing, which can reduce Amax and LAI. As shown in the seasonal patterns of Amax at these sites, Amax decreased after each cutting or grazing event and later gradually recovered, which added fluctuation into the unimodal seasonal pattern of Amax (Figure 1). For example, at the AT-Neu site, the meadow was cut three times a year, the three cuts taking place between DOY 153 ~ 167, DOY 204 ~ 224, and DOY 264 ~ 301 (Wohlfahrt et al., 2008). We can observe decreases and recovery of Amax and Amax_a around the three cutting periods in our estimated Amax and Amax_a time series (Figure 1). After each cutting, the grassland could not recover to the same level of Amax as before cutting. Therefore, Amax was decreasing with each of the three cuttings. A similar pattern was also observed in IT-MBo, CH-Cha, and CH-Fru.

Amax_a and LAI jointly influenced the seasonal variations of Amax (b_Amaxa and b_LAI ≠ 0) at all sites (Table 3). b_LAI > 0.5 at all managed sites, except CH-Cha, meaning that LAI contributed more to the seasonal variations in Amax than Amax_a at these sites, which might be because cutting and grazing caused more changes in LAI than in Amax_a.

We investigated the effect of air dryness on the relative contribution of LAI to Amax, and we found that b_LAI < 0.5 and b_Amaxa > 0.5 in all arid grasslands (mean annual VPD > 10 hPa), which indicated that Amax_a had more influence on Amax than LAI in these sites. In comparison, among moist grassland sites (mean annual VPD < 10 hPa), there were LAI-driven sites and Amax_a-driven sites (Figure 3).

To explain the difference between arid and moist grasslands in the relative contribution of LAI and Amax_a to Amax, we further investigated the relationship of b_Amaxa with the relative magnitude of the seasonal variability of Amax_a and LAI and the relationship of mean annual VPD and SWC with the seasonal variability of Amax_a and LAI (Figure 4). We found that at sites where CV_Amaxa > CV_LAI, Amax_a contributed more to the seasonal variations of Amax than LAI (b_Amaxa > 0.5), and b_Amaxa < 0.5 at sites where CV_Amaxa < CV_LAI (Figure 4). Sites with CV_Amaxa > CV_LAI tended to have higher mean annual VPD (Figure 4A). Therefore, Amax_a dominated the seasonal variation of Amax at sites with higher mean annual VPD. Mean annual SWC did not show any correlation with the relative magnitude of CV_Amaxa and CV_LAI, but grasslands with higher mean annual SWC had lower absolute CV_Amaxa (Figure 4B).

Among the environmental variables, air temperature (Ta) showed strong control effects for the seasonal variations of Amax_a, LAI, and Amax (Table 4), with Ta having the largest absolute regression coefficient with them in most sites. Air dryness (VPD) also greatly influenced the seasonal variations of Amax_a, LAI, and Amax, especially for grasslands with a mean annual VPD > 10 hPa. It was consistently negatively associated with them at all sites. However, SWC, representing soil water condition, had a much smaller impact on Amax_a, LAI, and Amax. It positively correlated with them at all arid sites, but negatively correlated with them at some moist sites. Solar radiation and photoperiod were positively correlated with Amax_a, LAI, and Amax at many sites.

Environmental variables explained on average 38, 70, and 62% of the seasonal variations in Amax_a, LAI, and Amax, respectively, which suggested environmental variables had a greater impact on LAI than on Amax_a. Therefore, environmental factors affected the seasonal variations of Amax mainly through its impact on LAI rather than through its impact on Amax_a.

Environmental variables had little explanatory power for Amax_a, except a few arid grasslands (i.e., AU-DaP, AU-TTE, US-Var, US-Wkg) where environmental factors can account for over 50% variance of Amax_a. The low average R² across sites suggested that environmental factors could affect Amax_a, but were not the main drivers for the seasonal variations of Amax_a. We further investigated whether vegetation greenness (represented by NDVI) was better correlated with Amax_a and found that the average R² for Amax_a with NDVI as one of the independent variables was 0.44, which was slightly greater than without NDVI. Coefficients of NDVI for Amax_a were significant at 11 of the 15 sites. These results suggested that the seasonal variations of Amax_a cannot be well explained by classical environmental factors, but may be more correlated with variations in plant traits such as pigment contents.

Our results suggested that Amax_a and LAI jointly determined the seasonal dynamics of Amax. The seasonal variation of Amax at all arid grasslands (mean annual VPD > 10 hPa) was driven more by changes in Amax_a than LAI, but the drivers of Amax for the moist grasslands (mean annual VPD < 10 hPa) were bimodal in that some sites were Amax_a-driven and others LAI-driven. The relative contribution of Amax_a and LAI to Amax was related to the relative magnitude of the temporal variability (CV) of Amax_a and LAI, which in turn, was influenced by mean annual VPD.

Annual air dryness influenced the relative magnitude of CV_Amaxa and CV_LAI from two aspects. On one hand, mean annual VPD was negatively correlated with LAI (the correlation coefficient between them was-0.77, p < 0.001), suggesting that drier grasslands tended to have lower LAI, and hence the magnitude of temporal variation in LAI would be smaller. On the other hand, VPD could affect leaf photosynthetic capacity directly by its control on leaf stomatal conductance (Fu et al., 2022), and indirectly by its impact on leaf traits. Moreover, environmental variables explained a large proportion of the variance in Amax_a in arid grasslands (Table 4) compared with wetter grasslands, which also suggested Amax_a in arid grasslands were more sensitive to changes in the environment. Therefore, lesser LAI and greater sensitivity of Amax_a to environmental variables in arid grasslands resulted in a greater contribution of Amax_a to Amax.

At the seasonal scale, air temperature and air dryness were the driving factors for the variation in Amax. Other studies also identified the two factors as the most important factors to regulate grassland C fluxes (Polley et al., 2010; Scott et al., 2015; Wang et al., 2019). We also found that VPD was more correlated with Amax than SWC at the seasonal scale, which was also observed in previous studies (Ryan et al., 2017; Fu et al., 2022). The possible reason for this higher correlation of VPD and Amax at the seasonal scale could be related to the non-linear relationship between SWC and Amax, as suggested by previous studies (Xu and Zhou, 2011; Fu et al., 2022). Xu and Zhou (2011) found that the Amax of two grass species was enhanced under moderate soil moisture with reductions under both severe water deficit and excessive water conditions. Fu et al. (2022) suggested that Amax was also non-linearly correlated with soil moisture at the ecosystem scale, while VPD was consistently negatively related to Amax.

The environmental factors used in our study can explain the variances in LAI better than Amax and Amax_a. The explanatory power (R²) of linear regression can be impacted by intensive management practices or other irregular natural events, such as fire and extreme temperatures (Moore et al., 2016). For some sites, such as AT-Neu, IT-MBo, CH-Cha, and CH-Fru, previous studies at these sites have demonstrated that management practices, like cutting, grazing, and fertilization, can obscure the response of photosynthesis to environmental variables (Wohlfahrt et al., 2008; Zeeman et al., 2010; Marcolla et al., 2011; Rogger et al., 2022).

In our study, environmental variables only explained a small proportion of variance in Amax_a on the 8-day scale, much less than that of Amax and LAI, which is probably because canopy photosynthetic capacity per leaf area (Amax_a) was controlled more by leaf traits, such as nitrogen (N) and chlorophyll (Chl) content, than by environmental factors, compared with Amax and LAI. According to its definition in our study, Amax_a is the average of the maximum photosynthetic rate per leaf area of individual plants in the ecosystem, or more precisely, all plants within the footprint of the flux tower, a variable independent of the ecosystem structure embedded in LAI and intrinsically related to individual plants. Reichstein et al. (2014) have summarized that environmental variables can account for only a fraction of the temporal and spatial variance of flux-derived ecosystem functional properties (EFP), e.g., seasonal carbon-use efficiency and light-saturated gross primary production (similar to Amax), and suggested that a substantial part of the variation in EFP would be explained by plant traits. Furthermore, Feng and Dietze (2013) examined the relationship between the seasonal variation of Amax and leaf traits on 25 grassland species. They found that leaf N and Chl content significantly correlated with leaf Amax in both mass-and area-based manner. Therefore, we thought that the large unexplained variance in Amax_a was related to plant traits.

There are still some uncertainties in our estimation and analysis. The LAI data and Amax estimates used in this study had some uncertainty, which could propagate into Amax_a. To reduce the uncertainty of LAI, we tested whether the radius we used to extract LAI from the MODIS LAI dataset could have impacted the temporal dynamics of LAI. We extracted the average LAI of pixels within a 250, 500, 1,000, and 3,000 m radius centered at the flux tower. We found that these buffer sizes did not have much impact on the temporal dynamics in LAI. Then we assessed the influence of the uncertainty of the LAI dataset on Amax_a. Following He et al. (2014), we calculated the relative uncertainty in Amax_a induced by the uncertainty in LAI, which was expressed as the standard deviation for each LAI observation provided by the MODIS LAI dataset. The result (Table 5) showed that on average, the uncertainty in LAI could introduce 27.28% of uncertainty in Amax_a.

The uncertainty of Amax estimates was assessed by the fitting quality flag provided by the REddyProc package. The results were shown in Figure 5. On average, each site had 52% of the Amax estimates that were of good quality, 48% of moderate quality, and no bad quality estimates.