Evaluation of JULES-crop performance against site observations of irrigated maize from Mead, Nebraska

. The JULES-crop model (Osborne et al., 2015) is a parametrisation of crops within the Joint UK Land Environment Simulator (JULES), which aims to simulate both the impact of weather and climate on crop productivity and the impact of croplands on weather and climate. In this evaluation paper, observations of maize at three FLUXNET sites in Nebraska (US-Ne1, US-Ne2 and US-Ne3) are used to test model assumptions and make appropriate input parameter choices. JULES runs are performed for the irrigated sites (US-Ne1 and US-Ne2) both with the crop model switched off (prescribing leaf area index (LAI) and canopy height) and with the crop model switched on. These are compared against GPP and carbon pool FLUXNET observations. We use the results to point to future priorities for model development and describe how our methodology can be adapted to set up model runs for other sites and crop varieties.


Introduction
The Joint UK Land Environment Simulator (JULES) Clark et al., 2011) is a process-based model that simulates the fluxes of carbon, water, energy and momentum between the land surface and the atmosphere. It is used in carbon cycle, climate change and impacts studies, and can be run on its own ("stand-alone" mode) or as a component of a coupled Earth system model. As described in the model description paper (Osborne et al., 2015), JULES-crop is a parametrisation of crops that has been added to JULES in order to improve land-atmosphere interactions in areas where crops are predominate in addition to enabling the simulation of the effect of weather and climate on food and water resources.
JULES treats each vegetation type as existing on a separate tile within a grid box. Energy and carbon flux calculations are performed separately for each tile, and prognostics, such as leaf area index (LAI) and canopy height, are calculated and stored for each tile separately. Each vegetation tile has a different set of input parameters and leaf-level carbon assimilation is calculated differently depending on whether the tile is modelling a plant with a C3 or a C4 plant photosynthetic pathway. JULES-crop introduces a distinction between natural plant functional types (PFTs) and crops. Crop tiles have their growth and development parametrised by a crop development index (DVI) and have different calculations for the allocation to plant carbon pools, leaf area index and height compared to natural PFTs. However, in most other respects, Published by Copernicus Publications on behalf of the European Geosciences Union.
such as the calculation of gross primary productivity (GPP) and respiration, natural PFTs and crops are modelled in the same way within the JULES code. In its current stage of implementation, JULES-crop is available only in offline JULES runs, although there are plans to extend it for use in coupled runs in the future.
Other land-surface models have also been extended include specific representations of key crops. For example, Community Land Model (CLM)-crop has been evaluated at the site level for several crop types (maize, soybean and spring wheat (Drewniak et al., 2013); winter wheat (Lu et al., 2016) and physiology parameters were calibrated to optimise productivity (Bilionis et al., 2015). ORCHIDEE-CROP has been evaluated for maize and winter wheat at a number of European sites (Wu et al., 2016) and was shown to reproduce the seasonality of leaf area index and carbon and energy fluxes. Similarly, the incorporation of a phenology scheme into the SImple Biosphere (SIB) model improved the prediction of both leaf area index and carbon fluxes for maize, soybean and wheat crops at a number of sites in North America (Lokupitiya et al., 2009). Song et al. (2013) implemented crop-specific phenology and carbon allocation schemes into the Integrated Science Assessment Model (ISAM) landsurface model and calibrated against observational data from a corn-soybean rotation at Mead and Bondville (US) sites. This model was able to reproduce the diurnal and seasonal variability of carbon, water and energy fluxes.
In Osborne et al. (2015), global runs using JULES-crop were carried out for four generic crop types -maize, soybean, wheat and rice -and the effect of including the new crop parametrisation was shown on sensible heat flux, moisture flux and net primary productivity (NPP) for some key countries. The model yield was also compared against global and country FAO crop yields. Site runs were performed at four FLUXNET sites with a maize-soybean rotation: Mead (US-Ne2 and US-Ne3), Bondville (US-Bo1) and Fermi (US-IB1). For input parameters that applied to both natural vegetation and crop tiles, C3 crops were given the parameter values of a standard C3 grass tile within JULES and C4 crops were given the values of a standard C4 grass tile. Osborne et al. (2015) speculated that an improved fit to observations could be obtained if these parameters were tuned to be more crop specific.
The other published study using JULES-crop to date, Williams and Falloon (2015), used the global set-up and the generic parametrisation of the four main crops from Osborne et al. (2015) to investigate the sensitivity of the yield from JULES-crop to the driving data variables, assessing both the relative importance of different variables and whether there is an advantage to using subdaily driving data rather than using daily driving data and performing an internal disaggregation to subdaily timescales. It also investigated the effect on the yield of initialising the model from climatology. No attempt was made to find more appropriate crop parameter values.
In this model evaluation paper, we use the observations available at the Mead FLUXNET sites US-Ne1, US-Ne2 and US-Ne3 to investigate how well each individual component of JULES performs for maize and how much of an improvement can be achieved by using more appropriate parameter values, taking into account advances in the JULES code since the Osborne et al. (2015) study. This investigation splits into three distinct parts. We initially look at which processes and parameters can be tuned directly to maize observations from the Mead sites, without running the model. Second, for parts of the code shared between natural PFTs and crops in the model (the calculation of gross primary productivity and respiration), we test the performance of the tuned parameters by running JULES with the crop model switched off and forcing with observed leaf area index (LAI) and canopy height, to remove the feedback between net primary productivity and LAI. Finally, we will use the tuned parameters in JULES runs for irrigated maize at Mead with the JULES-crop parametrisation switched on.
This paper is organised as follows. Section 2 gives information about the observations and the model set-up used for the JULES runs presented in this paper, both those with and without the JULES-crop parametrisation switched on. Particular attention is paid to the choice of input parameter values, which are tuned to the available observations. Section 3 compares the results from the model runs against the observations. Section 4 contains an overall assessment about the suitability of the model for modelling maize at these sites and discusses ways that the model could be improved. It also comments on the more general applicability of the parameters and methods used in this paper for tuning JULES for other sites and crop varieties. A summary of the JULES-crop parametrisation and the other relevant parts of the JULES code is given in Sect. A.
2 Experimental set-up

Observations
There are three FLUXNET sites at the University of Nebraska Agricultural Research and Development Center near Mead, Nebraska, which are located within 1.6 km of each other: US-Ne1, US-Ne2 and US-Ne3. Both US-Ne1 and US-Ne2 are irrigated with a central pivot system, whereas US-Ne3 is entirely rainfed Suyker et al., 2004Suyker et al., , 2005. US-Ne1 grows maize, whereas US-Ne2 and US-Ne3 are maize-soybean rotations. The observations span from 2001 to 2015 (although not all variables were available for this entire period).
The observations of the biomass of green leaves, yellow leaves, stem and reproductive parts of maize (kernel, cob, husk, ear shank, silk) were made after the plant material was dried to a constant temperature of 105 • C. In the observations, green leaves encompasses all green leaf material from Geosci. Model Dev., 10, 1291-1320, 2017 www.geosci-model-dev.net/10/1291/2017/ the collar to the leaf tip, yellow leaves are defined as greater than 50 % necrotic (or entirely yellow) leaf and the stem includes stem, leaf sheaths, immature or undeveloped ears and unfurled leaves. Hourly incident and absorbed Photosynthetically Active Radiation (PAR) (400 to 700 nm) observations are available from the Mead FLUXNET sites. Absorbed PAR was calculated using two point quantum sensors above the canopy, pointing up and down, and two line quantum sensors below the canopy, pointing up and down. The line quantum sensors below the canopy integrate over an area 1 cm by 1 m, in order to even out effects such as sunflecks.
The observations were used in three ways: to determine the input parameters to the JULES runs (air temperature, carbon pools, leaf nitrogen, absorbed PAR, canopy height, LAI), to drive the JULES runs themselves (meteorological variables, LAI, canopy height) and to compare the JULES run results against (GPP, carbon pools, LAI, canopy height). Observations from all three sites were considered in the input parameter tuning, whereas only observations from the irrigated sites were used to drive and validate the model runs.

Model set-up
The relevant features of the JULES land-surface model, including the JULES-crop parametrisation, are described in Sect. A. Two types of JULES runs were used in this study: 1. Maize is treated as a natural PFT tile (i.e. crop model is switched off), with LAI and crop height prescribed from observations (linearly interpolated to create a daily time series).
2. Maize is considered as a crop tile (i.e. crop model is switched on).
The runs were driven by hourly observations of downward shortwave radiation, downward longwave radiation, precipitation, air temperature, wind speed, pressure, specific humidity and diffuse radiation fraction. Each year and site was modelled as a separate run, each starting on 1 March. Annual global CO 2 atmospheric concentrations were taken from Dlugokencky and Tans (2016).
The following sections describe in more detail how the choice of input parameters was made. Observations from both the irrigated sites at Mead and the rainfed site at Mead were considered when tuning the model input parameters that were designed to take the same value whether irrigation is switched on or off in the model. However, in these cases, observations from the rainfed site are clearly denoted on the plots, in order to check for cases where these model approximations break down. It was assumed that there was no limitation from nitrogen availability. A summary of the model input parameters used in both types of runs are given in Tables 1, 2a, 2b, 3a, 3b and 4.

Crop development parameters
The cardinal temperatures T b , T o and T m in this analysis have been kept the same as Osborne et al. (2015), which were chosen based on the literature review in Sánchez et al. (2014). As in Osborne et al. (2015), there was the assumption of no dependence of thermal time on the photoperiod.
The thermal times were calculated using the available Mead data for the sowing date, the date at which 50 % of the plants had emerged 1 , the date at which 50 % of the plants were at the R1 or "estimated R1" growth stage (i.e. had begun the reproductive phase), the date at which 50 % of the plants had reached the R6 growth stage (maturity) and the harvest date, together with the observed hourly air temperature and Eq. (A1). These thermal times are given in Table 5. In the runs presented in Sect. 3, the thermal times for sowing to emergence, emergence to flowering and flowering to harvest for each year at a site are used in JULES-crop directly, to simulate the crop development as closely as possible for a finished crop season, where the harvest date is known. 2 The sowing date is prescribed (i.e. l_prescsow=T). An option for sowing date to be calculated dynamically using rate of change of day length and soil temperature and moisture does exist (l_prescsow=F), but this is not considered here as it is still under development and not recommended for use (Osborne et al., 2015).
Since harvest dates are available, T mort was set low enough that it did not trigger harvest.

Carbon partitioning
The carbon partitioning parameters α i , β i were tuned to observations of the biomass of green leaves, yellow leaves, stem and reproductive parts of maize. The ratio of carbon to biomass in each part of the plant was assumed to be the same and constant in time. The C leaf pool in the model contains green leaves only (since C leaf is directly linked to LAI and photosynthesis) and the C harv pool consists of both the reproductive parts of the plants and the yellow leaves. Stem carbon in the model is split between the C stem and C resv pools. The biomass observations were linearly interpolated to get a daily time series and then differentiated with respect  Figure 1. Top: ratio of rate of change of C leaf to rate of change of above-ground carbon C ag ; middle: ratio of rate of change of C stem + C resv to rate of change of above-ground carbon; bottom: ratio of rate of change of C harv to rate of change of above-ground carbon. Solid black line uses the original crop parameters from Osborne et al. (2015), dashed black line uses the tuned parameters. Blue, green and red lines are derived from US-Ne1, US-Ne2 and US-Ne3 observations respectively.
to time. Ratios of these rates were then plotted as a function of DVI (Fig. 1). Using these plots alongside the function for root carbon from de Vries et al. (1989) (since there were no direct measurements of root biomass available from the Mead sites), new, tuned values for α i , β i were found. These tuned parameters (dashed lines) show an improvement in the proportion of the increase in above-ground carbon that goes to the green leaves ( Fig. 1, top) and the proportion of the increase in above-ground carbon that goes to the stem ( Fig. 1, middle) for DVI < 0.8 as compared to the parameters used in Osborne et al. (2015) (solid line). However, note that, even after the tuning, the proportion of the increase in aboveground carbon that goes to the green leaves does not drop off sharply enough for DVI > 0.8 compared to the observations. The tuned partition fractions are shown more clearly in Fig. 2 (colours), together with the functions given in de Vries et al. (1989) (the α i , β i in Osborne et al., 2015 were fitted to these functions with minor adjustments as a result of global runs). It was not possible to fit p root accurately to the expression from de Vries et al. (1989) for approximately a DVI of 1.0 to 1.4 given the constraints above. In addition, in reality, water stress can also increase the fraction of NPP going to the roots (see discussion in, e.g., de Vries et al., 1989 andSong et al., 2013), but this effect is not taken into account in JULEScrop. However, we do not see a notable difference between the irrigated sites US-Ne1 and US-Ne2 (blue and green lines respectively) and rainfed site US-Ne3 (red lines) in Fig. 1.

Remobilisation of stem carbon
The stem biomass observations were used to tune the value for the stem reserve remobilisation constant τ . The relation governing the stem reserve remobilisation can be rearranged to where M stem is stem biomass (including reserves), M max stem is the maximum value of M stem in that site in that year and d max is the day since M max stem occurred.
Geosci. Model Dev., 10, 1291-1320, 2017 www.geosci-model-dev.net/10/1291/2017/  Therefore, plotting 1 − M stem M max stem against 1 − 0.9 d max should give a straight line with gradient τ . Using the assumption that the day with maximum stem biomass was approximately the same day as the day with the maximum stem biomass measurement, a straight line was fitted to the observations and an approximate value of τ = 0.12 was obtained. However, as can be seen in Fig. 3 (which displays both the new, tuned value τ = 0.12 (black, dashed line) and the value used in Osborne et al., 2015 of τ = 0.35 (black, solid line), which was obtained from de Vries et al., 1989), this parametrisation does not capture the large spread in the observations (blue, green and red lines). The uncertainty this introduces into the model is not critical, since there are no strong feedbacks involved (unlike, for example, uncertainty in specific leaf area (SLA) just after emergence), but it will affect the outputted yield.

Senescence
The observations of green leaf biomass and above-ground biomass were used to tune the senescence parameters µ, ν and DVI sen . The above-ground biomass measurements were combined with the partition fractions from Sect. 2.3.2, the carbon to biomass ratios from Sect. 2.3.7 and the senescence parametrisation from Eq. (A4) to get a time series for green leaf biomass (Fig. 4, centre and right plots, black lines), normalised to the maximum value in each year. This could then be compared to the normalised observed time series for green leaf biomass (Fig. 4, left, coloured lines). It is clear that, if the parametrisation from Osborne et al. (2015) is used (Fig. 4, centre plot, solid black lines), senescence starts late and then progresses too abruptly as compared to the observations. However, with the new parametrisation (with the new free parameters µ, ν and DVI sen ), it is possible to get a much better fit to the observations (Fig. 4, right plot, dashed black lines). Note that this tuning partially compensates for the bias in the proportion of carbon going to the leaves between DVI 0.8 and 1.0 in Fig. 1 (top). If this bias was not present, senescence could start more gradually, which would enable a better fit to leaf carbon at around DVI = 1.75. Also, the tuned lines underestimate the leaf biomass at around DVI = 1.75, which will help to compensate for the model being unable www.geosci-model-dev.net/10/1291/2017/ Geosci. Model Dev., 10, 1291-1320, 2017 to capture the drop in photosynthetic capacity in the green maize leaves towards the end of the season.

Crop height
Stem biomass measurements up until the maximum in each year and the corresponding crop height measurements from the Mead FLUXNET sites were used to fit the allometric constants κ and λ, through rearranging Eq. (A5) to h = κ M λ stem where κ = κ(1 − τ ) λ . For consistency, it is important that the τ used in this expression is the same value as the τ used in Eq. (A2). Figure 5 shows the observations (points), along with the fit using parameters from Osborne et al. (2015) (solid black line, λ = 0.4, κ = 3.06) and a tuned fit (dashed black line, λ = 0.38, κ = 3.43).

Specific leaf area
The allometric constants γ and δ relating specific leaf area to DVI (Eq. A7) are tuned using Fig. 6, which plots SLA observations against DVI (points), the tuned fit (dashed line) and the parameters used in Osborne et al. (2015) (solid line). The crop in the model is very sensitive to SLA for low values of DVI because of the feedback between leaf area index and leaf carbon. The model lines in Fig. 6 have the steepest gradient for low values of DVI, where there is also a greater spread of observations.

Carbon to biomass ratio in stem and leaves
The observations of carbon fraction of the green leaf biomass (canopy mean) against day after sowing is shown in Fig. 7.

Initial amount of carbon in crops
Assuming that, near emergence, approximately half of the plant carbon is above ground (Fig. 2), values for the parameters governing initialisation of C init = 8.0 × 10 −4 and DVI init = 0.1 can be derived from the above-ground biomass measurements plotted in Fig. 8 and the carbon to biomass ratios. Since there are no measurements below DVI = 0.1, and the model is very sensitive to these parameters, we do not attempt to set a DVI init below 0.1 and extrapolate. Note also that the initial value of carbon is very sensitive to the ther- Table 3. Values of the crop-specific JULES parameters used to represent maize. Units are given in brackets; (-) denotes dimensionless. These parameters are all specified in the JULES_CROPPARM namelist. Values of the crop-specific JULES parameters used to represent maize. Units are given in brackets; (-) denotes dimensionless. These parameters are all specified in the JULES_CROPPARM namelist.

JULES
Osborne et mal time for emergence. Figure 8 also shows that the value C init = 1.0 × 10 −2 , which was used in Osborne et al. (2015) to initialise the crop at DVI init = 0.0, is too high to be consistent with the above-ground biomass observations.

Yield fraction
As discussed above, the C harv pool in JULES contains both the reproductive parts of the maize crop (kernel, cob, husk, ear shank and silk) and the yellow leaf carbon and the proportion of this carbon pool that contributed to yield carbon f yield is set by the user. The value of f yield can be derived using the latest observations in each season of the biomass of the reproductive part of the crop, the proportion of this reproductive biomass which is composed of kernels, and the yellow leaf biomass. The yield fraction is then calculated as the kernel fraction of the sum of the reproductive part of the crop and the yellow leaves, leading to an approximate value of f yield = 0.74 ( Fig. 9). This assumes that there is no significant change in f yield between the last measurement of the season and the harvest and also that the carbon fraction of the biomass in the both the reproductive parts and the yellow leaves is the same. Typically, an accurate value of f yield is not important in impact studies, since this constant can be incorporated into a yield gap parameter.

Parameters required by natural PFT tiles only
To obtain the allometric parameters required to relate the plant carbon pools to plant height and LAI when the crop model is switched off, LAI bal was assumed to be approximately equal to LAI up to the maximum LAI at the site for each year. As discussed in Sect. A2, a ws is assumed to be equivalent to 1-τ , i.e. a ws = 0.88. The stem biomass observations can be used to obtain values for a wl , b wl and η sl , for a set ratio of carbon to biomass in the stem (see Sect. 2.3.7). First, a value for η sl of 0.017 kg C m −1 (m 2 leaf) −1 was obtained by plotting the stem biomass observations against LAI multiplied by crop height for points up until the maximum LAI for each site in a particular year (Fig. 10, left). Second, a wl and b wl were simultaneously fitted to (a) the stem biomass observations against LAI for points up until the maximum Geosci. Model Dev., 10, 1291-1320, 2017 www.geosci-model-dev.net/10/1291/2017/  LAI for each site in a particular year (Fig. 10, right), (b) crop height against stem biomass observations, up until the maximum stem biomass measurement for each site in a particular year (Fig. 11) and (c) LAI bal against LAI up until the maximum LAI for each site in a particular year (Fig. 12). This gave a wl = 9.5 × 10 −3 kg C m −2 and b wl = 1.767.
As we saw in Fig. 6, Eq. (A12) is not a good approximation for maize, particularly when DVI is less than 0.5. For the purpose of these runs, an approximate value at DVI = 1 was used.

Parameters required by both crop tiles and natural
PFT tiles 2.5.1 Canopy radiation scheme The JULES default C4 grass settings for the PAR leaf scattering coefficient ω PAR = 0.17 and the PAR leaf reflection coefficient α refl,PAR = 0.1 were used (these are very similar to the values quoted in Sellers (1985)  that were used in Osborne et al. (2015) to model maize. The soil albedo was set to 0.133, which was the value from the nearest grid box in the ancillary used in the HadGEM2-ES model (Collins et al., 2011;Jones et al., 2011), which was used in the Osborne et al. (2015) global runs.
The canopy clumping factor was tuned by comparing the fraction of incident PAR absorbed by the canopy (frac- tion of absorbed PAR (FAPAR)), using absorbed and incident PAR observations and interpolated LAI observations, to the model FAPAR, using observed diffuse radiation fraction and interpolated LAI observations up until flowering. The python package pySellersTwoStream (see code availability section) was used to calculate the model FAPAR since it is able to reproduce the results of the JULES radiation scheme exactly but can be called directly from our (python) Geosci. Model Dev., 10, 1291-1320, 2017 www.geosci-model-dev.net/10/1291/2017/  Figure 9. Yield fraction against the sum of the biomass in the reproductive parts of the maize crop (kernel, cob, husk, ear shank and silk) and the yellow leaf biomass, using the last measurement of the season. Dots, vertical crosses (+) and diagonal crosses (x) are US-Ne1, US-Ne2 and US-Ne3 observations respectively. Solid black line shows the value used implicitly in Osborne et al. (2015) and dashed black line shows the new, tuned value.
analysis scripts, without the need for extra JULES runs for each combination of parameters tested.
Absorbed PAR through the canopy in the model closely follows a exponential decay function. Calculating FAPAR involves integrating this exponential decay over the canopy; Fig. 13 (centre row) shows the resulting FAPAR distribution against total LAI for a uniform canopy (canopy clumping factor a = 1). For mostly direct radiation (diffuse radiation fraction 0.2-0.3), the rate of decay with layer LAI in the model shows a clear dependence on the zenith angle (Fig. 13, centre right), whereas for mostly diffuse radiation (diffuse radiation fraction 0.8-0.9), this zenith angle dependence is greatly reduced (Fig. 13, centre left). While the observations (Fig. 13, top row) also show a strong zenith angle dependence as the fraction of diffuse radiation decreases, the observations are, in general, consistent with a much lower effective decay constant (in particular, the model FAPAR values are higher than the observations at intermediate LAI values ∼ 2). The observed FAPAR values also have a much larger scatter than seen in the model FAPAR.
Decreasing the canopy clumping factor is equivalent to decreasing the effective decay constant in the model. Figure 14 shows the value of the clumping factor that would be needed to reproduce each FAPAR observation, given the observed LAI and diffuse radiation fraction. While there is a large spread in clumping values derived in this way, these re- sults appear to indicate that a clumping factor between 0.5 and 0.8 would be consistent with the majority of the observations. In this study, we therefore set a = 0.65. Figure 13 (bottom row) shows that using this clumping factor value to calculate model FAPAR gives a better fit to the observations, particularly for the intermediate LAI values. Erectile, vertical and horizontal leaf angle distributions (for a uniform canopy) were also investigated, but the spherical distribution gave the best fit to the FAPAR observations. The FAPAR observations can not be used to tune the model once green leaf area index has started to drop significantly, as the observations include PAR absorbed by any part of the plant, whereas the JULES canopy scheme models the PAR absorbed by photosynthesising leaves only. Whether the model canopy scheme needs to be extended to include the shading of green leaves by yellow leaves and other non-root biomass depends on the distribution of the remaining green leaves through the canopy (essentially, the model is roughly assuming that all the green LAI is at the top of the plant and so does not get shaded by other plant material). Differ-Geosci. Model Dev., 10, 1291-1320, 2017 www.geosci-model-dev.net/10/1291/2017/  Sellers (1985) modelled maize assuming that green and dead leaves are evenly distributed throughout the canopy, whereas de Vries et al. (1989) showed that "maximum leaf photosynthesis in a senescencing crop declines with time. The oldest leaves in the base of the canopy are affected first".

Photosynthesis light response curve
In the literature, the photosynthetic capacity of maize leaves (per leaf area) declines with age and the older leaves are lower in the canopy Dwyer and Stewart (1986); Stirling et al. (1994). As discussed in Sect. A4, change in photosynthetic capacity through the canopy can be modelled in JULES by a non-zero k nl , which we assume is due to change in nitrogen per unit leaf area through the canopy.
The nitrogen per unit leaf area as a function of layer LAI at anthesis (60 days after sowing) in Massignam et al. (2001) for the highest nitrogen availability level (150 kg N ha −1 ; residual soil nitrate 31 kg ha −1 ) was consistent with a k nl of approximately 0.07. Since this is low, in this study, the variation of nitrogen per unit leaf area through the canopy is neglected; i.e. k nl = 0.0. The inclusion of a non-zero k nl would have the effect of increasing GPP, as the plant would be able to make more efficient use of the incoming radiation.
In this study, trait-based physiology was switched off (i.e. l_trait_phys=F). However, the same results could be obtained by switching trait-based physiology on and choos- ing values for the new parameters that are equivalent to the ones used here. Figure 15 shows the observations of the nitrogen mass per unit carbon mass (left) and per unit leaf area (right) averaged over the canopy. In both plots, nitrogen rapidly decreases with time at the beginning and end of the season, which cannot be captured by JULES. The inclusion of a non-zero k nl would also not solve this problem, as this would simply increase the nitrogen per leaf area mid-season, as can be seen in Fig. 16 for k nl = 0.2.
In this study, the temperature dependence of V cmax is fixed by fitting Eq. (A14) to the expression given in de Vries et al. (1989) (Fig. 17). The default JULES C4 grass parametrisation of V cmax is more sharply peaked, has its maximum at a higher temperature and is more asymmetrical. Also plotted is the expression for the temperature dependence for maize V cmax from Massad et al. (2007). Puntel (2012) modelled V cmax for maize at the Mead site and fit the results with MaizeGro, using the default temperature dependence, which gives a peak at approximately 33 • C. Puntel (2012) verified this relation by successfully fitting the model to results from modern maize cultivars from Kim et al. (2007), Crafts-Brandner and Salvucci (2002) and Naidu et al. (2003), which all show the peak in V cmax at approximately the same temperature. Puntel (2012) related the normalisation of V cmax to the leaf nitrogen per biomass; for example, at 30 g N kg −1 at the V14 growth stage, maximum assimilation at 25 • C was 37 µmol m −2 s −1 . The temperature dependence of maize at high temperatures was looked at in more detail in Crafts-Brandner and Salvucci (2002), which included an investigation into the dependence on the rate of temperature change.
The experiment with the more gradual temperature change in Crafts-Brandner and Salvucci (2002) corresponds well to the high temperature dependence of the de Vries et al. (1989) expression.
The canopy average V cmax,norm was tuned using the value of V cmax at 25 • C at 340 vppm CO 2 at a specific leaf weight of 450 kg h −1 , which is the canopy average at DVI = 1 (for maize cv Pioneer) from de Vries et al. (1989). n l0 was set to the approximate value of the observations in Fig. 15 (left) at DVI = 1, which then constrains n e (since n e = V cmax,norm /n l0 when k nl = 0). The quantum efficiency α was set to the value from de Vries et al. (1989) of 0.055 µmol C m −2 s −1 (µmol photons m −2 s −1 ) −1 for maize, which was quoted for temperatures lower than 45 • C (above this temperature, it drops sharply -an effect which is not reproduced in JULES). This is consistent with values in the literature (e.g. Massad et al., 2007, and references therein) and consistent with the fitted values of α from Puntel (2012). The value of α for maize is not dependent on leaf age or position (Dwyer and Stewart, 1986). This method of tuning the JULES parameters has assumed that the two limiting rates are predominantly W c and W light , not W e .
Note, however, that the photosynthesis light response curve in de Vries et al. (1989) has an exponential dependence on the absorbed radiation, which causes the shape to vary slightly from the non-rectangular hyperbolae used in JULES (with hard-wired values of curvature from Collatz et al., 1992), leading to lower values of photosynthesis below approximately 1500 µmol photons m −2 s −1 .
The parameters involved in calculating the leaf internal carbon dioxide partial pressure, q crit and f 0 (in Eq. A19), were not expected to strongly limit the results since this current study focusses on carbon fluxes rather than water fluxes, the runs are irrigated and the rate W e is not expected to be limiting. q crit was left at its default C4 grass value (as in Osborne et al., 2015) and f 0 was set to 0.4 (consistent with the range of maize measurements quoted in de Vries et al., 1989).

Respiration
Values for µ rl and µ sl (from Eq. A22) were obtained for maize from de Vries et al. (1989) of µ rl = 0.39 and µ sl = 0.43 (note that this assumes one constant value for the nitrogen per carbon in leaves over the crop season and τ = 0.12). Fixing the value for the dark respiration coefficient f dr (used in Eq. A20) is complicated by the inclusion in the code of inhibition of leaf respiration in the light. Also, Atkin et al. (1997) demonstrated that the dark respiration in darkness decreases as the time the leaf has been in darkness increases. This complicates the use of the light response curves for fitting this parameter, since this means that parameters measured during the day will not necessarily correspond to those needed in JULES for modelling the average dark respiration over a 24 h period. Using de Vries et al. (1989) values for the maximum rate of leaf photosynthesis at 450 kg biomass per hectare and maintenance respiration at 25 • C for maize gives f 24 h dr = 0.0081 over the course of 24 h. Even with a correction for inhibition of dark respiration in the light, this is inconsistent with the spread of fitted values of dark respiration to maximum assimilation to light response curves measured at the site between 10:00 and 14:00 local time, presented in Puntel (2012) (leaf is exposed to ambient light pre-measurements), which are much higher, unless the dark respiration derived from the light curves is assumed to have a contribution from what JULES considers the "growth respiration". In general, the dark respiration coefficient estimated from light response curves for maize appears to be higher than the value derived from the maintenance respiration measurement in de Vries et al. (1989) (e.g. Collatz et al., 1992;Dohleman and Long, 2009), which is consistent with there being a component from growth respiration. In our JULES runs, we will use f dr derived from the maintenance respiration observation in de Vries et al. (1989), corrected assuming that in the day of measurement 50 % of leaves experienced inhibition of the dark respiration by light; i.e. f dr is set to 0.0081/0.85 = 0.0095 (this assumption was later tested, and found to be accurate to within 2 %).
de Vries et al. (1989) gave a growth respiration coefficient of 0.22, 0.18, 0.19 and 0.18 for maize leaves, stem, roots and cob/grain respectively. These values can not be used directly in JULES since, as described earlier, the growth respiration coefficient in JULES is a constant for each carbon pool. Here, we set r g to 0.25 for every PFT, as in the JULES Global Land (GL4.0) configuration (Walters et al., 2014) (note, however, that this approximation of a constant r g for each plant carbon pool would break down for other crops e.g. soybean).
It is also worth noting that Puntel (2012) found that the maximum assimilation rate had a much stronger relationship with leaf nitrogen than the leaf dark respiration rate. In addition, Stirling et al. (1994) shows a strong dependence in dark respiration in maize over time (using fits to light response curves), which can not be captured in JULES: at degree day 220 (roughly where the leaf area reaches a maximum), it is approximately twice as high at degree day 50. As we have discussed, maintenance respiration and V cmax covary in JULES, but the growth respiration is linked to net primary productivity, which increases in the crop up until approximately anthesis. Therefore, the total leaf respiration in   Osborne et al. (2015) and this study respectively. the model will vary in time, and will have a different dependence on time to V cmax . However, the issues we have already identified with the modelling of the evolution of V cmax over time will impact the accuracy of the modelling of the maintenance component of the leaf respiration over time. Leaf dark respiration rates also differ between different maize hybrids (Earl and Tollenaar, 1998). There is therefore a large uncertainty in the parameter f dr and the overall determination of growth respiration.

Results and Discussion
In this section we present the results from the JULES runs and compare with observations from the Mead sites. The runs with the crop model switched off and prescribed LAI and height are useful for evaluating the parameter choices for photosynthesis and respiration, without the additional complication of the feedback between LAI and NPP, as will be discussed first. The results from the full crop-model configuration will then be evaluated.

Gross primary productivity
Plots of modelled GPP (blue) against observed GPP (green) are shown in Fig. 18 for years in which irrigated maize was grown at the Mead FLUXNET sites US-Ne1 and US-Ne2. While the overall shape of the plots is good, it is clear that GPP in the model is significantly overestimated after the midseason peak in observed GPP (corresponding to where LAI declines as the crop leaves senesce). As discussed earlier, the model V cmax at a certain temperature stays constant, whereas in reality it would decline over the crop season. Implementing this decline into JULES would result in a much closer fit between the model GPP and observed GPP. To a lesser extent, there also appears to be an overestimation of GPP in the model before senescence. This was investigated in more detail by comparing plots of FLUXNET GPP against observed Absorbed Photosynthetically Active Radiation (APAR) with plots of model GPP against APAR, for hourly measurements before the crop reaches DVI = 1, for LAI bins of size 1. Figure 19 shows the LAI bin 3.5 to 4.5. There is a clustering of points due to the hourly resolution of the data, which is most clearly seen in the model output. Hours with high diffuse radiation fractions (red) are similar in both the FLUXNET data and the model output, although the scatter in the FLUXNET data is higher, as expected from the plots of observed FAPAR (Fig. 13). For lower diffuse radiation fractions in the model, GPP decreases due to a combination of the effect of sunflecks and an increase in the effective decay constant of absorbed PAR through the canopy at the beginning and end of the day. Even when the scatter in the FAPAR observations is taken into account, the decrease in GPP for lower diffuse radiation fractions does not appear to be as large in the model as in the GPP observations, and this is the source of the overestimation of GPP we saw in the model output in Fig. 18 before the onset of senescence.
This effect was investigated further by considering the dependence on air temperature and vapour pressure deficit in the FLUXNET GPP data. As expected, the lower temperature points (Fig. 20, top left) and lower vapour pressure deficit (VPD) points (Fig. 20, top right) are clustered at low values of APAR. However, there does not seem to be a dependence on temperature or VPD at a constant APAR across the range of GPP observations. Soil moisture stress is a factor that we have neglected in our runs (since we have assumed perfect irrigation), which could, if implemented, reduce GPP when the soil moisture is low. However, as Fig. 20 shows for soil moisture content at a depth of 10 cm (bottom left) and 25 cm (bottom right), at higher APAR values, points below a threshold of 30 % appear to be distributed evenly across the range of GPP observations for a constant APAR.
Including a decrease in leaf nitrogen concentration through the canopy (while keeping the total amount of nitrogen constant) would have the effect of making the light use of the plant more efficient, which would increase model GPP still further. Decreasing V cmax,norm would have the effect of decreasing model GPP at higher APAR values, but this would not solve the issue at mid-range APAR points ∼ 800 µmol photons (m 2 ground) −1 s −1 and would also worsen the fit of the points with high diffuse radiation fractions. It is therefore difficult to see a clear way in which the model parameter settings or processes should be improved. It would be possible to improve the validation against observations by decreasing α or changing the curvature parameter in the non-rectangular hyperbola implemented for light response within JULES (currently hard wired) but it is difficult to justify this theoretically.

Respiration
The results from the model runs without the crop model can also be used to test the parametrisation of respiration.
Using a number of assumptions, the measurements from Mead can be used to get an approximate value for leaf main-tenance respiration. First, approximate values for NPP were obtained by linearly interpolating the Mead above-ground biomass measurements to get a daily time series, and then differentiating. The fraction of NPP directed to the roots at each DVI was calculated from the expression for maize in de Vries et al. (1989) (plotted in Fig. 2) and then used to obtain the total NPP. Combining these NPP values with the GPP observations and assuming a value for the growth respiration coefficient of r g = 0.25 and summing over the crop season leads to an estimation of the plant maintenance respiration R pm . It is necessary to sum over the whole season, since the NPP and GPP calculated in this way appear to be slightly out of step with each other, and this effect dominates the daily time series of derived maintenance respiration.
The interpolated carbon pool observations were used to calculate the factor 1 + µ rl that converts between the leaf maintenance respiration and the total plant maintenance respiration. Note that the stem carbon observations had to be corrected using τ to get C stem . This factor was used to convert the leaf maintenance respiration outputted by the model to the total plant maintenance respiration. Figure 21 shows the R pm derived from observed GPP against R pm /f dr derived from the outputted model leaf maintenance respiration. The x axis therefore is independent of f dr , which can be obtained from the gradient. Data from 2010 is not included (since the crop was damaged by hail). Both the default JULES C4 grass f dr (solid line) and the f dr used in our maize configuration (dashed line) are shown. It can clearly be seen that the new maize f dr is a better fit than the default C4 grass value. While there are many model and parameter assumptions (r g , µ rl , µ sl , β = 1, C root , τ ) that have gone into this plot, this is still an important consistency check of our parameters.  Tables 1, 2, 3, and 4 Figure 22 compares the model GPP and the observations, and shows very close agreement. This is influenced by a cancellation of two effects: as identified in the previous section, the GPP per APAR in the model is biased high, whereas the outputted LAI is biased low, shown in Fig. 23. In part the reduction in modelled LAI compared to observations was deliberately introduced when tuning the senescence parameters so that a quicker decrease in LAI compensates partially for the model not including a decrease in leaf photosynthetic capacity. However, it is also clear that the interannual variability of LAI is not reproduced by the model; in particular, in some years (2006,2010,2011 for US-Ne1 and 2011 for US-Ne2), the LAI is too small in the crop season up to anthesis. This is due to the high sensitivity of the plant in its early life to parameter settings, due to the feedback between NPP and LAI. In these site and year combinations (2006,2010,2011 for US-Ne1 and 2011 for US-Ne2), temperatures between DVI 0.1 and DVI 0.2 are higher on average, and so DVI is increasing more rapidly, which gives the plant less time to accumulate NPP, leading to a reduced rate of increase of LAI with respect to DVI in the model runs at this growth stage. On the other hand, the SLA observations for these years in the early crop season are particularly high compared to the rest of the distribution, which means that the observations do not show this reduced rate of increase of LAI at this growth stage. Fitting γ and δ to the SLA observations in just these site and year combinations (2006,2010,2011 for US-Ne1 and 2011 for US-Ne2) gives 18.0 and −0.45 respectively. Using these parameters in JULES runs with the crop model gives much better agreement with LAI observations (Fig. 24). This is also consistent with the result from US-Ne2 in 2010: since the crop emerges 9 days after the crop in US-Ne1, the period of relatively high temperatures mostly falls before the crop is initialised. It is possible that parametrising SLA with day after emergence rather than with DVI might improve the fit between model and observed LAI by reducing the sensitivity of the SLA parametrisation to temperature.

Results from JULES runs with the crop model
The canopy height is well represented in the runs (Fig. 25). The above-ground carbon in the model also fits the observations well (Fig. 26). The harvest carbon pool (which includes the reproductive parts of the plant and the yellow leaves) is overestimated in the model, which is consistent with the overestimation of GPP during the senescence period.

Conclusions
The JULES-crop parametrisation of crops within JULES was introduced to improve the carbon and energy fluxes in the model over croplands and to investigate the effect of weather and climate on food and water resources, at global, regional and local scales. In this evaluation paper, we have looked in detail at how the input parameters in this pre-existing model can be tuned for one crop (maize) at one location (Mead, US), where there are a wide variety of observations to probe how the model components perform, both separately and in combination.
In previous analyses with JULES-crop, it has been assumed that model photosynthesis and respiration parameters can be set to the default C3 grass values for C3 crops and the default C4 grass values for C4 crops. We have used literature results and the observations available at this site to improve the maize parameters required in both the crop-model part of H JULES (such as partition fractions and allometric constants) and the generic vegetation code.
With the new parameters, there is good agreement between modelled GPP and observed GPP up until anthesis if the feedback between NPP and LAI is removed by switching the crop model off and prescribing LAI (and canopy height) when the skies are mostly overcast. The model tends to overestimate GPP for clearer skies. After anthesis, there is a much greater overestimation of GPP, due to the model being unable to capture the decrease in photosynthetic capability at the leaf level over time in the crop. The respiration parameters were more difficult to test in isolation, but integrating model respi-ration over the entire crop season produced results that were consistent with the GPP and carbon pool observations.
Running the full crop model, including all the new parameters, produced GPP time series that were very close to the observations. This was helped partially by a cancellation of two biases -the model GPP for a certain LAI was biased high, as we have just discussed, and the LAI in the model was biased low compared to the observations. There were a few anomalous years in which the peak LAI in the model was approximately two-thirds that of the peak LAI in the observations, which may imply oversensitivity to initial conditions. The amount of above-ground carbon was reproduced well, although the amount of carbon in the harvest pool was overestimated in most cases.
There should be three main priorities for extending this work to improve the representation of maize at these sites. First, work should be done to tune the parametrisation of soil moisture stress of maize so that the water balance of the irrigated sites could be accurately modelled and runs for the non-irrigated site could also be included. Second, a parametrisation of the maximum rate of carboxylation of Rubisco V cmax should be added that allows it to vary over the course of the crop season. Third, these runs have been tightly constrained by using observed sowing, emergence, flowering and harvest dates to generate the thermal times needed as input to JULES. For most regions, and for any climate projections, this sort of data will not be available. Therefore, it would be a useful test of the model to investigate the performance at the Mead sites if the model is given generic values for the thermal time parameters.
While this study has focussed on modelling one crop variety at one site, it also provides a demonstration of how knowledge of the structure of the model can be used to tease apart different components of the model so that they can be tuned or evaluated against observations. This ranged from the tuning of parameters in simple allometric relations such as that relating stem carbon to canopy height, to tuning the canopy parameters using the external representation of the  Tables 1, 2, 3, and 4. canopy scheme in pySellersTwoStream, up to running JULES with the crop model switched off and prescribed LAI and canopy height, in order to tune GPP without the complication of the feedback between GPP and LAI. It therefore provides a case study, which can be used when setting up and evaluating the model for other crop varieties and sites.
Code availability. This study uses JULES revision 5061, which is between the 4.6 and 4.7 releases. The code can be downloaded from the JULES FCM repository at https://code.metoffice.gov.uk/ trac/jules/ (JULES collaboration, 2017) (registration required). The pySellersTwoStream package is available at https://github.com/ tquaife/pySellersTwoStream (Quaife, 2016). The version used in this study was downloaded on 15 September 2016.
Data availability. Unless otherwise noted, all site observations discussed in this paper were obtained from the Site Information pages of the AmeriFlux website hosted by Oak Ridge National Laboratory (http://public.ornl.gov/ameriflux/, AmeriFlux collaboration, 2016) or by personal communication with the Mead sites Research Technologist.
Note: these data are currently being transitioned to a new location: http://fluxnet.fluxdata.org/.

Appendix A: Model description
In this section, we will summarise the relevant features of JULES and the JULES-crop parametrisation within it, paying particular attention to new model features available since the Osborne et al. (2015) study (i.e. post version 4.0). These new options are indicated in Tables 1, 2, 3, and 4.

A1 Crop model
In JULES-crop, the development status of each crop within a grid box is parametrised by a crop development index (DVI). DVI is −2 before sowing, −1 at sowing, 0 at emergence and 1 at flowering. Under favourable conditions, harvest occurs at a DVI of 2. The DVI has three main functions within the JULES-crop model: it determines the harvest date, the partitioning of NPP between the crop carbon pools and the dependence of the specific leaf area on leaf carbon.
The increase in DVI over the course of the crop's lifetime is determined by crop-specific thermal time parameters, set by the user. If the dependence on photoperiod length is neglected (as in Osborne et al., 2015), thermal time becomes an accumulation of effective temperature between one development stage and the next, where effective temperature is defined by i.e. a triangular function, peaking at an optimal temperature T o , which is zero below a base temperature T b and above a maximum temperature T m . T o , T b and T m are parameters specified by the user for each crop. T o , T b and T m are given in Kelvin and thermal time in units of degree days. Crop growth is modelled by accumulating net primary productivity over the course of a day (NPP acc ) and splitting this carbon between the crop root, stem, leaf, harvest and reserve carbon pools for that tile (C root , C leaf , C stem , C harv and C resv respectively) according to where τ is the fraction of stem carbon that is partitioned into the stem reserve pool (containing the remobilisable carbohydrates) and p i (for i = root, stem, leaf, harv) are the partition coefficients defined by where j = root, stem, leaf, harv. α i and β i are numerical constants that are tuned to observational data. α harv and β harv are both set to zero. All other α i and β i are set by the user for each crop. Note that j p j = 1.
The crop carbon pools are initialised at DVI init , which is at or just after emergence. At initialisation, the crops are given a certain amount of carbon C init , which is distributed between the carbon pools according to the values of p i at DVI = DVI init .
Once p stem drops below 0.01, carbon from the stem reserve pool is mobilised to the harvest pool, by reducing C resv by 10 % each day and adding this carbon to the harvest pool (as proposed in de Vries et al., 1989). Similarly, once the DVI is above a threshold value DVI sen , carbon from the leaf pool is mobilised to the harvest pool to simulate leaf senescence, by reducing C leaf by a fraction, each day when DVI > DVI sen . ν and µ are numerical constants that are tuned to observational data. After DVI init and if the sowing date is prescribed, the model harvests the crop and resets the crop tile if any of the following conditions are satisfied: 1. DVI reaches 2 (i.e. the desired harvest condition); 2. LAI > 15, since once the model reaches such large LAI it is clearly unrealistic; 3. the temperature of the second soil layer from the top falls below a user-defined temperature T mort at any time after DVI = 1; 4. DVI > 1.0, the carbon in the roots, leaves, stem and stem reserve pool of the crop falls below C init and the amount of carbon in the harvest pool is greater than zero; 5. the crop age reaches 1 year, so that a new crop can be sown each year.
The crop height h is calculated from the C stem pool using where κ and λ are allometric constants and f C,stem is the fraction of carbon in the dried stem (excluding the stem reserves), all given as input by the user.
The green (i.e. photosynthesising) leaf area index (LAI) is calculated from the leaf carbon and the specific leaf area (SLA) by where f C,leaf is the carbon fraction of the dry leaves. The SLA depends on the DVI via where γ and δ are allometric constants which are set by the user. JULES-crop outputs water-limited potential yield if irrigation is switched off and potential yield if irrigation is on, expressed in kg C m −2 . This yield is calculated by multiplying the value of C harv on the day of harvest by a parameter f yield supplied by the user, which represents the fraction of C harv that is economically valuable, i.e. the maize kernel in our runs.

A2 Relationship between LAI, canopy height and plant carbon for natural vegetation
When the crop model is switched off, different allometric functions are used to approximate the carbon in the leaf, stem and root pools based on the prognostics LAI and canopy height h. These allometric functions make use of a "balanced" leaf area index (LAI bal ), which is calculated from canopy height using LAI bal = a ws η sl a wl h where a ws , a wl , η sl and b wl are all allometric constants, defined in relation to the respiring stem carbon S and the total stem carbon W: We assume here that S is equivalent to C stem and W is equivalent to C stem + C resv in the crop model. Therefore, a ws is equivalent to 1 − τ in the crop model and these equations can be compared directly to Eq. (A5) until the start of the remobilisation of the crop stem reserve pool. The size of the leaf carbon pool C leaf is calculated by multiplying the LAI by the canopy-averaged specific leaf density σ l (in kg C (m 2 leaf) −1 ), which is assumed to be constant, i.e.
The root carbon C root is approximated by

A3 Canopy
JULES has a number of options for calculating the photosynthetically active radiation (PAR) available to leaves at different depths in the plant canopy. In this discussion, we will focus on the canopy radiation scheme used in Osborne et al. (2015) (can_rad_mod 5) and the canopy radiation scheme currently recommended for layered canopies in JULES (can_rad_mod 6), which both treat the direct and diffuse components of the incident radiation separately (as in Sellers, 1985) and include sunflecks. We also assume a zenith angle dependence (l_cosz=T). JULES assumes that the incident PAR is half of the incident shortwave radiation. The amount of incident PAR composed of diffuse radiation is given as part of the driving data. The canopy is split into 10 equal layers of green leaf area index (LAI). The equations for absorption and scattering at each layer for the incident diffuse beam and the incident direct beam are solved separately, taking into account the distribution of leaf angles and the zenith angle. The sunlit fraction of the leaf is also calculated, and absorbs light from the direct component of the direct beam radiation ("sunflecks"), in addition to the diffuse light from the direct beam and light from the diffuse beam. The shaded fraction of the leaf absorbs light scattered from the direct beam and light from diffuse beam only (i.e. no direct sunlight). JULES has two leaf angle distributions currently implemented -spherical and horizontal. As of JULES version 4.6, JULES also includes a canopy clumping factor a, which scales LAI within the canopy radiation scheme and represents variation within and across canopy structures.

A4 Modelling C4 photosynthesis
In JULES, potential leaf-level photosynthesis (unstressed by water availability and ozone effects) is calculated as the smoothed minimum of three rates, following Collatz et al. (1991Collatz et al. ( , 1992: (a) the Rubisco-limited rate W c , which depends on the maximum rate of carboxylation of Rubisco, (b) the light-limited rate W light and (c) the rate associated with the transport of photosynthetic products for C3 plants or PEP (phosphoenolpyruvate) carboxylase limitation for C4 plants W e . For C4 plants, W c is set to the maximum rate of carboxylation of Rubisco, V cmax . V cmax is calculated using and V cmax (T c = 25 • C) are within 5 % of each other. T upp and T low are used to give the leaf an optimum temperature range, which is superimposed on the Q 10 dependence in f T . If trait-based physiology is switched off in JULES (l_trait_phys=F) V cmax,norm = n e n l , where n l is the mass of nitrogen per mass of carbon in the leaf (with units kg N (kg C) −1 ), which varies through the canopy, and n e is a normalisation constant, fitted to data. The input parameters specified by the user are n 0 l (n l at the top of the canopy) and n e .
In the JULES canopy radiation scheme can_rad_mod 5, V cmax,norm is assumed to vary through the canopy according to exp(−k n LAI layer /LAI). In can_rad_mod 6, V cmax,norm varies through the canopy according to exp(−k nl LAI layer ). k n and k nl are PFT-dependent parameters set by the user.
The light-limited rate of leaf photosynthesis for C4 plants is calculated in JULES using where α is the quantum efficiency in mol CO 2 (mol PAR photons) −1 and I APAR is the absorbed photosynthetically active radiation (APAR) in mol PAR photons m −2 s −1 . As discussed, can_rad_mod 5 and can_rad_mod 6 include the effect of sunflecks by spitting the leaf into a sunlight and a shaded part, which have different values of I APAR and therefore different W light . The rate associated with PEP carboxylase limitation W e in JULES is where P * is the surface air pressure and c i is the leaf internal carbon dioxide partial pressure, which is calculated for C4 plants using where is the photorespiration point (zero for C4 plants) and c a is canopy CO 2 pressure. q is the canopy level specific humidity deficit, q crit is the critical specific humidity deficit and f 0 is the ratio of c i to c a at which the canopy level specific humidity deficit is zero. c a is calculated from R CO 2 P * / , where R CO 2 is the atmospheric CO 2 mass mixing ratio and = 1.5194 is the ratio of molecular weights of CO 2 and dry air. As an example, for zero specific humidity deficit, an atmospheric CO 2 mass mixing ratio of 5.6×10 −4 (2003 global average; Dlugokencky and Tans, 2016), f 0 = 0.8 (JULES C4 grass default), the value of W e is 5.9V cmax .
The rate of gross leaf photosynthesis W is the smoothed minimum of W c , W light and W e (calculated using nonrectangular hyperbolic functions with the curvature parameters hard wired). The net potential (i.e. unstressed) leaf photosynthetic carbon uptake A p is the gross leaf photosynthesis minus the dark leaf respiration R d . The potential leaf photosynthesis is converted to a net photosynthesis by multiplying by a soil water stress parameter β. Stomata at points with negative or zero net photosynthesis or where the leaf resistance exceeds its maximum value are closed (i.e. leaf gross photosynthesis is zero). Leaf resistance is calculated from the net (i.e. water-limited) rate of photosynthesis, (c a − c i ), the leaf temperature and the ratio of leaf resistance for CO 2 to leaf resistance for H 2 O (= 1.6).

A5 Respiration
In JULES, the (non-water-limited) leaf dark respiration R d (in mol CO 2 (m 2 leaf) −1 s −1 ) is calculated by R d = 0.7f dr V cmax for I APAR LAI > 10 µmol CO 2 (m 2 ground) −1 s −1 f dr V cmax otherwise (A20) to allow for the inhibition of dark respiration during daylight. R d is summed over the canopy levels for sunlit and shaded leaves to get R dc , the canopy dark respiration in (in mol CO 2 (m 2 ground) −1 s −1 ). The plant maintenance respiration in kg C (m 2 ground) −1 s −1 is calculated (for the setting l_scale_resp_pm=T) using = 0.012R dc β 1 + µ rl where N root , N stem and N leaf are the nitrogen in the roots, stems and leaves respectively. µ rl is the mass ratio of nitrogen to carbon in the roots divided by the ratio of nitrogen to carbon in the leaves. µ sl is the mass ratio of nitrogen to carbon in the stem (not including stem reserves) divided by the ratio of nitrogen to carbon in the leaves. The factor 0.012 relates mol CO 2 to kg C. If the option l_scale_resp_pm=F is set, the root and stem terms do not depend on β.
In JULES, plant growth respiration R pg is a fixed fraction r g (the growth respiration coefficient) of the gross primary productivity ( G ) minus the plant maintenance respiration: Note that this relation results in the correct growth respiration on timescales of the order of a day or longer (on the model time step scale, R pg will be negative in the night, which is misleading if taken in isolation). The net primary productivity N is therefore = (1 − r g )( G − R pm ).

A6 Irrigation
In JULES, irrigation is implemented such that the water in the top two soil layers is continuously topped up to a critical level (often the field capacity) during the "irrigation season", if sufficient irrigation water is available. We will consider the irrigation season to last all year (irr_crop=0) and treat the supply of irrigation as unlimited (l_irrig_limit=F). With these settings, the soil water stress parameter β stays approximately equal to 1; i.e. the plant is not water stressed. When irrigation is on, the root distribution has a negligible influence on model performance.

A7 Nitrogen limitation
Although JULES has a nitrogen cycle implemented (as of version 4.4), it can not yet be used in conjunction with the crop model. We therefore make the assumption here that the crops are not nitrogen limited.