2.1 Simulation 1: repro-cox-1998

Our first simulation, repro-cox-1998, closely reproduces the optimal configuration presented in the Cox et al. (1998) study. Cox et al. (1998) modelled the fluxes for FIFE site 4439 (situated at 39^∘03^′ N, 96^∘32^′ W; 445 m above mean sea level). This tallgrass prairie site is roughly central within the 15 km × 15 km FIFE study area. It had been lightly grazed by domestic livestock but was ungrazed in 1986 and 1987 and was burned on 16 April 1987 (Kim and Verma, 1990a, 1991b). At the flowering stage in 1987, more than 80 % of the vegetation was composed of C₄ grasses (Kim and Verma, 1990a).

For their analysis, Cox et al. (1998) selected daylight hours that were both after 10:00 LT, to exclude dew evaporation, and from days with no rainfall during that day or the preceding day. This minimised the effect of evaporation of rainfall from the canopy and soil surface and let them focus on modelling transpiration and net canopy assimilation. We will also restrict our analysis to these same time periods. The model was spun up by repeating the entire run 10 times, and the output from the 11th run was analysed.

For driving data, we use a site-averaged product of the FIFE portable Automatic Meteorological Station (AMS) data at 30 min resolution (Betts and Ball, 1998). We prescribe both LAI and soil moisture from observations (Stewart and Verma, 1992) rather than calculating these variables internally using the JULES phenology or soil hydrology schemes. We use a “bucket approach” to calculate the soil moisture stress factor from the average soil moisture in the top 1.4 m (this option has been available from JULES 4.6 onwards), again to mimic the Cox et al. (1998) analysis. The wilting soil moisture θ_wilt was set to 0.205 m³ m⁻³ and the critical soil moisture θ_crit was set to 0.387 m³ m⁻³, taken directly from Cox et al. (1998). The resulting stress factor is plotted in Fig. 2 and clearly shows the dry period during late July and early August.

Figure 2Daily mean soil moisture stress factor β for each JULES simulation at FIFE site 4439 in 1987.

Download

JULES and the Cox et al. (1998) optimal configuration both use the Collatz et al. (1992) C₄ photosynthesis scheme. They also both use the same stomatal conductance parameterisation: Jacobs (1994), which is in turn a simplified version of the Leuning (1995) scheme. We select the “big leaf” option from the available canopy schemes in JULES, again to mimic Cox et al. (1998).

In this way, we are able to closely reproduce the Cox et al. (1998) calculation of daytime net canopy carbon assimilation and daytime canopy conductance with a modern version of JULES. Any remaining differences are minor. For example, in Cox et al. (1998), leaf temperature is calculated from the air temperature and observed sensible heat flux, whereas, in JULES, the full energy balance is modelled. There are also differences in the calculation of evaporation from soil and canopy, which are not the focus of this study. The calculation of aerodynamic resistance also differs. For example, in this run, canopy height is prescribed using the data from Verma et al. (1992) for this site in 1987 (see Sect. A5 for more information), whereas it was not modelled explicitly as part of the Cox et al. (1998) analysis.

Many of the key FIFE datasets used in this run have large uncertainties, despite being comprehensively measured by multiple teams. LAI measurements have an error of approximately 75 % due to the inherent variability of prairie vegetation. LAI measurements are also affected by leaf curling or folding as the leaves pass through the detector. There are therefore significant differences between datasets (for a more detailed description, see Sect. A2). For example, at the beginning of August, LAI measurements vary from 2.5 (Stewart and Verma, 1992) to 0.7 (the FIFE_VEG_BIOP_135 dataset). Soil moisture was also comprehensively measured across the FIFE area by multiple groups (see Sect. A3). While these observations are qualitatively consistent, one of the datasets shows a bias in the lower soil levels at site 4439 in 1987 compared to the other datasets. Within-site variability in soil moisture is also large. Soil properties were similarly well studied: there are four different datasets which can be used to calculate the wilting and critical soil moisture values, plus the values from two additional published studies (described in Sect. A4). However, measurements differ from each other by more than 0.15 m³ m⁻³ in some cases. There also appear to be differences between layers, with the top 10 cm having consistently lower wilting and critical thresholds than soil at a depth of about 30 cm, for example. It is therefore vital that we consider the implications of the spread in observed LAI, soil moisture and soil properties at this site when drawing our conclusions.

2.2 Simulation 2: global-C4-grass

In our second simulation, we use a recent JULES configuration, presented in Harper et al. (2016). This study introduced a trait-based approach to calculating leaf physiology in JULES and tuned plant parameters to observations in the TRY database (Kattge et al., 2011). Global vegetation was split into nine plant functional types (PFTs), including one to represent all C₄ grasses. The developments introduced in Harper et al. (2016) resulted in improved site-scale and global simulations of plant productivity and global vegetation distributions (Harper et al., 2018). Our global-C4-grass configuration is based on the representation of C₄ grasses in Harper et al. (2016) and takes advantage of many of the modern features of JULES. This includes a layered canopy scheme that treats the direct and diffuse components of the incident radiation separately (as in Sellers, 1985) and includes sunflecks (Dai et al., 2004; Mercado et al., 2007, 2009). It also calculates the overall soil moisture stress factor β from the sum of the stress factors in each layer, weighted by the root mass distribution. Since we are focusing specifically on the parameterisation of water stress, we continue to prescribe LAI and soil moisture, rather than calculate these parameters dynamically with the JULES phenology and soil hydrology schemes.

The driving data were taken from the site-averaged Betts and Ball (1998) product. The diffuse radiation fraction was calculated from shortwave radiation using the method in Weiss and Norman (1985) (see Sect. A1 for more information). A spherical leaf angle distribution was used, as in Harper et al. (2016). LAI was prescribed using the Stewart and Verma (1992) observations and the vegetation was set to generic C₄ grass.

The Stewart and Verma (1992) soil moisture observations were partitioned into the four JULES soil layers (thicknesses of 0.1, 0.25, 0.75 and 2.0 m) using an offline version of the soil hydrology scheme in JULES, assuming the same root distribution as natural C₄ grass in Harper et al. (2016). This is described in more detail in Sect. A3.1. The wilting and critical volumetric soil moisture values and the soil albedo were set to the same values as the repro-cox-1998 run. As Fig. 2 shows, the resulting soil moisture stress factor is almost identical to the simulation repro-cox-1998. Canopy height was also prescribed using the same observations as the repro-cox-1998 configuration, and the run was initialised from the spun-up repro-cox-1998 run.

2.3 Simulation 3: tune-leaf

For the third configuration, tune-leaf, we calibrate the JULES parameters to measurements of the tallgrass prairie vegetation at this particular site. At the flowering stage in 1987, the vegetation at FIFE site 4439 was dominated by three C₄ grass species: 27.1 % Andropogon gerardii (big bluestem), 22.2 % Sorghastrum nutans (Indiangrass) and 16.6  % Panicum virgatum (switchgrass) (Kim and Verma, 1990a). Since individual LAI observations for each species (as used in, e.g. Kim and Verma, 1991b) were not available, we continue to model this site with a single plant tile. We tune the leaf parameters of this tile to be approximately representative of the dominant species at this site, A. gerardii.

2.3.1 Leaf properties prior to the application of water stress in the model

As discussed above, JULES uses the Collatz et al. (1992) C₄ photosynthesis scheme to calculate the unstressed net leaf photosynthetic carbon uptake and the Jacobs (1994) relation to calculate stomatal conductance. In this section, we calibrate these parameterisations to the available in situ observations. A brief description of each of the model parameters fitted in this section is given in Table A2, and they are defined in full in Clark et al. (2011) and Best et al. (2011). Throughout this calibration work, the model points/lines are calculated with the Leaf Simulator package (Williams et al., 2019). This package exactly reproduces the way that JULES calculates leaf carbon uptake and stomatal conductance but allows leaf-level observations to be used as input.

Knapp (1985) compared leaf-level measurements of A. gerardii and P. virgatum in burned and unburned ungrazed plots on the Konza Prairie Research Natural Area in 1983 and the response of these two species to different water stress conditions. Their plots were located at 39^∘05^′ N, 96^∘35^′ W, which is within what subsequently became the FIFE study area. The burning occurred in April 1983, prior to initiation of growth of the warm-season grasses. They found significant differences between vegetation in the burned plot and unburned plots during the May–September period. The particular FIFE site we are modelling in our simulations, site 4439, was also burned prior to the start of the experiment (15 April 1987; Kim and Verma, 1990a) and was ungrazed throughout the FIFE period. Therefore, we use the observations from the burned plot in Knapp (1985) during May–June 1983, when they describe water availability as “not limiting” (we will investigate this claim in more detail in Sect. 2.3.2), to constrain our unstressed leaf photosynthesis parameters in the tune-leaf configuration. First, we set specific leaf area and the ratio of leaf nitrogen to leaf dry mass for A. gerardii and P. virgatum to Knapp (1985) observations taken between 25 May and 10 June 1983. Once these parameters are fixed, we then fit the other parameters in the model light response curve by comparison with the light curve presented in Knapp (1985), which was compiled from observations taken May–June 1983 at 35±2 ^∘C (Fig. 3).

Figure 3Mean observations from Fig. 1 in Knapp (1985) from the burned plot, early season (May–June 1983) for A. gerardii (cyan diagonal crosses) and P. virgatum (yellow vertical crosses) for net CO₂ assimilation rate against incident PAR, at 35±2 ^∘C. JULES parameters are fitted to the A. gerardii observations (cyan dashed line), P. virgatum (yellow dashed line) and a combination of both (green solid line). Also shown are the relations from the repro-cox-1998 (red dotted line) and global-C4-grass runs (blue dot-dashed line) at 35 ^∘C. Fitted lines assume no water stress (i.e. β=1) and c_i=200 µ mol CO₂ (mol air)⁻¹. Model lines have been created using the Leaf Simulator package, which reproduces the internal JULES calculations.

Download

Knapp (1985) also investigated the temperature dependence of net leaf photosynthesis by artificially altering the temperature of leaves of A. gerardii and P. virgatum. Their observations showed that the peaks in both species occurred at approximately the same temperatures but that the peak was significantly broader in A. gerardii than P. virgatum. In JULES, the temperature dependence of net leaf assimilation for C₄ plants is introduced through a temperature-dependent parameterisation of the maximum rate of carboxylation of RuBisCO, V_cmax. This enters the calculation of both the gross rate of photosynthesis and the dark leaf respiration R_d (since model R_d is proportional to model V_cmax). Therefore, we can use the relation between net leaf assimilation and temperature presented in Knapp (1985) to calibrate the JULES parameters governing the temperature dependence of V_cmax in the model. The result is illustrated in Fig. 4, alongside the parameterisations used in the repro-cox-1998 and global-C4-grass runs. The lines calibrated to the Knapp (1985) observations peak at approximately 38 ^∘C, whereas the repro-cox-1998 and global-C4-grass parameterisations peak at approximately 32 and 41 ^∘C, respectively. This leads to very different model behaviour in the temperature range 32–42 ^∘C, where the repro-cox-1998 parameterisation shows a dramatic decline in V_cmax, which contrasts sharply with the increase shown in the global-C4-grass parameterisation and the more stable lines calibrated to the Knapp (1985) observations. Note also that Polley et al. (1992) found “no apparent relationship” between leaf temperature and net leaf carbon assimilation in measurements of A. gerardii, S. nutans and P. virgatum, taken at ambient temperatures between 24.1 and 47.8 ^∘C. They speculate that the difference between their results and the temperature relations found by Knapp (1985) is due to seasonal acclimatisation.

Figure 4V_cmax against leaf temperature for A. gerardii (cyan diagonal crosses) and P. virgatum (yellow vertical crosses), using the normalised observations from Fig. 2 in Knapp (1985), scaled using the fitted light response curves of A. gerardii and P. virgatum at 35 ^∘C shown in Fig. 3. JULES parameters are fitted to these derived A. gerardii observations (cyan dashed line) and P. virgatum observations (yellow dashed line) and a combination of both (green solid line). Also shown are the relations from the repro-cox-1998 (red dotted line) and global-C4-grass runs (blue dot-dashed line). Model lines have been created using the Leaf Simulator package.

Download

As already stated, for the tune-leaf configuration, we use JULES parameters fit to the A. gerardii data from Knapp (1985), since A. gerardii is the dominant species at this site. However, to investigate the uncertainty introduced by the variation between species, we repeat the runs using parameters fitted to the approximate midpoint of A. gerardii and P. virgatum light response curves and V_cmax temperature relations. We would expect the best parameter set to lie between these two parameterisations. However, note that Knapp (1985) does not have data for Sorghastrum nutans, the second-most dominant plant species at FIFE site 4439, so we were not able to take this species into account in this part of the calibration.

It should also be noted that Knapp (1985) reported a drop in the ratio of leaf nitrogen to leaf dry mass over the course of the 1982 season of more than 50 % in the burned plots. This could be a contributing factor to the drop in leaf assimilation they observed over the course of 1983. We were not able to incorporate a time-varying ratio of leaf nitrogen to leaf dry mass into our simulations, which could lead to an overestimation of leaf assimilation in the senescence period.

There were also gas exchange measurements on individual leaves of A. gerardii, S. nutans and P. virgatum taken as part of the FIFE intensive field campaigns in 1987 (Polley et al., 1992). These observations were taken on upper canopy leaves perpendicular to the direct beam of the Sun, with varying absorbed PAR and internal CO₂ concentrations (FIFE_PHO_LEAF_46). This includes observations taken before, during and after the dry spell. Therefore, if we are to use these observations to calibrate the unstressed model parameters, we have to process them in such as way as to minimise the influence of the parameterisation of water stress in the model.

Figure 5Black crosses: A_l–c_i curves for Andropogon gerardii (left) and Panicum virgatum (right) from FIFE_PHO_LEAF_46 (Polley et al., 1992), normalised by the mean A_l of the data points with c_i>150 µmol CO₂ (mol air)⁻¹ in that curve. Only curves with mean incident PAR greater than 1200 µmol PAR m⁻² s⁻¹ have been used. Coloured points: normalised A_l calculated from observed c_i and incident PAR for each data point in the curve and the mean T_leaf observation for each curve, using the JULES relations. The JULES parameters are taken from the repro-cox-1998 configuration (red triangles), the global-C4-grass configuration (blue circles) and fits to A. g. data (tune-leaf default configuration; cyan diamonds) and P. v. data (yellow diamonds). Model points have been calculated using the Leaf Simulator package.

Download

To achieve this, we identified individual net leaf assimilation (A_l) versus leaf internal CO₂ concentration (c_i) curves from the FIFE_PHO_LEAF_46 dataset for A. gerardii and P. virgatum (using the observation time and leaf area). We normalised each A_l–c_i curve using the mean A_l for c_i>150 µmol CO₂ (mol air)⁻¹ for that curve. We then selected A_l–c_i curves with mean incident radiation greater than 1200 µmol PAR m⁻² s⁻¹. This procedure minimises the dependence on water stress or individual leaf nitrogen levels, since these factors approximately cancel out in the relations used internally in JULES when they are manipulated in this way. We can then use these normalised curves to calibrate the model A_l–c_i response at low c_i. For A. gerardii and, to a lesser extent, P. virgatum, this leads to a decrease in the initial slope of the A_l–c_i curve (Fig. 5).

We also attempted to use the A_l–c_i curves identified in the FIFE_PHO_LEAF_46 dataset to calibrate the parameters in the JULES relationship between internal leaf CO₂ concentration and external CO₂ concentration c_a. Each individual A_l–c_i curve was taken at approximately constant humidity, and c_a is also provided for each point on the curve. JULES uses the Jacobs (1994) parameterisation:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{c_{\mathrm{i}}-\mathrm{\Gamma }}{c_{\mathrm{a}}-\mathrm{\Gamma }}}=f_{\mathrm{0}}\left(\mathrm{1}-{\displaystyle \frac{dq}{d{q}_{\mathrm{crit}}}}\right),\end{array}$$

where Γ is the photorespiration compensation point (Γ=0 for C₄), and dq is specific humidity deficit at the leaf surface. f₀ and dq_crit are plant-dependent parameters: f₀ is a scaling factor on c_i and dq_crit governs the strength of humidity dependence of c_i. This parameterisation predicts that plotting c_i against c_a at constant humidity would give a straight line, with gradient $f_{\mathrm{0}}\left(\mathrm{1}-\frac{dq}{d{q}_{\mathrm{crit}}}\right)$. However, when plotting observations from FIFE_PHO_LEAF_46, we found that the slope of the c_i:c_a ratio changed as c_a increased (see Fig. S8 in the Supplement). Therefore, we were unable to calibrate the JULES c_i:c_a ratio to these data.

Figure 6Ratio of internal to external CO₂ against specific humidity deficit dq. Crosses are derived from leaf measurements of Andropogon gerardii (cyan) and other C₄ grasses (black), taken in the Konza Prairie (Jesse Nippert and Troy Ocheltree, published in Lin et al., 2015). Straight lines show the Jacobs model for C₄ plants, i.e. $c_{\mathrm{i}}:c_{\mathrm{a}}=f_{\mathrm{0}}\left(\mathrm{1}-\frac{dq}{d{q}_{\mathrm{crit}}}\right)$. Red dotted line: repro-cox-1998; blue dot-dashed line: global-C4-grass; green solid line: tune-leaf. Green dashed lines: varying dq_crit and setting f₀ to the best fit to the Lin et al. (2015) data for this dq_crit. Black dotted lines: Medlyn model using $g_{\mathrm{1}}^{\mathrm{fit}}$, $g_{\mathrm{1}}^{\mathrm{fit}}/\mathrm{2}$, $g_{\mathrm{1}}^{\mathrm{fit}}/\mathrm{4}$, where $g_{\mathrm{1}}^{\mathrm{fit}}$ is the value of the Medlyn model parameter g₁ fitted in Lin et al. (2015) to their Konza Prairie C₄ grass measurements. The green solid line (the tune-leaf configuration) is a good approximation to the Medlyn model with $g_{\mathrm{1}}=g_{\mathrm{1}}^{\mathrm{fit}}$ (because they have both been fit to the same dataset). The green dot-dashed and green dotted lines have been tuned to be close to the Medlyn model lines with $g_{\mathrm{1}}=g_{\mathrm{1}}^{\mathrm{fit}}/\mathrm{2}$ and $g_{\mathrm{1}}=g_{\mathrm{1}}^{\mathrm{fit}}/\mathrm{4}$, respectively.

Download

Instead, we use leaf measurements of C₄ grass in the Konza Prairie, collected in 2008 and published as part of Lin et al. (2015). These were taken at ambient CO₂ levels, under unstressed conditions. We can derive the c_i:c_a ratio from the supplied stomatal conductance, net assimilation and internal CO₂ observations, and plot this against specific humidity deficit at the leaf surface, calculated from chamber vapour pressure deficit (VPD), neglecting the effect of the leaf boundary layer (Fig. 6). We calibrate the Jacobs model parameters f₀ and dq_crit to these data (green solid line). Given the large scatter of the data and resulting poor fit (R²=0.04), we will also explore the effect of varying dq_crit (green dashed lines a, b, c). In each case, f₀ is set to best fit this dataset for this dq_crit (the parameter values are given in Table A3).

Both Knapp (1985) and Polley et al. (1992) found that leaf stomatal conductance g_s is proportional to the net leaf assimilation at this site. Their results are approximately consistent with the Lin et al. (2015) observations, given the difference in ambient CO₂ levels and the weak dependence on VPD.

As discussed above, in JULES, dark leaf respiration R_d is calculated from model V_cmax, scaled by a constant. For the tune-leaf simulation, we tune this constant such that the model dark leaf respiration at 30 ^∘C matches the dark leaf respiration from Polley et al. (1992) at 30 ^∘C (Fig. 7). This is roughly double the dark leaf respiration at 30 ^∘C in the repro-cox-1998 and global-C4-grass configurations. The Polley et al. (1992) relation was fitted to observations made at leaf temperatures of approximately 14–46 ^∘C. While our tuned model parameterisation of dark leaf respiration compares reasonably well in the range 25–35 ^∘C, it rapidly diverges from the Polley et al. (1992) observations beyond this range. This is particularly true for the higher temperature values, where the observations in Polley et al. (1992) show an increase with temperature, whereas the tune-leaf JULES configuration shows a decrease.

Figure 7Comparison of leaf dark respiration against leaf temperature relations from Polley et al. (1992) (black solid line) Kim and Verma (1991a) (black dotted line), repro-cox-1998 (red dotted line), global-C4-grass (blue dot-dashed line), tuned to A. g. (cyan dashed line), tuned to P. v. (yellow dashed line) and tuned to both A. g. and P. v. (green solid line). All lines assume no light inhibition of respiration. All JULES lines are top of the canopy (TOC) values without water stress. The lines that reproduce JULES configurations have been calculated using the Leaf Simulator package.

Download

Polley et al. (1992) found no significant difference between A. gerardii, S. nutans and P. virgatum for a variety of leaf properties: net leaf assimilation under ambient conditions, maximum assimilation under high light and CO₂ saturation, temperature response of net assimilation and relationship between assimilation and stomatal conductance under ambient conditions. This implies that the uncertainty we have introduced by not considering S. nutans data throughout most of this calibration is relatively minor.

2.3.2 Onset of water stress and relationship between water stress and leaf water potential

In this section, we calibrate the parameter governing the onset of soil water stress in the model, p₀. In the repro-cox-1998 and global-C4-grass simulations, p₀ is set to 1, meaning that the model vegetation starts to experience soil water stress at a volumetric soil moisture $\mathit{\theta }={\mathit{\theta }}_{\mathrm{crit}}=\mathrm{0.387}$ m³ m⁻³ (Fig. 1). This leads to a soil moisture stress factor β of 0.75–0.55 during the first 10 days of June 1987, i.e. a reduction of 25–45 % compared to the case where model vegetation is not limited by water availability (Fig. 2).

We can investigate this in more detail using leaf water potential observations as an indicator of the stress levels of the vegetation. Leaf water potential is affected by both the soil water content and the atmospheric water content, as well as other factors affecting transpiration. Both Polley et al. (1992) and Knapp (1985) found a relationship between leaf water potential and net leaf assimilation in their measurements of grasses in the FIFE study area. Polley et al. (1992) measured leaves of A. gerardii and S. nutans throughout the 1988 growing season. These observations showed a drop in net leaf carbon assimilation as the leaf water potential declined through the season: leaf water potentials from −0.34 to −1.5 MPa were consistent with net leaf carbon assimilation rates of 16.2 to 41.5 µmol m² s⁻¹, whereas lower leaf water potentials of −1.5 to −2.45 MPa were consistent with lower rates of 3.9 to 15.5 µmol m² s⁻¹ (at internal CO₂ concentrations of 200 µmol mol⁻¹ and absorbed PAR of 1600 µmol absorbed quanta m² s⁻¹). Knapp (1985) carried out weekly leaf water potential measurements of A. gerardii and P. virgatum in 1983 for late May to early October, which showed midday leaf water potential dropping from −0.4 MPa in late May to less than −6.6 MPa (the pressure chamber limit) at the end of July. During this period, net leaf assimilation dropped from approximately 40 µmol m² s⁻¹ to less than 10 µmol m² s⁻¹.

Kim and Verma (1991b) proposed a model which considers the prairie vegetation to be completely unstressed until the leaf water potential drops below −1 MPa. This was partially motivated by the Polley et al. (1992) measurements and evaluated using observations of FIFE site 4439 in 1987, i.e. the same site and time period we use in this study. Kim and Verma (1991a) proposed an alternative water stress model, also based on data in Polley et al. (1992), where both the maximum rate of carboxylation of RuBisCO V_cmax and the maximum rate of carboxylation allowed by electron transport J_max had a dependence on leaf water potential. According to this parameterisation, a leaf water potential of −0.4 MPa introduces a factor of 0.97 into V_cmax, for example, and a leaf water potential of −0.8 MPa introduces a factor of 0.91.

Midday leaf water potential for A. gerardii in the burned plot was approximately −0.4 MPa during the Knapp (1985) “early season” measurement period. Therefore, according to both the Kim and Verma (1991b) and Kim and Verma (1991a) models, considering this period “unstressed” is a very good approximation (i.e. β=1, to within 3 %) and agrees with their statement that “water was not limiting” the vegetation during this period. This validates our use of the Knapp (1985) dataset to tune the “unstressed” JULES parameters in the previous section.

Figure 8Leaf water potential observations for 4 days taken at FIFE site 4439 in 1987, published in Kim and Verma (1991a).

Download

We can now use the same arguments to determine how much water stress the vegetation should be experiencing at the beginning of June in our runs at FIFE site 4439 in 1987. Kim and Verma (1991a) present hourly leaf water potential measurements for A. gerardii leaves at this site, for a selection of days in 1987 (Fig. 8). On 5 June 1987, they measured a minimum leaf water potential of approximately −0.8 MPa at 14:00 LT. According to the Kim and Verma (1991b) model, vegetation at this leaf water potential would not be water stressed, and according to the Kim and Verma (1991a) model, V_cmax would be reduced by approximately 9 %. This contrasts sharply with the reduction in net assimilation throughout the day of 39 %, due to water stress (i.e. β=0.61), experienced in both the repro-cox-1998 and global-C4-grass simulations on this day.

For the tune-leaf configuration, we therefore reduce the early season water stress, to be more consistent with Kim and Verma (1991a) and Kim and Verma (1991b). This can be achieved by introducing a non-zero p₀ value in the stress factor β. This reduces the soil moisture threshold at which the plant becomes completely unstressed (β=1) from θ_crit to ${\mathit{\theta }}_{\mathrm{wilt}}+({\mathit{\theta }}_{\mathrm{crit}}-{\mathit{\theta }}_{\mathrm{wilt}})(\mathrm{1}-p_{\mathrm{0}})$, as illustrated in Fig. 1. Assuming that the stress factor β is 0.9 on 5 June 1987 leads to p₀=0.3. The effect of different values of p₀ will be shown in more detail in Sect. 3.

We now examine whether any previous modelling studies at this site support or conflict with this reduction in the soil moisture threshold at which the plant becomes completely unstressed. Crucially, the maximum soil moisture stress factor considered in the original Cox et al. (1998) study was 0.7; therefore, a setup with a p₀ of $\mathrm{1}-\mathrm{0.7}=\mathrm{0.3}$ and parameters re-tuned to give a 30 % reduction in unstressed net leaf assimilation would have given the same fit to the data. Similarly, a stress function with p₀=0.3 fits the plot of the ratio of actual to potential evapotranspiration to available water in Verma et al. (1992) (when corrected for their different soil properties) at least as well as a stress function with p₀=0. An increase in p₀ can also be considered a proxy for decreasing θ_crit (which, as we have already noted, has a large uncertainty; see Sect. A4). A p₀ of 0.2, for example, can be used to mimic the impact of changing θ_crit from 0.387, as used in this study and in Cox et al. (1998), to 0.348, as used in Verma et al. (1992).

Kim and Verma (1991a) present hourly water potential measurement of A. gerardii leaves at FIFE site 4439 for 3 other days (in addition to 5 June 1987): 2 July (peak growth period), 30 July (dry period) and 20 August 1987 (early senescence). These show a minimum of −1.2, −2.6 and −1.7 MPa, respectively (Fig. 8). Given the relationships between leaf water potential and net leaf assimilation described above, these leaf water potential measurements imply a drop in leaf assimilation during the middle of the day in the dry period. In contrast, Polley et al. (1992) found “no evident seasonal trend” in the maximum leaf assimilation rate or carboxylation efficiency, despite taking observations throughout the day before, during and after the dry spell in 1987¹. We were unable to reconcile these results satisfactorily using the associated data in the FIFE_PHO_LEAF_46 dataset (chamber vapour pressure, leaf and chamber CO₂ concentrations, leaf and chamber temperatures).

3.1 repro-cox-1998 and global-C4-grass  simulations

GPP in the repro-cox-1998 simulation after 10:00 LT compares very well to GPP derived from the flux tower data (Fig. 9), for all growth stages. This is expected, given that this simulation is designed to reproduce the model from Cox et al. (1998), which was tuned to this flux dataset. The global-C4-grass simulation reproduces the carbon fluxes reasonably well outside the dry period, although GPP is underestimated during the growth stages. For example, GPP is underestimated by approximately 30 % during the middle of the day on 5 June. During the dry period, however, the global-C4-grass simulation poorly captures the early morning peak and subsequent decline in GPP indicated by the carbon flux observations. The repro-cox-1998 run captures this behaviour through its response to leaf temperature. The diurnal cycle of air temperature on these days is shown in Fig. S5 and modelled leaf temperature in Fig. S6. Recall that V_cmax in the repro-cox-1998 simulation declines at leaf temperatures above 32 ^∘C. This causes a decline in modelled carbon assimilation during the hottest parts of the day (this is demonstrated explicitly in additional runs in the Supplement). However, as discussed in Sect. 2, the temperature response in the repro-cox-1998 configuration is not supported by observations in Knapp (1985) or Polley et al. (1992). Therefore, it appears that, while the model is successfully capturing the shape of diurnal cycle during the dry period, it is not achieving this with the correct physical process.

Similarly, net canopy assimilation in the repro-cox-1998 simulation compares well to the time series derived from the flux tower observations, although it has lower leaf respiration, particularly on 23 and 30 July (Fig. 10). As discussed in Sect. A7, the leaf respiration values assumed when processing the flux measurements were based on observations of leaf respiration in Polley et al. (1992). In Sect. 2.3, we showed that the repro-cox-1998 simulation underestimates leaf respiration compared to the Polley et al. (1992) dataset, particularly at the higher temperatures experienced during the middle of the day in the dry period. While the global-C4-grass configuration also simulates lower leaf respiration values than seen in Polley et al. (1992), a combination of a low bias in the GPP and a peak in V_cmax at higher temperatures (compared to the repro-cox-1998 simulation) reduces the impact on net canopy assimilation.

The latent heat flux is reasonably well modelled in general in both the repro-cox-1998 and global-C4-grass simulations outside the dry period (errors in the peak of the diurnal cycle of less than 20 %). However, both simulations overestimate the latent heat flux during the dry period (Fig. 11). This is expected, given that we have already shown that the canopy carbon assimilation is overestimated, and stomatal conductance is proportional to the net leaf assimilation in the model.

3.2 tune-leaf simulations

The tune-leaf configuration generally overestimates both GPP (Fig. 9) and net canopy assimilation (Fig. 10) compared to the observations and the repro-cox-1998 and global-C4-grass simulations. On days during the dry period, the tune-leaf simulation behaves characteristically similarly to the global-C4-grass simulation in that it also does not capture the midmorning peak and subsequent decline in GPP and assimilation. When fitting the tune-leaf configuration in Sect. 2, we highlighted uncertainties in some of the key parameters, and we will now look at the effect of these in turn.

Firstly, the tune-leaf configuration is based on observations of the dominant grass species at this site, A. gerardii. In Sect. 2, we also fitted parameters to another grass species at this site: P. virgatum and a “combined” set fitted to both species. Since A. gerardii is almost twice as abundant at this site in 1987 as P. virgatum, and in the absence of parameter fits to the other grass species at this site, we would estimate that the most representative parameters lie somewhere between these two parameter sets. Using this combined A. g./P. v. parameter set increases GPP and net canopy assimilation on the order of roughly 10 % compared to using the set fitted solely to A. g. (Figs. 9, 10), from which we conclude that the error introduced from using the dominant grass species is relatively minor.

Figure 12The diurnal cycle of GPP at site 4439 in the FIFE area for 8 days in 1987: 5 June (early growth), 2 and 11 July (peak growth), 23 July, 30 July and 11 August (dry period), and 17 and 20 August (early senescence). Green band shows how the tune-leaf simulation would vary for p₀ in the range 0 to 0.4 (lower limit corresponds to p₀=0, upper limit to p₀=0.4).

Download

A key difference between the tune-leaf configuration and the other configurations is the introduction of a non-zero p₀. Figure 12 illustrates that varying p₀ from 0 (as in the repro-cox-1998 and global-C4-grass simulations) to 0.4 has a strong effect on GPP, as expected. It demonstrates the importance of ensuring that the threshold for water stress is consistent with the “unstressed” leaf observations we calibrated against. Continuing to use p₀=0 with the newly tuned unstressed parameters would have resulted in much-too-low GPP during the early growth period. Recall also that changing p₀ can be considered a proxy for changing the critical soil moisture. Therefore, these runs also demonstrate the sensitivity to uncertainty in the soil properties.

Figure 13The diurnal cycle of GPP at site 4439 in the FIFE area for 8 days in 1987: 5 June (early growth), 2 and 11 July (peak growth), 23 July, 30 July and 11 August (dry period), and 17 and 20 August (early senescence). Green band shows how the tune-leaf simulation would vary for a canopy structure factor a in the range 0.3 to 1 (upper limit corresponds to a=1, lower limit to a=0.3).

Download

The effect of varying the canopy structure factor on GPP can be seen in Fig. 13. This can also be seen as a proxy for examining the effect of reducing LAI as it changes the effective LAI seen by the model radiation scheme. Varying the canopy structure factor in the range 0.8–1.0 has a negligible effect on GPP on these days. However, reducing the canopy structure factor from 0.8 to 0.3 has a large negative impact on GPP. As discussed in Sect. 2, this range is inside the error given in the LAI dataset documentation. The error in LAI for this site therefore has a large impact on the modelled canopy carbon fluxes.

Figure 14The diurnal cycle of GPP at site 4439 in the FIFE area for 4 days in during the dry period of 1987. Solid green lines use the tune-leaf configuration. Dashed dashed green lines show how GPP varies if dq_crit is increased, while f₀ is changed to maintain the best fit to the Konza Prairie C₄ grass observations in Lin et al. (2015) (upper, middle and lower dashed lines correspond to parameter combinations a, b and c, respectively, as defined in Table A3).

Download

Less straightforward to investigate is the effect of the uncertainty in the calibration of the JULES c_i humidity response. Recall that the observational dataset used in Sect. 2 had a large spread in c_i compared to its range of specific humidity deficit values. This made it difficult to tune the parameter dq_crit separately to the overall scaling factor f₀. We therefore take the approach of systematically varying dq_crit (while setting f₀ to keep the best fit to the observations in Fig. 6), to show qualitatively that a different humidity calibration cannot improve the agreement with the GPP observations during the dry spell. Figure 14 compares modelled GPP for three different dq_crit, f₀ combinations: dq_crit=0.048, f₀=0.59 (upper green dashed line), dq_crit=0.040, f₀=0.64 (central green dashed line) and dq_crit=0.035, f₀=0.68 (lower green dashed line) for 4 days during the dry spell. Plots of specific humidity deficit on these days are given in Fig. S7. None of these parameter combinations are able to fit the steady but low rate of GPP during the middle period of the day: they transition from almost no humidity-induced effect on GPP to a sudden decline. The timing of this decline varies across the 4 days shown. This demonstrates that, while lower c_i values in these runs during the day in the dry period can reduce GPP, the magnitude of the slope of c_i:c_a against dq is too large. These two effects cannot be reconciled while still maintaining consistency with the unstressed observations in Lin et al. (2015). This implies that the Jacobs parameterisation used in JULES, where the relationship between c_i:c_a and specific humidity deficit does not vary over the course of the run, does not have the flexibility needed to capture the behaviour of GPP at this site.

3.3 What potential model developments could improve the diurnal cycle of JULES GPP at this site?

As we have seen, the global-C4-grass configuration, which is typical of how this site would be modelled in a global JULES run, is unable to capture the diurnal cycle of GPP (and also net canopy assimilation and latent heat flux) at this site during the dry period in 1987. Replacing the generic C₄ grass tile parameters with parameters that are calibrated to observations taken of vegetation at this particular site (the tune-leaf configuration) does not improve ability of the model to capture the diurnal cycle in these fluxes. We have demonstrated that this conclusion is robust to uncertainties in LAI, soil moisture, leaf parameters, canopy parameters and soil parameters.

We will now explore a number of possible options for improving the standard representation of the dry period diurnal GPP cycle at this site. Firstly, the model diurnal cycle can be greatly improved via the careful selection of parameters in the existing leaf temperature-dependent calculation of V_cmax. This was demonstrated in the model runs in Cox et al. (1998), which we have closely reproduced with the repro-cox-1998 configuration. This method has the advantage that it provides a close fit to data and does not require any changes to the model code. A disadvantage of this method is that the V_cmax model parameterisation becomes an effective parameterisation which no longer has a clear biological interpretation. It therefore becomes more difficult to constrain from results in the literature. The numerical success of this method is due to high leaf temperatures acting as a proxy for high atmospheric demand during the middle of the day in the dry period (Figs. S6 and S7). While these temperature parameters provide a good approximation at this site in this particular year, it does not follow that these same temperature parameter values would be appropriate for other locations or at this location under a changing climate.

Secondly, the model could be extended to include a soil moisture effect on the internal leaf CO₂ concentration c_i. As we demonstrated in Sect. 3, the current expression for c_i in JULES cannot simultaneously fit the unstressed observations and be able to reduce c_i to the required levels to affect GPP during the dry season without also increasing the strength of the response to specific humidity deficit. This results in the humidity-induced stomatal closure occurring too suddenly on days during the dry period. Introducing a soil moisture dependence in c_i would allow c_i to be lower on days where soil water was limiting for all humidity levels, while maintaining the higher values on unstressed days. Zhou et al. (2013) and De Kauwe et al. (2015) both achieve this by adding a soil moisture dependence to the VPD term in the Medlyn conductance model (Medlyn et al., 2011). The Medlyn model is based on the theoretical argument that stomata should act to minimise the amount of water used per unit carbon gained, leading to a stomatal conductance $g_{\mathrm{0}}+\mathrm{1.6}\left(\mathrm{1}+\frac{g_{\mathrm{1}}}{\sqrt{D}}\right)\frac{A}{c_{\mathrm{a}}}$, where g₀ and g₁ are free parameters.

As demonstrated in De Kauwe et al. (2015), the parameters in the Jacobs model (f₀, dq_crit) can be chosen so that the resulting c_i:c_a ratio approximates the Medlyn model, for mid-range VPD values. The unstressed Konza Prairie C₄ grass measurements used in Sect. 2 to calibrate the c_i:c_a ratio in the tune-leaf configuration were actually provided in Lin et al. (2015) as part of a comprehensive study to tune the g₁ parameter in the Medlyn model for different vegetation types (with g₀=0). Using the Medlyn model with their calibrated g₁ value ($g_{\mathrm{1}}^{\mathrm{fit}}=\mathrm{1.04}$) does indeed give a similar c_i:c_a to our tune-leaf configuration (Fig. 6, solid green line).

Figure 15The diurnal cycle of GPP at site 4439 in the FIFE area for 4 days in during the dry period of 1987. Solid green lines use the tune-leaf configuration. The c_i:c_a ratio in this configuration closely corresponds to the c_i:c_a ratio for C₄ grasses in the Konza Prairie in Lin et al. (2015), using the Medlyn model and fitting the Medlyn model parameter g₁ to measurements taken in 2008 ($g_{\mathrm{1}}=g_{\mathrm{1}}^{\mathrm{fit}}=\mathrm{1.04}$ kPa^−0.5). The dot-dashed lines and dotted lines show the results from fitting the JULES parameters dq_crit and f₀ to approximate the Medlyn model when $g_{\mathrm{1}}=g_{\mathrm{1}}^{\mathrm{fit}}/\mathrm{2}$ and $g_{\mathrm{1}}=g_{\mathrm{1}}^{\mathrm{fit}}/\mathrm{4}$ (the parameter values are given in Table A3).

Download

Therefore, to investigate the effect of a soil moisture-dependent g₁ on GPP, we can set the JULES c_i:c_a ratio to mimic a lower g₁ and try this out on days with low soil moisture. For this test, we choose JULES parameter values that provide a rough approximation to the Medlyn model with $g_{\mathrm{1}}=g_{\mathrm{1}}^{\mathrm{fit}}/\mathrm{2}$ and $g_{\mathrm{1}}=g_{\mathrm{1}}^{\mathrm{fit}}/\mathrm{4}$ (Fig. 6, dot-dashed and dotted green lines). These reductions in g₁ are well within the range observed in Zhou et al. (2013) for a range of different vegetation types under water-limited conditions. The resulting JULES parameter values are given in Table A3. Figure 15 demonstrates that lowering c_i:c_a in this way is able to qualitatively reproduce the shape of the diurnal cycle of GPP in the dry period. The run mimicking $g_{\mathrm{1}}=g_{\mathrm{1}}^{\mathrm{fit}}/\mathrm{2}$, in particular, is a very good match to the observations. This shows the potential value of extending JULES to allow interaction between the plant response to soil moisture dependence and VPD.

Another way to implement this interaction in JULES would be to add a dependence on leaf water potential, since leaf water potential is affected by both soil moisture (water supply) and VPD (atmospheric water demand). As discussed in Sect. 2.3, there is an observed relationship between leaf water potential and leaf assimilation in grass species at this site.

Previous studies have demonstrated that models with an explicit dependence on leaf water potential can successfully capture the dry period diurnal cycle at this site. Kim and Verma (1991a) were able to qualitatively capture the midmorning peak and subsequent decline in net canopy photosynthesis on 30 July at this site, using a model in which both V_cmax and J_max had a dependence on their leaf water potential measurements. Furthermore, Kim and Verma (1991b) were able to reproduce this behaviour in canopy conductance at this site on 30 July and 11 August 1987 using a model that included an explicit dependence on observed leaf water potential, in addition to a direct dependence on VPD.

Leaf water potential is not currently modelled explicitly within JULES. Typically, in plant hydraulic models, leaf water potential is calculated assuming a steady-state water balance, using the soil water potential, transpiration, and leaf-to-root and root-to-soil resistance terms (as in, e.g. Newman, 1969). Adding this to the JULES code is technically non-trivial, as water stress is currently applied to leaf-level processes before transpiration is calculated. Also, modelling the plant resistances would require additional input parameters, which would need to be constrained from observations.

Stress parameterisations involving leaf water potential come in a range of complexities. The simplest involve inserting a leaf-water-potential-dependent stress factor into an existing part of the model, e.g. the limiting photosynthesis rates as in Kim and Verma (1991a), or stomatal conductance, as in Kim and Verma (1991b) and Tuzet et al. (2003). More sophisticated models include the plant hydraulics as part of schemes incorporating risk–benefit analysis (e.g. Sperry et al., 2017; Eller et al., 2018) and/or chemical signalling (e.g. Tardieu and Davies, 1992; Dewar, 2002; Huntingford et al., 2015).

Finally, another way to improve the diurnal cycle of GPP in the dry period would be to incorporate a parameterisation of leaf rolling. For example, effective leaf area available to the radiation scheme could be decreased during hot, dry weather. Kim and Verma (1991a) attribute the residual overestimation of net canopy carbon assimilation on days during the dry period of their leaf water potential-based model to this effect. It would therefore be interesting to investigate the contribution that leaf rolling makes to the overall plant water use strategy. However, while the occurrence of leaf rolling/folding at the FIFE site has been recorded, the effect has not been quantified. This would be a necessary first step for modelling this process at this site.

3.4 Can the FIFE dataset make a useful contribution to current-day JULES evaluation and development work?

A global land-surface model such as JULES needs to perform well for a wide range of climate regimes, timescales, spatial scales and vegetation types. Model evaluation or development work needs to represent this variety. The availability of comprehensive databases, such as FLUXNET (Baldocchi et al., 2001) and TRY (Kattge et al., 2011), have revolutionised land-surface science by giving easy access to observations from a wide variety of sources, in a common format. Given this context, why would a modeller consider also using the FIFE dataset?

Firstly, FIFE provides an ideal case study for improving the model representation of water stress on carbon and water fluxes in JULES in tallgrass prairie. While, at one time, tallgrass prairie extended over 10 % of the contiguous United States (Fierer et al., 2013), it has declined 82 %–99 % since the 1830s due to agricultural use (Sampson and Knopf, 1994; Blair et al., 2014a). However, grasslands in general (including other grass- and graminoid-dominated habitats, such as savanna, open and closed shrubland, tundra) cover more terrestrial area than any other single biome type (up to 40 % of Earth's land surface; Blair et al., 2014a). It is therefore important to include lots of examples of grasslands in any global analyses of vegetation responses to changing conditions. The Konza Prairie LTER site, where FIFE was based, has been used extensively to investigate the dynamics and trajectories of change in temperate grassland ecosystems, including drivers such as fire, grazing, climate and nutrient enrichment (see Blair et al., 2014b for a review).

FIFE looked at the processes for representing water stress in detail and intensively studied the relevant factors. This has led to a wide variety of complementary observations and literature specifically focusing on how these data can be used to inform models. LAI is a good illustration of this advantage. As we have discussed, LAI is an important parameter for modelling canopy water and carbon fluxes. LAI was measured by multiple groups at FIFE, directly and indirectly, and the large differences found between the different attempts was fully explored at the time. We can use their results to inform our own use of these datasets.

When adding a new process to a global land-surface model, it is important to tune new parameters to a comprehensive range of datasets. For example, as mentioned in Sect. 3.3, Lin et al. (2015) use data for 314 species from 56 sites across the world to tune the new g₁ parameter introduced in the Medlyn model of stomatal conductance for key plant functional types. This breadth of sites and vegetation types is essential. Each site contributed leaf gas exchange observations taken under similar protocols to allow a carefully controlled common analysis.

Access to individual experiments, which have investigated the combined effect of a wide range of processes, such as FIFE, can play a complementary role in land-surface evaluation and development. For example, FIFE provides cases where improving an individual process in isolation degrades overall model performance. As we have shown, calibrating unstressed model V_cmax (T=25 ^∘C) from leaf observations without also calibrating when the model is considering the vegetation to be unstressed significantly underestimates early season GPP. Similarly, tuning the model parameters to improve the fit to canopy GPP and evapotranspiration can result in an unrealistic temperature dependence of V_cmax. Looking at sites in a holistic way can also highlight complications or influences that might not a priori have been considered, such as leaf rolling in our case.

There are two main disadvantages to the use of FIFE in evaluation and model development studies. The first is the limited time period: observations are available for a period of up to 3 years, with some key measurements only undertaken during the intensive field campaigns. Where long-term effects are being studied, alternative datasets would need to be used.

The second disadvantage is that it is relatively more time consuming to add FIFE to an evaluation study, compared to adding an extra site from one of the large, standardised databases such as FLUXNET. This is partly because FIFE provides a choice of different datasets to use for forcing, calibrating parameters and evaluation, which takes time to investigate. It is also partly because, although the data are easily downloadable, well documented and in common file formats, they still need to be manipulated into a format that can be used in JULES runs. We aim to address this issue by providing a suite that can be used to preprocess the FIFE data and run JULES with the configurations described in this paper (see the “code and data availability” section).

This aim is central to the provision of this paper. FIFE is the first “JULES golden site”, a concept that was launched at the annual JULES meeting in 2018. A JULES golden site is a site targeted by the JULES community because it can help address one of the key science questions facing JULES and has high-quality observational data that can be used to drive JULES and evaluate the output. It creates a network of researchers within the JULES community with experience of how this site can be exploited for JULES development, with input from site investigators. A key component is the provision of shared runs and evaluation datasets, which can be gradually expanded and improved.

In our study, we have focused on the contribution that FIFE can make to the development of water stress in JULES. This has governed the choices we have made when setting up our configurations, e.g. choosing to prescribe LAI and soil moisture. However, we note here that FIFE could also be used to investigate other processes, such as plant and soil respiration (Sect. A7), the seasonal decline in leaf nitrogen (Knapp, 1985) and the modelled energy balance (Kim and Verma, 1990a; Colello et al., 1998).