Asymmetric responses of primary productivity to altered precipitation simulated by ecosystem models across three long-term grassland sites

Field measurements of aboveground net primary productivity (ANPP) in temperate grasslands suggest that both positive and negative asymmetric responses to changes in precipitation (P) may occur. Under normal range of precipitation variability, wet years typically result in ANPP gains being larger than ANPP declines in dry years (positive asymmetry), whereas increases in ANPP are lower in magnitude in extreme wet years compared to reductions during extreme drought (negative asymmetry). Whether the current generation of ecosystem models with a coupled carbon–water system in grasslands are capable of simulating these asymmetric ANPP responses is an unresolved question. In this study, we evaluated the simulated responses of temperate grassland primary productivity to scenarios of altered precipitation with 14 ecosystem models at three sites: Shortgrass steppe (SGS), Konza Prairie (KNZ) and Stubai Valley meadow (STU), spanning a rainfall gradient from dry to moist. We found that (1) the spatial slopes derived from modeled primary productivity and precipitation across sites were steeper than the temporal slopes obtained from inter-annual variations, which was consistent with empirical data; (2) the asymmetry of the responses of modeled primary productivity under normal inter-annual precipitation variability differed among models, and the mean of the model ensemble suggested a negative asymmetry across the three sites, which was contrary to empirical evidence based on filed observations; (3) the mean sensitivity of modeled productivity to rainfall suggested greater negative response with reduced precipitation than positive response to an increased precipitation under extreme conditions at the three sites; and (4) gross primary productivity (GPP), net primary productivity (NPP), aboveground NPP (ANPP) and belowground NPP (BNPP) all showed concave-down nonlinear responses to altered precipitation in all the models, but with different curvatures and mean values. Our results indicated that most models overestimate the negative drought effects and/or underestimate the positive effects of increased precipitation on primary productivity under normal climate conditions, highlighting the need for improving eco-hydrological processes in those models in the future.

1 Introduction

Precipitation (P) is a key climatic determinant of ecosystem productivity, especially in arid and semi-arid grasslands (Lambers et al., 2008; Sala et al., 1988; Hsu et al., 2012; Beer et al., 2010). Climate models project substantial changes in amounts and frequencies of precipitation regimes worldwide, and this is supported by observational data (Karl and Trenberth, 2003; Donat et al., 2016; Fischer and Knutti, 2016). Potential for increasing occurrence and severity of droughts and increased heavy rainfall events related to global warming will likely affect grassland growth (Knapp et al., 2008, 2017a; Gherardi and Sala, 2015; Lau et al., 2013; Reichstein et al., 2013). As a consequence, better understanding of the responses of grassland productivity to altered precipitation is needed to project future climate–carbon interactions, changes in ecosystem states, and to gain better insights on the role of grasslands in supporting crucial ecosystem services (e.g., livestock production).

Gross primary productivity (GPP) of ecosystems is controlled by environmental conditions, in particular water availability (Jung et al., 2017), and by biotic factors affecting leaf photosynthetic rates and stomatal conductance, which scale up to canopy-level functioning (Chapin III et al., 2011). About half of GPP is respired while the remainder, net primary productivity (NPP), is primarily invested in plant biomass production, including photosynthetic and structural pools aboveground (foliage and stem) and belowground (roots) (Waring et al., 1998; Chapin III et al., 2011). NPP responses to precipitation have been observed using multi-year, multi-site observations (Hsu et al., 2012; Estiarte et al., 2016; Knapp and Smith, 2001; Wilcox et al., 2015). Positive empirical relationships between grassland aboveground NPP (ANPP) and precipitation (P) have been found in spatial gradients across sites (Sala et al., 1988) and from temporal variability at individual sites (Huxman et al., 2004; Knapp and Smith, 2001; Roy et al., 2001; Hsu et al., 2012). The ANPP–P sensitivities obtained from spatial relationships are usually higher than those obtained by temporal relationships (Estiarte et al., 2016; Fatichi and Ivanov, 2014; Sala et al., 2012). Possible mechanisms behind the steeper spatial relationship may be (1) a “vegetation constraint” reflecting the adaptation of plant communities over long timescales in such a way that grasslands make the best use of the typical water received from rainfall for growth (Knapp et al., 2017b) and (2) the spatial variation in structural and functional traits of ecosystems (soil properties, nutrient pools, plant and microbial community composition) that constrain local ANPP–P sensitivities (Lauenroth and Sala, 1992; Smith et al., 2009; Wilcox et al., 2016). For projecting the effect of climate change on grassland productivity in the near to mid-term (coming decades), inter-annual relationships are arguably more informative than spatial relationships because spatial relationships reflect long-term adaptation of ecosystems and because ANPP–P relationships from spatial gradients are confounded by the covariation of gradients in other environmental variables (e.g., temperature and radiation) and soil properties (Estiarte et al., 2016; Knapp et al., 2017b).

In temporal ANPP–P relationships, an important observation is the asymmetric responses of productivity in grasslands to altered precipitation (Knapp et al., 2017b; Wilcox et al., 2017). Compared to negative anomalies of ANPP from years with decreased precipitation, positive anomalies of ANPP during years with increased precipitation were usually found to have a larger absolute magnitude, suggesting a convex positive response (positive asymmetry) (Bai et al., 2008; Knapp and Smith, 2001; Yang et al., 2008). Yet, when grasslands are subject to extreme precipitation anomalies that fall beyond the range of normal inter-annual variability, an extreme dry year is associated with a larger absolute ANPP loss than the gain found during an extreme wet year. This suggests a convex negative response (negative asymmetry) when considering a larger range of rainfall anomalies than the current inter-annual regime (Knapp et al., 2017b). This is also supported by current dynamical global vegetation models, which suggest a stronger response to extreme dry conditions compared to extreme wet conditions (Zscheischler et al., 2014). The sign of the asymmetric response of grassland productivity to altered rainfall thus depends on the magnitude of rainfall anomalies, the size distribution of rainfall events and ecosystem mean state (Gherardi and Sala, 2015; Hoover and Rogers, 2016; Parolari et al., 2015; Peng et al., 2013).

Table 1 Key plant, soil and climate characteristics of the three grassland sites. MAT, mean annual temperature; and MAP, mean annual precipitation. MAT and MAP are based on the periods for the three sites with ANPP measurements.

Relationships between precipitation and grassland productivity have previously been studied with site observations (Hsu et al., 2012; Knapp et al., 2017b; Luo et al., 2017; Wilcox et al., 2017; Estiarte et al., 2016), but they remain to be quantified and characterized in ecosystem models used for diagnostic and future projections of the coupled carbon–water system in grasslands, in particular grid-based models used as the land surface component of Earth system models. In this study, we aim to evaluate the responses of simulated productivity to altered precipitation from 14 ecosystem models at three sites representing dry (304 ± 118 mm yr⁻¹), mesic (827 ± 175 mm yr⁻¹) and moist (1429 ± 198 mm yr⁻¹) rainfall regimes. The specific objectives of this study are to (1) test if the productivity–P sensitivities of spatial relationships are greater than the temporal ones in the models such as those found in the observations; (2) test if models reproduce the observed asymmetric responses under inter-annual precipitation conditions; (3) assess the simulated productivity–P sensitivities related to different precipitation regimes including normal and extreme conditions, and to test in particular if sensitivities for extreme drought conditions are stronger than those for high-rainfall conditions; (4) analyze the simulated curvilinear productivity–P relationships for a large range of altered precipitation amounts across the three sites.

2 Materials and methods

2.1 Experimental sites

We conducted model simulations using three sites: the Shortgrass steppe (SGS) site at the Central Plains Experimental Range, the Konza Prairie Biological Station (KNZ) site and the Stubai Valley meadow (STU) site. These sites represent three grassland types spanning a productivity gradient from dry to moist climatic conditions. The dry SGS site is located in northern Colorado, USA (Knapp et al., 2015; Wilcox et al., 2015). The KNZ site is a native C₄-dominated mesic tallgrass prairie in the Flint Hills of northeastern Kansas, USA (Heisler-White et al., 2009; Hoover et al., 2014). The moist site of STU is a subalpine meadow located in the Austrian Central Alps near the village of Neustift (Bahn et al., 2006, 2008; Schmitt et al., 2010). Experimental measurements of annual ANPP were carried out spanning different time ranges. Estimated mean ANPP for SGS, KNZ and STU sites are 91 ± 36, 387 ± 82 and 525 ± 210 g DM (dry mass) m⁻² yr⁻¹. Details of the ecological and environmental factors are summarized in Table 1.

These three grasslands were selected because they lie along a mean annual precipitation (MAP) gradient and have detailed meteorological data to force the models. While two are “natural” grasslands (KNZ and SGS) and one (STU) is not, global land surface models do not typically differentiate regarding the origin of ecosystem types and heavily managed grasslands and pastures represent a significant fraction of mesic grasslands globally. Semi-natural subalpine grasslands in the Alps were created several centuries ago, are very lightly managed and should be in equilibrium concerning soil physical conditions. It should be noted though that the grassland at STU is cut once a year and lightly fertilized every 2–4 years and in consequence differs in plant composition and soil fungi : bacteria ratio, which leads to different drought responses compared to abandoned grassland (Ingrisch et al., 2017; Karlowsky et al., 2018). Further, it is worth noting that the mesic grassland in the USA would also be forested if human-initiated prescribed fires were to be removed from the system (Briggs et al., 2005). Thus, these grassland sites lie along a continuum of dry natural grassland, mesic natural grassland maintained by human management and anthropogenic moist grassland maintained by human management.

2.2 Ecosystem model simulations

In order to test the hypothesis of an asymmetric response of productivity to variable rainfall (Knapp et al., 2017b), simulations were conducted with 14 ecosystem models – CABLE, CLM45-ORNL, DLEM, DOS-TEM, JSBACH, JULES, LPJ-GUESS, LPJmL-V3.5, ORCHIDEE-2, ORCHIDEE-11, T&C, TECO, TRIPLEX-GHG and VISIT – all using the same protocol defined by the precipitation subgroup of the model–experiment interaction study (Table 2). At all three grassland sites, observed and altered multi-annual hourly rainfall forcing time series were combined with observations of other climate variables. These variables were air temperature, incoming solar radiation, air humidity, wind speed and surface pressure. Model simulations were carried out using soil texture properties measured at each site as reported in Table 1. Simulated productivity during the observational period is influenced at least in some models (for instance those having C–N interactions) by historical climate change and CO₂ changes since the preindustrial period. Thus, instead of assuming that productivity was in equilibrium with current climate, historical reconstructions of meteorological variables from gridded CRUNCEP data at half-hourly time step (Wei et al., 2014) were combined and bias corrected with site observations to provide bias corrected historical forcing time series from 1901 to 2013 (CRUNCEP-BC). In addition to the observed current climate defining the ambient simulation, nine altered rainfall forcing datasets were constructed by decreasing or increasing the amount of precipitation in each precipitation event by −80, −70, −60, −50, −20, +20, +50, +100 and +200 % during the time span of productivity observations at each site, leaving all other meteorological variables unchanged and equal to the observed values. Modelers performed all simulations described below based on the same protocol (see below) and the model output was compared with measured ecosystem productivities (GPP; NPP; ANPP; and BNPP, belowground NPP), whenever available.

Table 2 Summary of ecosystem models used in this study, including model name, nitrogen (N) cycle and relevant references. Also see Tables S1–S14 in the Supplement for details of the simulated processes for grasslands in the ecosystem models, including the N cycle, phosphorus (P) cycle, carbon (C) allocation scheme, carbohydrate reserves, leaf photosynthesis and stomatal conductance including treatment of water stress, scaling of photosynthesis from leaf to canopy, phenology, mortality, soil hydrology, surface energy budget, root profile and dynamics, and grassland species.

Simulation S0 spin-up: models simulated an initial steady state spin-up run for water and biomass pools under preindustrial conditions using the 1901–1910 CRUNCEP-BC climate forcing in a loop and applying fixed atmospheric CO₂ concentration at the 1850 level.

Simulation S1 historical simulation from 1850 until the first year of measurement (1986 for SGS, 1982 for KNZ and 2009 for STU): starting from the spin-up state, models were prescribed with increasing atmospheric CO₂ concentrations and dynamic historical climate from CRUNCEP-BC. Because there is no CRUNCEP-BC data for 1850–1900, the CRUNCEP-BC climate data from 1901 to 1910 was repeated in a loop instead.

Simulation SC1 ambient simulation for the measurement periods (1986–2009 for SGS, 1982–2012 for KNZ and 2009–2013 for STU): starting from the initial state in the start year of the period and run with observed CO₂ concentrations and meteorological data corresponding to site observations at the hourly or half-hourly scale.

Simulations SP1–SP9 altered precipitation simulations for the measurement periods (1986–2009 for SGS, 1982–2012 for KNZ and 2009–2013 for STU): starting from the initial state in the start year of the period and run using the nine altered rainfall forcing datasets with observed CO₂ concentration.

2.3 Metrics of the response of productivity to precipitation changes

In the analysis, we begin with testing our first specific objective, i.e., if the productivity–P sensitivities of spatial relationships are greater than the temporal ones in the models as found in the observations. We calculated the temporal slopes and spatial slopes between productivities and precipitation from multi-year ambient simulations (SC1). Temporal slopes are site based and relate inter-annual variability in precipitation to inter-annual variability in the productivities using linear regression analysis. Spatial slopes relate mean annual precipitation to mean annual productivity across the three sites.

We then calculated two indices to analyze the asymmetric responses of primary productivity to precipitation simulated by ecosystem models and derived by observations whenever data were available. The two indices are (1) the asymmetry of productivity–P for current inter-annual variability, based on SC1 where observations for ANPP are also available; and (2) the sensitivity of productivity to P for simulations where mean precipitation was altered, based on SP results. With these metrics, we test our second and third specific objectives, i.e., whether models could reproduce the observed asymmetric responses of productivity in grasslands to altered precipitation under normal and extreme conditions.

Finally, we analyze the nonlinearity of modeled response of productivity to precipitation, which is described by the parameters of the curvilinear productivity–P relationships across the full range of altered precipitation scenarios, based on fits to model output for the ambient (SC1) and altered (SP) simulations. Detailed methods for the two indices used to analyze the asymmetric responses of primary productivity to altered precipitation and the curvilinear productivity–P relationships are introduced in the following.

2.3.1 Asymmetry index from inter-annual productivity and precipitation

In order to characterize the asymmetry of productivity to precipitation, we define the asymmetry index (AI) from inter-annual productivity and precipitation data as follows:

$$\begin{array}{}\text{(1)}& \mathrm{AI}=R_{\mathrm{p}}-R_{\mathrm{d}},\end{array}$$

where R_p is the relative productivity pulse in wet years and R_d is the relative productivity decline in dry years defined by

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{p}}=\left(\mathrm{med}\right(f_{\mathrm{p}\mathrm{90}})-\stackrel{\mathrm{‾}}{f})/\stackrel{\mathrm{‾}}{f},\text{(3)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{d}}=(\stackrel{\mathrm{‾}}{f}-\mathrm{med}(f_{\mathrm{p}\mathrm{10}}\left)\right)/\stackrel{\mathrm{‾}}{f},\end{array}$$

where f is the inter-annual productivity, being a function of environmental factors from models or observation; $\stackrel{\mathrm{‾}}{f}$ is mean annual productivity in the period of measurements (Table 1); med(f_p90) is the median value of productivities in wet years with annual precipitation higher than the 90th percentile level; and med(f_p10) is median value of productivities in all the dry years when annual precipitation is lower than the 10th percentile level.

In general, R_p > 0 indicates that the median value of productivities in wet years is higher than the mean annual productivity in the period of measurements; and R_d > 0 indicates that the median value of productivities in dry years is smaller than the mean annual productivity in the period of measurements. Therefore, AI > 0, i.e., a positive asymmetry, means that there is a greater increase of productivity in wet years than decline in dry years; and AI < 0, i.e., a negative asymmetry, means that there is a greater decline of productivity in dry years than increase in wet years.

Furthermore, uncertainty ranges of R_p, R_d and AI were estimated as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{R}_{\mathrm{p}}\in \left[R_{{\mathrm{p}}_{\mathrm{low}}},R_{{\mathrm{p}}_{\mathrm{up}}}\right]=[{\displaystyle \frac{(\mathrm{med}\left(f_{\mathrm{p}\mathrm{90}}\right)-\mathrm{mad}(f_{\mathrm{p}\mathrm{90}}\left)\right)-\stackrel{\mathrm{‾}}{f}}{\stackrel{\mathrm{‾}}{f}}},\\ \text{(4)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{(\mathrm{med}\left(f_{\mathrm{p}\mathrm{90}}\right)+\mathrm{mad}(f_{\mathrm{p}\mathrm{90}}\left)\right)-\stackrel{\mathrm{‾}}{f}}{\stackrel{\mathrm{‾}}{f}}}],{\displaystyle }& {\displaystyle }{R}_{\mathrm{d}}\in \left[R_{{\mathrm{d}}_{\mathrm{low}}},R_{{\mathrm{d}}_{\mathrm{up}}}\right]=[{\displaystyle \frac{\stackrel{\mathrm{‾}}{f}-(\mathrm{med}\left(f_{\mathrm{p}\mathrm{10}}\right)+\mathrm{mad}(f_{\mathrm{p}\mathrm{10}}\left)\right)}{\stackrel{\mathrm{‾}}{f}}},\\ \text{(5)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\stackrel{\mathrm{‾}}{f}-(\mathrm{med}\left(f_{\mathrm{p}\mathrm{10}}\right)-\mathrm{mad}(f_{\mathrm{p}\mathrm{10}}\left)\right)}{\stackrel{\mathrm{‾}}{f}}}],\text{(6)}& {\displaystyle }& {\displaystyle }\mathrm{AI}\in \left[{\mathrm{AI}}_{\mathrm{low}},{\mathrm{AI}}_{\mathrm{up}}\right]=[R_{{\mathrm{p}}_{\mathrm{low}}}-R_{{\mathrm{d}}_{\mathrm{up}}},R_{{\mathrm{p}}_{\mathrm{up}}}-R_{{\mathrm{d}}_{\mathrm{low}}}],\end{array}$$

where $R_{{\mathrm{p}}_{\mathrm{low}}}$ and $R_{{\mathrm{p}}_{\mathrm{up}}}$ are the lower and upper bounds of R_p using one median absolute deviation, i.e., mad(f_p90); $R_{{\mathrm{d}}_{\mathrm{low}}}$ and $R_{{\mathrm{d}}_{\mathrm{up}}}$ are the lower and upper bounds of R_d using one median absolute deviation, i.e., mad(f_p10); and AI_low and AI_up are the lower and upper bounds of AI corresponding to estimated R_p and R_d ranges.

2.3.2 Sensitivity of productivity to altered versus inter-annual precipitation variability

For altered precipitation, in particular for the extreme SP simulations where mean precipitation was altered and annual precipitation of a few years was outside the range of observed precipitation variation, we tested the hypothesis of whether the asymmetry response becomes negative – that is the impacts of extreme dry conditions on productivity are much greater than the positive effects of extreme wet scenarios (Knapp et al., 2017b). Thus, we tested the mean change in productivity imposed by the change in precipitation, and we defined the sensitivity of productivity to altered rainfall conditions (S) as

$$\begin{array}{}\text{(7)}& {\displaystyle }S=(\stackrel{\mathrm{‾}}{f_{P_{\mathrm{a}}}}-\stackrel{\mathrm{‾}}{f_{P_{\mathrm{c}}}})/\left(\left|\stackrel{\mathrm{‾}}{P_{\mathrm{a}}}-\stackrel{\mathrm{‾}}{P_{\mathrm{c}}}\right|\right),\end{array}$$

where $\stackrel{\mathrm{‾}}{f_{P_{\mathrm{a}}}}$ and $\stackrel{\mathrm{‾}}{f_{P_{\mathrm{c}}}}$ are the mean productivities of altered and ambient simulations; $\stackrel{\mathrm{‾}}{P_{\mathrm{a}}}$ and $\stackrel{\mathrm{‾}}{P_{\mathrm{c}}}$ are the mean annual precipitation amounts in altered and ambient simulations. It should be noted that the sensitivity of productivity to altered rainfall conditions could present the asymmetry response from normal to extreme conditions.

Figure 1 Relationships between GPP (a), NPP (b), ANPP (c), and BNPP (d) and precipitation (P) derived from multi-year ambient simulations (SC1) in two ways. Temporal slopes are site based and relate inter-annual variability in P to inter-annual variability in the productivities using linear regression analysis. Spatial slopes relate mean annual P to mean annual productivity across three sites. In each panel, SGS, KNZ and STU are from dry to moist, given from left to right. The red lines are the ensemble mean of modeled temporal slopes, and the red shading represents the model uncertainty range using the interquartile spread of the temporal slopes between individual simulations (10th and 90th percentiles). The blue line is the ensemble mean of modeled productivities, and the blue error bar represents the model uncertainty range using the interquartile spread of the productivities between individual simulations (10th and 90th percentiles). In (c), the grey lines are the observed temporal slopes, and the black line shows the observed spatial slope. The grey shading represents the observed uncertainty range using the bootstrap sampling method (10th and 90th percentiles), and the black error bar represents the observed uncertainty range using the interquartile spread of the inter-annual productivities (10th and 90th percentiles). Note that we simply converted observed ANPP from dry mass (g DM m⁻² yr⁻¹) to carbon mass (g C m⁻² yr⁻¹) with a factor of 0.5.

2.3.3 Curvilinear productivity–P relationships across the entire range of altered P

In general, plant productivity increases with increasing precipitation and saturates when photosynthesis becomes less limited by water scarcity. We fitted the response of simulated productivity to altered precipitation using the Eq. (8):

$$\begin{array}{}\text{(8)}& y=a\left(\mathrm{1}-e^{-bx}\right),\end{array}$$

where the independent variable x is the mean annual precipitation (mm) and the dependent variable y one of the productivities (GPP, NPP, ANPP and BNPP). Parameter a (g C m⁻² yr⁻¹) is the maximum value of productivity at high precipitation and parameter b (mm⁻¹) is the curvature of modeled productivity to altered precipitation.

Figure 2 Asymmetry responses of inter-annual GPP (a), NPP (b), ANPP (c) and BNPP (d) to precipitation in ambient simulations at the three sites SGS, KNZ and STU. The asymmetry index was calculated as the difference between the relative productivity pulses (R_p) and declines (R_d) in wet years and dry years (see Eqs. 1–3). Black pentagrams in (c) represent asymmetry indices from observations. The corresponding black error bars represent the observed uncertainty ranges using Eqs. (4)–(6). A black asterisk at the bottom of a panel indicates a significant asymmetry response of the model ensemble at a 0.1 significance level by a non-parametric statistical hypothesis test (Wilcoxon signed-rank test).

3 Results

3.1 Temporal versus spatial slopes of productivity–P

From the ambient simulations, ensemble model results indicate that the slopes of the spatial relationships were steeper than the temporal slopes for GPP, NPP and ANPP for the subset of models that simulated this flux, while these differences in slopes were less obvious for BNPP (Fig. 1). We compared model results with site observations for ANPP–P temporal slopes of the ambient simulation across the three sites (Fig. 1c). Observed and modeled temporal slopes decreased from the dry (SGS) to moist (STU) site, from 0.10 g C m⁻² mm⁻¹ (0.05 to 0.14 for the 10th and 90th percentiles) to 0.05 g C m⁻² mm⁻¹ (−0.14 to 0.55 for the 10th and 90th percentiles) in the observations, and from 0.14 g C m⁻² mm⁻¹ (0.02 to 0.36 for the 10th and 90th percentiles) to 0.03 g C m⁻² mm⁻¹ (−0.04 to 0.29 for the 10th and 90th percentiles) for the model ensemble mean. Although there were some discrepancies in the range of spatial and temporal slopes across models (Fig. S1 in the Supplement), the multi-model ensemble mean captured the key observation of spatial slopes steeper than temporal slopes for ANPP (Fig. 1).

Figure 3 Sensitivity of GPP (a), NPP (b), ANPP (c) and BNPP (d) for altered precipitation simulations at the three sites SGS, KNZ and STU. Curves show the ensemble mean of models, and the shading represents the model uncertainty range using the interquartile spread of the sensitivities between individual simulations (10th and 90th percentiles). Curves above the zero line represent responses under increasing precipitation conditions relative to the control, and curves below the zero line show responses under decreasing precipitation conditions relative to the control. Vertical dashed lines represent precipitation variations of 1 standard deviation (1σ), 2 standard deviations (2σ) and 3 standard deviations (3σ), which were derived from long-term annual precipitation at the three sites respectively.

3.2 Asymmetry of the inter-annual primary productivity response to precipitation

The asymmetry of each model was diagnosed using the asymmetry index (Eq. 1), which showed large variation across models (Figs. 2, S2). Considering all the models as independent ensemble members, the mean AI of GPP and NPP showed significantly negative values at p < 0.1 level for SGS (ensemble value of −0.11${}_{-\mathrm{0.31}}^{\mathrm{0.12}}$ and −0.20${}_{-\mathrm{0.48}}^{\mathrm{0.11}}$ respectively with 10th and 90th percentiles). Hence, for SGS simulated declines of GPP and NPP in dry years were larger than the increases in wet years. For STU, the mean AI values were only slightly negative (ensemble value for GPP −0.03${}_{-\mathrm{0.07}}^{\mathrm{0.02}}$ and for NPP −0.04${}_{-\mathrm{0.09}}^{\mathrm{0.01}}$ with 10th and 90th percentiles), while AI was very close to zero at KNZ. By contrast, observation-based AI values, estimated from long-term inter-annual ANPP measurements, suggest a decrease from positive (${\mathrm{0.32}}_{\mathrm{0.14}}^{\mathrm{0.49}}$ for SGS and 0.20${}_{\mathrm{0.04}}^{\mathrm{0.37}}$ for KNZ) to negative (−0.21 for STU). At the dry (SGS) and mesic (KNZ) sites (Fig. S2), most of the model simulations overestimated the extent of negative drought effects in dry years (R_d) and/or underestimated the positive impacts on ANPP in wet years (R_p). For example, CABLE and ORCHIDEE-2 overestimated the drought effects in dry years at both of the two sites, and CLM45-ORNL and VISIT underestimated the positive impacts in wet years at both of the two sites (Fig. S2). At the moist site (STU), models agreed with observations regarding the negative sign of AI (negative asymmetry) but AI magnitude is not well captured.

3.3 Sensitivities of primary productivity to altered precipitation

The model-derived sensitivities given by Eq. (7) generally presented greater negative impacts of reduced precipitation than positive effects of increased precipitation under both normal (inter-annual) and extreme conditions (Fig. 3). The results also indicated that models represented a constant asymmetry pattern (negative asymmetry under normal and extreme conditions) across the full range of altered precipitation rather than a double asymmetry pattern (positive asymmetry under normal condition and negative asymmetry under extreme condition) established by Knapp et al. (2017b), which confirmed that models did not capture the positive asymmetric responses of productivities to altered precipitation under normal conditions for the dry (SGS) and mesic (KNZ) sites.

Primary productivity at the dry site (SGS) was more sensitive to precipitation changes compared to the moist site (STU). Along with increases in precipitation, the largest sensitivity values were found for SGS (ensemble mean of ${\mathrm{1.35}}_{\mathrm{0.42}}^{\mathrm{2.49}}$ g C m⁻² mm⁻¹ for GPP with 10th and 90th percentiles, ${\mathrm{0.68}}_{\mathrm{0.24}}^{\mathrm{1.47}}$ g C m⁻² mm⁻¹ for NPP, ${\mathrm{0.24}}_{\mathrm{0.08}}^{\mathrm{0.61}}$ g C m⁻² mm⁻¹ for ANPP and ${\mathrm{0.16}}_{\mathrm{0.14}}^{\mathrm{0.18}}$ g C m⁻² mm⁻¹ for BNPP) and then KNZ (${\mathrm{0.32}}_{-\mathrm{0.09}}^{\mathrm{1.23}}$ g C m⁻² mm⁻¹ for GPP, ${\mathrm{0.20}}_{-\mathrm{0.05}}^{\mathrm{0.72}}$ g C m⁻² mm⁻¹ for NPP, ${\mathrm{0.13}}_{\mathrm{0.01}}^{\mathrm{0.21}}$ g C m⁻² mm⁻¹ ANPP and ${\mathrm{0.06}}_{\mathrm{0.01}}^{\mathrm{0.28}}$ g C m⁻² mm⁻¹ for BNPP with 10th and 90th percentiles) when precipitation was altered by +20 %. The values of S decreased with further increased precipitation, indicating that additional water does not increase productivity in the same proportion exceeding a certain threshold. In contrast to SGS, the values of sensitivity for both GPP and NPP at STU are close to zero in response to added precipitation conditions, implying that the precipitation above ambient was not a limiting factor for grassland production in the models at this site.

Figure 4 Responses of simulated annual GPP (a, d, g), NPP (b, e, h) and CUE (NPP/GPP; c, f, i) to altered and ambient precipitation (P) levels at the three sites STU, KNZ and SGS. The fitted equation is Eq. (8) for GPP and NPP (see Fig. S3 for fitted a and b). The grey dashed line represents ambient precipitation. It should be noted that the x-axis scales are different between the sites.

The values of sensitivity decreased with reduced precipitation at KNZ and SGS, indicating larger negative impacts on primary productivity when conditions become drier. For the moist site of STU, primary productivities showed less sensitivity to moderately dry conditions, and sensitivity only increased with more extreme rainfall alterations out of 3σ (∼ 40 % precipitation change). Additionally, the values of S for ANPP were smaller than those of BNPP at KNZ and SGS, while there were no differences between ANPP and BNPP at STU (Fig. 3). Thus, model results suggest that the dry site (SGS) can be particularly vulnerable to altered rainfall compared to the moist site (STU), which was more robust in response to altered rainfall.

Figure 5 Responses of simulated annual ANPP (a, d, g), BNPP (b, e, h), and the ratio of ANPP and NPP (c, f, i) to altered and ambient precipitation (P) levels at the three sites STU, KNZ and SGS. The fitted equation is Eq. (8) for ANPP and BNPP (see Fig. S4 for fitted a and b). The grey dashed line represents ambient precipitation. It should be noted that the x-axis scales are different between the sites.

3.4 Curvilinear responses of productivity to altered precipitation

At SGS and KNZ, simulated GPP and NPP increased with increasing precipitation. In contrast, at the moist STU, most models showed saturation in productivity for precipitation above ambient values (Fig. 4). Along with increasing precipitation, GPP and NPP showed nonlinear concave-down response curves in all models, with different curvatures b and maximum productivity a (Fig. S3). The ensemble mean values of the curvature parameter b fitted from Eq. (8) to each modeled GPP across the full range of altered precipitation are ${\mathrm{5.1}}_{\mathrm{2.7}}^{\mathrm{9.2}}$ × 10⁻³ mm⁻¹ at STU, ${\mathrm{3.3}}_{\mathrm{0.9}}^{\mathrm{8.0}}$ × 10⁻³ mm⁻¹ at KNZ and ${\mathrm{1.4}}_{\mathrm{0.0}}^{\mathrm{2.3}}$ × 10⁻³ mm⁻¹ at SGS with 10th and 90th percentiles (Fig. S3).

The responses of GPP and NPP to altered precipitation were proportional to each other for each model, and as a result changes in carbon use efficiency (CUE) were very small compared to the background CUE differences diagnosed in the ambient simulation (Fig. 4c, f, i). However, JSBACH and LPJmL-V3.5 produced a sharp decline of CUE below ambient precipitation at SGS and KNZ.

Only seven models simulated ANPP and BNPP separately (Fig. 5). The responses of ANPP and BNPP to altered precipitation were similar to those of GPP and NPP. When fitting Eq. (8) to ANPP–P (Fig. S4), the curvatures b ranged from 3.0 × 10⁻³ mm⁻¹ (ORCHIDEE-11) to 9.2 × 10⁻³ mm⁻¹ (TECO) at STU, from 0.7 × 10⁻³ mm⁻¹ (TRIPLEX-GHG) to 6.1 × 10⁻³ mm⁻¹ (VISIT) at KNZ, and from 0.9 × 10⁻³ mm⁻¹ (T&C) to 2.3 × 10⁻³ mm⁻¹ (CLM45-ORNL) at SGS; the modeled maximum values a for ANPP ranged between 173 g C m⁻² yr⁻¹ (VISIT) and 827 g C m⁻² yr⁻¹ (TECO) at STU, 49 g C m⁻² yr⁻¹ (CLM45-ORNL) and 557 g C m⁻² yr⁻¹ (ORCHIDEE-2) at KNZ, and 94 g C m⁻² yr⁻¹ (CLM45-ORNL) and 523 g C m⁻² yr⁻¹ (ORCHIDEE-2) at SGS.

The ANPP : NPP ratio, i.e., aboveground carbon allocation, showed a nonlinear increase (concave-down) with increasing precipitation in ORCHIDEE-2 and ORCHIDEE-11, a nonlinear decrease (concave-up) in T&C due to translocation of C reserves from roots and only minor changes in other models (Fig. 5c, f, i).

4 Discussion

4.1 Comparison of modeled and observed responses of productivity to altered precipitation

Spatial slopes steeper than temporal slopes of ANPP to precipitation are usually explained by two hypotheses: (1) vegetation constraint effects on ANPP responses to precipitation play a more important role in the temporal as compared to the spatial domain (Knapp et al., 2017b; Estiarte et al., 2016); (2) biogeochemistry (mainly resource limitations) and confounding factors (e.g., temperature and radiation), rather than species attributes, constrain community-level ANPP in response to precipitation (Huxman et al., 2004). Thus, the former theory stresses more long-term intrinsic ecosystem properties, while the latter underlines the effects of external environmental factors. The current models tested here captured the relative magnitude of the difference between temporal and spatial slopes (Fig. 1c), which suggested that the models adequately considered the key processes underlying carbon–water interactions across different grassland sites. Only few grassland experiments have assessed BNPP (Luo et al., 2017), leaving the question open of whether the minor differences between temporal and spatial slopes for BNPP responses to precipitation as simulated by the models correspond to experimental observations (Fig. 1d).

The asymmetry index obtained from available long-term ANPP and precipitation observations reported positive values at SGS and KNZ (Fig. 2c), which suggested greater declines of ANPP in dry years than increases in wet years (Knapp and Smith, 2001). Knapp et al. (2017b) proposed the following underlying mechanisms. (1) In dry years, the carryover effects of soil moisture from previous years alleviate strong declines of ANPP (Sala et al., 2012), which is usually treated as a time-lag effect (Petrie et al., 2018; Wu et al., 2015). Additionally, rain use efficiency also increases with water scarcity, meaning that less water is lost through runoff (Gutschick and BassiriRad, 2003; Huxman et al., 2004). (2) In wet years, other resources such as nutrient availability may increase with increasing precipitation, contributing to a supplementary increase of ANPP (Knapp et al., 2017b; Seastedt and Knapp, 1993). In contrast, the negative asymmetry index derived from observations at the moist STU suggests that this process is not dominant for this site, while temperature and/or light limitations that are associated with rainy periods may become important during wet years and neutralize the effect of increased precipitation on ANPP (Fig. S4) (Nemani et al., 2003; Wu et al., 2015; Wohlfahrt et al., 2008).

In our results, most models did not capture the sign of observed asymmetry indices across the three sites (Fig. 2c), which suggests that some of the underlying processes (combined carbon–nutrient interactions, time-lag effects, dynamic root growth allowing variation in accessible soil water) are not accurately represented in the models. For example, grassland root depth affects ecosystem resilience to environmental stress such as drought, and arid and semi-arid grasses that have extensive lateral roots or possibly deep roots show relatively strong resistance (Fan et al., 2017). However, most models currently consider only two types of grasslands – C₃ and C₄ (Table S14), with fixed root fractions in each prescribed soil layers (Table S13). This is potentially unrealistic for semi-arid grass roots and can lead to underestimating the amount of soil water available to plants and their resistance to drought. The latter is a key candidate especially for explaining the negative asymmetry index at the dry SGS.

The sensitivity of productivity to increased and decreased precipitation for simulations where mean precipitation was normally altered generally suggested negative asymmetric responses at dry (SGS) and mesic (KNZ) sites (Fig. 3c). This contrasts with a meta-analysis of grassland precipitation manipulation experiments (Wilcox et al., 2017) and with the ANPP–P conceptual model (Knapp et al., 2017b), which suggest a positive asymmetry response in the range of normal rainfall variation. This emphasizes the finding that most models overestimate drought effects and/or underestimate wet year impacts on primary productivity of dry and mesic sites for current precipitation variability. Under extreme conditions with modified precipitation, models were in line with the hypothesis and the data showing that ANPP saturates in very wet conditions but declines strongly in very dry conditions (Knapp et al., 2017b). For BNPP sensitivities to altered precipitation, meta-analysis of previous experiments indicated symmetric responses to increasing and decreasing rainfall (Luo et al., 2017; Wilcox et al., 2017), which may be regulated by allocation controls on the ratio of ANPP and BNPP to total NPP in response to altered precipitation. However, in the participating models, BNPP shows a negative asymmetric response to altered rainfall (Fig. 3d), which may reflect a shortcoming of carbon–water interactions in the belowground ecosystems.

4.2 Curvilinear responses of productivities to altered precipitation by models

In general, precipitation in ecosystem models is distributed through three pathways (N. G. Smith et al., 2014): (1) intercepted by vegetation and subsequently evaporated or falling on the ground; (2) infiltrated into the upper soil layers with subsequent evaporation, root water uptake and plant transpiration, or percolated down to deeper layers to form ground water; (3) runoff from the soil surface if the intensity of precipitation exceeds infiltration rates. In reality as well as in models, soil moisture rather than precipitation is the variable regulating vegetation growth, and biological responses to changes in precipitation are manifested as functions of soil moisture in different soil layers (Sitch et al., 2003; N. G. Smith et al., 2014; Vicca et al., 2012). We calculated the surface soil water content (SSWC, 0–20 cm depth converted from reported soil layers) and total soil water content (TSWC) under ambient and altered precipitation as simulated by the 14 models, and we found different patterns with parabolic, asymptotic and threshold-like nonlinear curves, which is similar to the response curves of primary productivity at the three sites (Figs. S5, S6). For the moist STU, SSWC and TWSC did not show obvious changes in response to increased precipitation since soil moisture at this site is often relatively near field capacity, while the SSWC and TSWC quickly decreased with decreasing in precipitation (Figs. S5, S6). In contrast, SSWC and TSWC at SGS showed significant increases in response to altered increased precipitation and slow decreases for decreased precipitation, because the soil was already very dry under average ambient conditions. Thus, changes of SWC in response to precipitation contribute to driving the different response patterns of simulated primary productivity across the grassland sites.

The responses of primary productivity to precipitation in models might also be driven by the intrinsic structure and parameterizations of vegetation functioning besides changes of soil moisture (Gerten et al., 2008), which account for the large spread in the values of b and a among models at the three sites (Figs. 4, 5, S3, S4). For example, carbon–nitrogen cycle coupling in ecosystem models reduced the simulated vegetation productivity relative to a carbon-only counterpart model (Thornton et al., 2007; Zaehle et al., 2010). Of those models used in this study, only five of the 14 models include carbon–nitrogen–water interactions (Tables 2, S1, S2). We calculated the ensemble mean of productivity for this group of carbon–nitrogen models (CLM45-ORNL, DLEM, DOS-TEM, LPJ-GUESS and TRIPLEX-GHG) and carbon-only models (CABLE, JSBACH, JULES, LPJmL-V3.5, ORCHIDEE-2, ORCHIDEE-11, T&C, TECO and VISIT) across altered and ambient precipitation simulations at the three sites, and then fitted the productivity–P responses with Eq. (8) (Figs. S7, S8, S9). We found that ensemble mean of carbon–nitrogen models generally produce a weaker GPP, NPP and ANPP response to precipitation than ensemble mean of carbon-only models and similar responses for BNPP. The latter may be explained by fixed root profiles in most models (Table S13). Our findings suggest that N interactions in ecosystem models reduced the productivity–P sensitivities, but should be confirmed using the same model prescribed with different N availability. In addition to the influence of nutrient cycling, different definitions of vegetation compositions (C₃/C₄) (Table S14), root profiles (Table S13), phenology (Table S9) and carbon allocation (Table S4) at the three sites may also contribute to the large variations of modeled productivity–P responses and demands for more accurate calibration of models to the specificity of the local sites in future model intercomparison studies.

4.3 Uncertainties, knowledge gaps and suggestions of further work

In this work, we applied two indices to characterize the asymmetry responses in the normal precipitation range using inter-annual variability of present conditions and forcing models with continuously modified precipitation amounts. Asymmetry indices from the inter-annual gross and net primary productivities suggest large uncertainties (Fig. 2), while the sensitivity analysis to changes in mean precipitation reported clear responses (Fig. 3). This can be explained by the differences in other climatic factors (for example, temperature, radiation and vapor pressure), or timing and frequency of precipitation between dry and wet years. All these uncontrolled factors may contribute to the large uncertainties of asymmetric responses from inter-annual variations (Chou et al., 2008; Peng et al., 2013; Robertson et al., 2009).

Although the carbon–water interactions in current models have been improved during the last decades, there still exist large gaps for accurately diagnosing the errors in the representation of key processes and parameterizations. Suggestions that should be considered in future studies aimed at model–data interaction include the following. (1) Models should report SWC at the same depth of experiments and experimental data should be made available for better comparisons in following studies. This can provide insights into the bias of modeled sensitivities to precipitation and check explicitly the sensitivity of vegetation productivity to change in SWC. (2) More experiments are needed that assess also BNPP in order to evaluate the corresponding processes in models (Luo et al., 2017; Wilcox et al., 2017). (3) There still exist large gaps between changes of precipitation occurrence and intensity in reality and how we simulated them in the current work, i.e., the altered rainfall forcing datasets were constructed by decreasing or increasing the amount of precipitation in each precipitation event by a fixed percentage during the time span of productivity observations at each site and not by modifying precipitation structure or reproducing the real treatment. Further studies need to consider better different scenarios of precipitation occurrence and intensity under climate change (Lauenroth and Bradford, 2012), which will likely help to better understand the responses of productivities to altered precipitation in the next decades. In addition, modelers will need to simulate the control experiments corresponding to the real local precipitation manipulations applied by field scientists, e.g., considering the observed time series of modified precipitation and vegetation composition, root profiles, nutrient cycling, phenology and carbon allocation as close as possible to local conditions. This should be a priority for future model–experiment interaction studies.

5 Conclusions

This is the first study where a large group of modelers simulated the response of grassland primary productivity to precipitation using long-term observations for evaluating the asymmetry responses to altered precipitation. Our results demonstrated that the multi-model ensemble mean captured the key observation of spatial slopes steeper than temporal slopes for ANPP. On the other hand, our analyses revealed that most models do not capture the observed positive asymmetry responses for the dry (SGS) and mesic (KNZ) sites under the normal precipitation conditions, suggesting an overestimation of the drought effects and/or underestimation of the watering impacts on primary productivity in the normal state. In general, current models represented a constant asymmetry pattern (negative asymmetry under normal and extreme conditions) across the full range of altered precipitation rather than a double asymmetry pattern (positive asymmetry under normal condition and negative asymmetry under extreme condition) established by Knapp et al. (2017b).

This study paves the path for further analyses where collaboration between modelers and site investigators needs to be strengthened such that also data other than ANPP can be considered and to identify which specific processes in ecosystem models are responsible for the observed discrepancies. This will eventually allow us to produce more reliable carbon-climate projections when facing different precipitation patterns in the future.