JULES-crop: a parametrisation of crops in the Joint UK Land Environment Simulator

. Studies of climate change impacts on the terrestrial biosphere have been completed without recognition of the integrated nature of the biosphere. Improved assessment of the impacts of climate change on food and water security requires the development and use of models not only representing each component but also their interactions. To meet this requirement the Joint UK Land Environment Simulator (JULES) land surface model has been modiﬁed to include a generic parametrisation of annual crops. The new model, JULES-crop, is described and evaluation at global and site levels for the four globally important crops; wheat, soybean, maize and rice. JULES-crop demonstrates skill in simulating the inter-annual variations of yield for maize and soybean at the global and country levels, and for wheat for major spring wheat producing countries. The impact of the new parametrisation, compared to the standard conﬁguration, on the simulation of surface heat ﬂuxes is largely an alteration of the partitioning between latent and sensible heat ﬂuxes during the later part of the growing season. Further evaluation at the site level shows the model captures the seasonality of leaf area index, gross primary production and canopy height better than in the standard JULES. However, this does not lead to an improvement in the simulation of sensible and latent heat ﬂuxes. The performance of JULES-crop from both an Earth system and crop yield model perspective is encouraging. However, more effort is to develop the parametrisation of the model for speciﬁc applications.


Introduction
Understanding how climate variability and change will impact upon crop production systems is a research challenge of utmost importance to society. To date, studies of climate change impacts on the terrestrial biosphere have been completed without recognition of the integrated nature of the biosphere. Crop simulation models are widely utilised as they incorporate many known effects of how changes in atmospheric conditions can impact upon crop growth, development and yield. However, they do not simulate the wider interactions of crops and the environment. For example, climate change will impact upon water resources which will in turn impact upon the water available for irrigation of crops. Betts et al. (2015) used the Hadley Centre Earth System Model (HadGEM2-ES) to evaluate climate impacts on the terrestrial biosphere under a range of emission scenarios. By doing so they were able to assess several elements of the terrestrial system in a way that was fully integrated and consistent with the climate projections. However, they were only able to include natural systems as crops are not yet included in the model. Including a representation of crops within land surface models will facilitate a more comprehensive, integrated and internally consistent simulation of the impacts of climate change and variability on the full land system, ac-counting for interactions between different components and processes. This will ultimately enable improved projections of the impacts of climate change on food and water security, including interactions between the two. There is increasing evidence that the cultivation of crops affects weather and climate on local scales. Croplands now occupy 12 % of Earth's ice-free land surface and in several regions of the world are the dominant vegetation type on the land surface (e.g. midwest USA, Indo-Gangetic Plain). This extensification of agriculture has altered the biophysical characteristics of the land surface potentially altering regional climate. Therefore, there is reasoning to consider crops and climate as a truly coupled system and hence motivation to develop models which can fully represent the coupled feedbacks between them.
Efforts to simulate the environmental impacts on crop production are commonly thought to have begun in the 1960s at Wageningen (van Ittersum et al., 2003). Since then crop modelling has grown and there are now many models available in the research and agronomic domains. Such models have been deployed both as decision support tools and to research the impacts of climate change on future crop production. Recent advances in crop modelling include the application of crop models, traditionally developed at the field level, to cover the globe on a gridded basis (Deryng et al., 2011;Osborne et al., 2013) and inter-comparison of many crop models in simulating the same crop and the same set of conditions (Asseng et al., 2013).
The investigation of how croplands affect weather and climate is much less mature. The initial expansion of cropland area came at the expense of forests and the impact of this deforestation has received considerable research attention. However, croplands have also replaced more similar native grasslands. For example, McPherson et al. (2004) showed that the near-surface climate over the now intensively cultivated winter wheat belt in Oklahoma, USA, is significantly different to that over adjacent grasslands. McPherson et al. (2004) identify the differences in phenology between managed croplands and natural grasslands as the determinant of the differences.
The increase in understanding of how croplands might differentially impact the climate compared to natural vegetation has led to a recent surge in model development whereby land surface or global vegetation models have been extended to include explicit parametrisations of crops, in place of the use of grasslands as a surrogate (see review of Levis, 2010). Some developments have been motivated by improving the carbon and water budget of land surface modelling (Bondeau et al., 2007), others to include croplands in global or regional climate models to better represent their impact on the atmosphere (Lokupitiya et al., 2009;Chen and Xie, 2012;Levis et al., 2012), while others have been motivated to consistently simulate both yield and environmental impacts (Kucharik and Brye, 2003).
The aim of this model development was to develop a combined land surface and crop model capable of simulating both the impacts of climate variability on crop productivity, as well as the impact of croplands on the climate. To achieve this we have added a crop-specific parametrisation to the Joint UK Land Environment Land Surface (JULES) land surface model. JULES is the land surface scheme of the UK Met Office Unified Model and the next generation UK Earth System Model (UKESM) and, therefore, can be in time coupled to a state-of-the-art climate model. A full description of JULES can be found in Best et al. (2011) andClark et al. (2011). JULES does not currently include an explicit parametrisation of crops; instead, over cropped regions, the C 3 or C 4 grass plant functional types are used. Previous work has included crops in the model. Osborne et al. (2007) included a crop parametrisation in MOSES (i.e. in the fully coupled land surface-climate model) based on the groundnut version of the crop model GLAM. More recently, Van den Hoof et al. (2011) extended JULES to include a parametrisation of wheat based on the crop model SUCROS. Neither Osborne et al. (2007) nor Van den Hoof et al. (2011) developed a generic representation of crops suitable for the examination of different crops throughout the globe, something that is important from an Earth system modelling perspective. Therefore, the objective of this study was to develop a generic parametrisation of crops applicable to many crop types and at the global scale. However, the model has been designed to be flexible, meaning users can reparametrise the model depending on requirements (e.g. to represent different crop cultivars).
The following section describes the model development, Sects. 3 and 4 present an evaluation of the new model when applied at global and site levels, respectively, followed by a Discussion (Sect. 5).

Model description
The essence of JULES-crop is illustrated in Fig. 1. The additional model equations required to simulate crops essentially partition the carbon uptake of vegetation already simulated by JULES in to several crop organs and the size of the crop, important for land surface-atmosphere feedbacks, is derived from the organ biomass using allometric equations. The pattern of partitioning of assimilated carbon to the crop organs is affected by the crop development rate which itself is influenced by temperature. In addition to the new equations describing crop growth and development, changes to the model structure were also required to accommodate the additional plant functional types. New equations describing crop growth and development were added to the model. Each crop is considered as an additional plant functional type and a distinction is made between natural and crop plant functional types within the model, with the crop plant functional types requiring extra parameters to be specified. The detailed description of the crop parametrisation is split in to three parts. Firstly, the equations that determine the start and dura-  tion of the crop growing season are described. Secondly, the equations determining the rate of crop growth are described. Lastly, the changes to model structure are outlined. A full listing of new model parameters and variables can be found in Tables 1 and 2, respectively.

Growing season and development
The crop growing season begins when the crop is sown. This date can either be prescribed (i.e. if it is known) or calculated dynamically based on environmental criteria. In the latter case, sowing only occurs when the soil is wet enough (θ 2 >θ c,2 , where θ 2 is the soil moisture content in the second layer and θ c,2 is the critical soil moisture content in the sec-ond layer), which is warm enough (T soil,3 >T b + 2 K, where T soil,3 is the temperature in the third soil layer and T b is the base temperature), and when days are not rapidly shortening (dP /dt> − 0.02 h d −1 , where P is the day length). We wish to make users aware of this sowing option; however, we feel it needs further optimising and so results using the dynamic sowing date will not be included here. The use of subsurface soil moisture and temperature variables prevents sowing occurring too early in response to short-term fluctuations in weather. The rate of day length criterium ensures that crops are not sown too late in the year when conditions for growth are deteriorating. Once sown, the crop develops through three stages: sowing to emergence, emergence to flowering, and flowering to maturity. Harvest is assumed to occur at crop maturity. The rate of crop development is related to thermal time. Given the 1.5 m tile temperature (T ), an effective temperature (T eff ) is calculated based upon the crop-specific cardinal temperatures (T b , T o , T m -see Table 1 for description).
T eff is greatest and hence development is fastest at T = T o . As temperature falls below or rises above T o the rate of development linearly decreases until no development occurs when either T ≤ T b or T ≥ T m . For the sowing to emergence phase, T eff is not affected by T m or T o (i.e. T eff = T −T b ). This equa-  net primary productivity A c kg C m −2 Eq. (4) net carbon assimilation R dc kg C m −2 Eq. (4) canopy dark respiration tion is a "standard" way of calculating effective temperature (Challinor et al., 2004). An important difference to other available models is that JULES-crop simulates a decline of T eff above the maximum temperature, whereas others keep T eff at the maximum value no matter how high temperatures get. For some crops, progress towards flowering is slowed if the day length (P ) is less than (greater than) a cropspecific critical photoperiod (P crit ) for long-day (short-day) crop types. The degree of sensitivity to the photoperiod is represented by the parameter P sens which is positive for short-day plants and negative for long-day plants. This conceptual approach was motivated by Loomis (1992). Therefore, to slow development T eff is multiplied by the relative photoperiod effect (RPE), which is defined as follows: The status of crop development is represented by the development index (DVI) which takes the value of −1 upon sowing, increasing to 0 on emergence, 1 at the end of vegetative stage and 2 at crop maturity. The rate of increase of DVI is calculated as follows, where TT emr is the thermal time between sowing and emergence, TT veg is the thermal time between emergence and flowering and TT rep is the thermal time between flowering and harvest: The growing season ends when DVI = 2 at which time the prognostic variables related to crop growth (L, h, C root , C harv , C resv ) are reset to minimal values close to 0. To prevent growing seasons continuing indefinitely when conditions are no longer suitable, the crop is also harvested if the soil temperature in the second soil layer falls below T b at any time after DVI = 1 or if LAI > 15 (leaf area index). Vernalisation, a cold temperature requirement for development in some crops, is not included in this model version.

Growth
To simulate crop growth, net primary productivity ( ) is accumulated over a day and then partitioned between five carbon pools: root (C root ), structural stem (C stem ), stem reserves (C resv ), leaves (C leaf ), and harvested organs (C harv ). The original formulation for in JULES includes assumptions about the sizes of the leaf, stem and root carbon pools in order to estimate respiration loses. Stem carbon is a function of leaf area index (Eq. 42 of Clark et al., 2011) and root carbon is set to equal leaf carbon. Because these carbon pools are now explicitly simulated, is recalculated for the crop types with the following equation based on an algebraic reduction of the set of equations used in JULES: where r g is the fraction of gross primary productivity less maintenance respiration that is assigned to growth respiration, A c is the net canopy photosynthesis, and R dc is the rate of non-moisture-stressed canopy dark respiration. C leaf , C stem and C root are the carbon content of leaf, stem and root, respectively. The carbon in is accumulated over a day and then divided into five crop components according to "partition coefficients", one for each of the four root, stem, leaf and harvest pools defined above and a reserve pool. These components are added to the (state variable) pools of carbon describing the crop.
where τ is the fraction of stem carbon that is partitioned in to the reserve pool. p root + p leaf + p stem + p harv = 1.0. Partition coefficients for a given crop are typically predefined in process-based crop models according to either the length of time since emergence or to crop development stage (DVI; i.e. a function of thermal time since emergence). They are represented by fixed values for a given period of time (or thermal time) since emergence, and these values are listed in a look-up table and referenced for each iteration of the model (e.g. WOFOST, van Ittersum et al., 2003).
Here we define the partition coefficients as a function of thermal time using six parameters to describe continuously varying partition coefficients over the duration of the crop cycle. We use a multinomial logistic to define this function: p root = e α root +(β root DVI) e α root +(β root DVI) + e α stem +(β stem DVI) + e α leaf +(β leaf DVI) + 1 , p stem = e α stem +(β stem DVI) e α root +(β root DVI) + e α stem +(β stem DVI) + e α leaf +(β leaf DVI) + 1 , where α and β are empirically derived parameters describing the shape of the thermal time-varying partition coefficient for leaves, roots and stems, and DVI is the development index. Thus, for only six parameters (which is also the absolute minimum number of parameters needed to define partition coefficients for four carbon pools) we can define a much wider range of shapes of thermal time varying partition coefficients. Furthermore, these six parameters can be more feasibly calibrated than a larger number of "look-up" partition coefficients. This parametrisation is illustrated in Fig. 2 overlaid with example observed partitioning fractions from de Vries et al. (1989). Following the formulation of de Vries et al. (1989), once carbon is no longer partitioned to stems, carbon from the stem reserve pool is mobilised to the harvest pool at a rate of 10 % a day: Leaf senescence is treated simplistically by mobilising carbon from the leaf to the harvest pool at a rate of 0.05 d −1 once DVI has reached 1.5. This equation was inspired by Eq. (7), but based the period for which senescence starts on a specific DVI value (1.5) rather than waiting for partitioning of leaves to cease since for some crop types this does not happen.
At the end of each growth time step (24 h), the amount of carbon in the leaves is related to leaf area index (L) by where SLA = γ (DVI + 0.06) δ .
The values of γ and δ were determined by fitting the relationship to the paired values of DVI and SLA (specific leaf area) reported in de Vries et al. (1989). The amount of carbon in the stem is related to the crop height by (Hunt, 1990) The values of κ and λ were determined by fitting the relationship to the paired values of h and C stem at the Mead FLUXNET site (Verma et al., 2005).
Equations (9) and (11) are rearranged to derive the carbon content of leaves and stems, respectively, before each growth time step.
Because root biomass increases during the crop growing season the fraction of roots in each JULES soil layer varies according to the equation of Arora and Boer (2003) which defines the fraction of roots at depth z as where where d r is 0.5 for all crop types, and r dir is a crop-specific parameter.
To ensure crop establishment, the growing season is curtailed if the sum of root, leaf, stem and reserve carbon falls below the initial seed carbon content (or zero carbon content) if the sowing date is determined dynamically.

Changes to JULES code structure
The standard version of JULES represents the land surface as a combination of up to nine surface types including five plant functional types: broadleaf trees, needleleaf trees, C 3 grass, C 4 grass, shrubs, bare soil, inland lakes, snow and ice. Surface fluxes of heat, moisture and momentum are determined independently for each tile before being combined to a single set of fluxes according to the relative fractions of each tile. Each crop type is considered as a different tile. Therefore, it is possible to simulate many crops or crop varieties at a site or grid box in a single integration of JULES, in addition to the standard five plant functional types. The parameters required to represent vegetation within JULES were extended to the crop tile(s). The values were copied across from the JULES default parameters for C 3 and C 4 grass, depending on the crop photosynthetic capacity (see Table 3).
The values of the parameters required in Eqs.
(1)-(13) determine which crops are being simulated and can be varied according to different user requirements, e.g. crop species (e.g. maize or wheat), generic crop type (e.g. C 3 cereals) or cultivar (e.g. soybean PS123121 or soybean 21h321). Each parameter is described in Table 1. Values for each parameter can be determined by calibration against relevant observational data such as leaf area index, biomass, and yield from agricultural field stations. For this study such an exercise was not performed. Instead, suitable values were determined from either the literature or by tuning to fit site-level data in order to establish a model version that could be evaluated at site and global scales.  Table 4 and were either taken from the crop science literature or calibrated as described below. Specifically, the values for the partition parameters α root, stem, leaf and β root, stem, leaf and the specific leaf area coefficients γ and δ were calibrated against data in de Vries et al. (1989). The allometric coefficients κ and λ were determined by calibration against paired crop height and stem biomass data from FLUXNET sites. The cardinal temperatures (T b , T o , and T m ) were specified values in line with the range of values reported in the literature (see Gawith, 1999, andSanchez et al., 2014). The effect of photoperiod was not included (by setting P crit to 24) due to our method of determining thermal time between emerging and flowering (TT veg ) and thermal time between flowering and harvest (TT rep ) (see below).
The parameter r dir was set to 0 for all crop types, which effectively removes the effect of increasing root carbon on the vertical distribution. Early tests of the model revealed that including an effect of increasing root carbon led to high levels of water stress at the start of the crop growing season lead- ing to poor crop growth. Therefore, the effect was "turned off". The parametrisation was left in the model to allow other model users to experiment further with dynamic root growth. The global model runs were driven by the CRU-NCEP v4 climate data extended to include 2012 (N. Viovy, personal communication, 2013) as used by the Global Carbon Project (Le Quéré et al., 2013). This was regridded to a N96 grid (1.875 • longitude × 1.25 • latitude) and used with ancillaries from HadGEM2-ES (Collins et al., 2011;Jones et al., 2011) to evaluate the performance of the model in a Earth system model set-up. A multi-layer canopy radiation scheme was used, accounting for direct/diffuse radiation components including sunflecks (can_ran_mod = 5). The main run was from 1960 to 2010 and the spin-up consisted of repeating the first 10 years 5 times. The sowing dates were taken from Sacks et al. (2010), and a value for each land grid box was obtained using nearest-neighbour extrapolation. The values of TT veg and TT rep were allowed to vary spatially and determined such that, when used with the CRU-NCEP tempera- (TT veg +TT rep ) = 0.5, 0.45, 0.6, 0.6 for soybean, maize, wheat, and rice, respectively. Photoperiod sensitivity was not considered.This is because including it would have made calculating TT veg and TT rep almost impossible, because three variables would need calibrating at each grid cell (total TT, critical photoperiod, and sensitivity to photoperiod) from one observation (growing season duration). For comparison a control run was completed using the same model set-up but with the crop code switched off. This run is used to assess performance against the standard land surface scheme in the Met Office Hadley Centre Earth System Model -HadGEM2-ES. Figure 3 shows the planting date of Sacks et al. (2010) and the derived maps of TT veg and TT rep . Sacks et al. (2010) derived gridded planting dates from national-or district-levelreported planting dates which are given in months rather than days. Therefore, there is little spatial or temporal variation in the sowing date which might well be expected due to variations in local climate and management practices. However, the data serves a purpose in global modelling studies by providing an approximate start point for the growing season at the right time of year. Our method of calculating the crop thermal time requirements produces considerable spatial variability which is determined in reality by variation in study) to the annual accumulated thermal time and then used that relationship to determine thermal time requirements under future climate. The approach in this study was chosen as the simplest and most likely to achieve growing seasons of lengths close to observed. Due to the absence of a vernalisation parametrisation in the model only spring wheat was considered. The crop fractions were taken from Monfreda et al. (2008) and regridded to the N96 HadGEM2-ES resolution. Monfreda et al. (2008) provide observations in the year 2000 which were used to describe the crop coverages for the whole integration period due to a lack of available data sets covering this time period.

Evaluation
The simulated grid box annual yield for each crop averaged over the 50 years is shown in Fig. 4 alongside global gridded observations for circa 2000 (Monfreda et al., 2008). Figure 4 shows that in general the model is underestimating yields in arid, irrigated regions and overestimating them in tropical regions. In particular, simulated maize yields are significantly larger than observations in tropical regions. Given that the model does not include any information on the yield gap (the difference between actual farm-level yield and potential yield) or important land management such as irrigation the spatial variability of model output should not be too closely compared to that of observed yield. Instead, a greater appreciation of model performance can be gained from examining the year to year fluctuations in yield, given that the effects of changes in management and technology materialise over several years. Figures 5 and 6 show the simulated global-and countrylevel yields for wheat, soybean, maize and rice between 1960   and 2008 compared to the reported yields of FAO (2014). Simulated global yield was determined by multiplying the simulated annual maximum yield at each grid cell by the observed harvested area from Monfreda et al. (2008) regridded to the HadGEM2-ES spatial resolution. This grid cell estimate of production was summed over all grid cells to produce an estimate of global production which was then divided by the total harvested area to provide an estimate of global yield. Grid cell yields were determined from the annual maximum value of C harv which was multiplied by 2 to convert from carbon mass to total biomass, by 1.16 to account for grain moisture content, and by 10 to convert from kilograms per squared metre to megagrams per hectare. Not all grid cells were included in the analysis. Cells were excluded if the annual maximum DVI was less than 1.5 which was possible if the growing season was curtailed if LAI > 15 or t soil,2 <T bse . A similar analysis was conducted to determine country-level yields with averages taken over all grid cells within a particular country. Country yield observations were de-trended for comparison with model output. This is because the increasing trend in yield observations over the last 50 years is due to improvements in agricultural technology and management and increasing atmospheric carbon dioxide. This trend is not reproduced by the model as these processes were not represented in our set-up.
The average simulated yield for maize is overestimated; however, the model does a reasonable job of reproducing the inter-annual variability at the global (r = 0.48) and country scales (Fig. 6a), although there is a tendency to simulate larger variability than observed. For soybean, average yield is again much greater than observed but year on year variability is correlated with observations (r = 0.37) providing some confidence in the model's ability to simulate the observed response of soybean yield to climate. Regionally, in countries such as the USA (r = 0.39) and India (r = 0.52) JULES-crop is able to reasonably capture inter-annual variability of yields (Fig. 6b). For rice, yield levels are higher than reported, variability is overestimated and not correlated with observations (r = 0.24). At the country level, model simulations in India (r = 0.57) correlate with observations (Fig. 6c). The average simulated yield level for wheat is similar to the most recent observations; however, when comparing the year to year fluctuations in yield, the correlation between simulated and observed yields is low (r = 0.019). Because JULES-crop only simulates spring wheat then the comparison to reported wheat yields is slightly unfair given that the majority of wheat produced globally is from winter varieties. It is encouraging that the best agreement between simulated and observed yield fluctuations at the national level is for Turkey (r = 0.46) and Australia (r = 0.53), in which spring wheat varieties dominate.
For all crops there is a tendency for JULES-crop to simulate larger variability than observed. This may in part be explained by the lack of certain processes in the model (particularly those to do with land management). For example, not including a representation of irrigation in the model may explain why the model predicts lower yields than observations as irrigation would act to reduce the extent of crop failure in drought years. The model also does not include the impacts of pests and disease which may reduce overall yields in some years. Importantly, the model does not as yet include a nitrogen cycle which may reduce overall GPP (gross primary production), bringing the simulations in line with observations. This may also explain why there are strong deviations between the magnitude of observed and modelled yield in tropical countries where climatic conditions for growth are good in the absence of the limitations described above.
For some countries simulations of yield capture the magnitude and variability of observations. In other countries the model reproduces the variability in yield but over-predicts the magnitude. There are also countries where the model performs poorly in simulating both variability and magnitude. This variety of results is due in part to the use of generic parametrizations for global model runs which is a limitation of this type of Earth system model set-up. By using parameters that do not vary spatially we can not fully represent the range of crop varieties that are found globally.
To evaluate the impact of including the crop parametrisation on JULES, output from the simulation with crops included is compared to a control simulation of the standard JULES configuration with grass plant functional types taking the land fraction of crops. Impacts on the land surface will be mostly mediated via direct changes to the vegetation structure and also via indirect effects on state variables, most obviously the soil moisture content. To begin to examine the potential for impact, the changes to a key vegetation variable LAI values are shown in Fig. 7 for four major crop producing countries. To produce the country averages, grid cell LAI are combined by weighting by the grid cell contribution to total country crop area. In the USA and China each crop growing season occupies the similar set of summer months, whereas for India and Brazil the wheat cropping season is distinct from the other three crops. Peak LAI is greatest in Brazil and lowest in China, which is most likely a reflection of the absence of irrigation in the model and the relative abundance of rainfall in each country. In comparison to the standard JULES configuration the addition of crops adds a seasonality to LAI as there is no default seasonality to vegetation characteristics in JULES. The annual variation of crop LAI is dampened when aggregated with the other plant functional types, which explains the non-zero LAI in the non-growing season in the JULES-crop simulation. Figure 7 shows that the inclusion of crops alters the grid box net primary production (NPP) in terms of the timing of peak fluxes. There are also lower fluxes in winter due to the more realistic treatment of LAI at this time. Therefore, including a representation of crops in JULES may help improve the seasonality of LAI, which affects carbon fluxes. Figure 8 shows that the impact of these differences in vegetation size during the year is greatest for the surface moisture flux and sensible heat flux rather than the components of the radiation balance. The largest impacts are on the sensible heat flux towards the end of the crop growing season which is higher with the inclusion of crops. For India, there is a concomitant decrease in the surface moisture flux, implying that the total available energy at the surface is unaltered but is partitioned differently between sensible and latent heat fluxes. The impact of JULES-crop on the energy balance is however minimal. In this configuration the model is forced by prescribed meteorology at screen height. This has the tendency to dampen the model in comparison to a full atmospheric simulation in which the boundary layer state is able to evolve. It may, therefore, be expected that a GCM (global climate model) may be more sensitive to changes in the surface state.

Model set-up
To further understand the impact of adding crops to JULES, site-level simulations were also performed. Sites were selected by the vegetation cover (only croplands were considered) and by the availability of meteorological and biological data required to force the model and evaluate model results. The sites selected were are all in the USA: Mead in Nebraska (Verma et al., 2005) and Bondville and Fermi in Illinois. For each site, three simulations were performed: the standard configuration of JULES, standard JULES with the existing phenology parametrisation turned on, and the full JULES-crop parametrisation. For the JULES-crop simulation, the fractional coverage of the relevant crop type was set to 1 with all other functional types set to 0. For the JULES (non-crop) simulations, the fractional coverage of the relevant grass functional type (i.e. C 3 grass for soybean, C 4 for maize) was set to 1. All crop parameters were prescribed the same value as in the global simulations. The sowing date and thermal time requirements were taken from the relevant grid cell for each site.

Evaluation
Figures 9 and 10 compare JULES-crop simulations for the soybean crop type with standard JULES C 3 grass plant functional type with and without phenology, and with observations where available. The crop parametrisation captures the evolution of LAI and canopy height across the season, although the model underestimates these growth variables. The model also simulates lower GPP fluxes compared to observations; however, crop yields are comparable. The standard C 3 grass with phenology configuration of JULES also simu-lates growth and decay of vegetation cover but over a longer period of time than the observed growing season. Without the phenology routine the LAI is set to the default for C 3 grass of 2.0 all year. Interestingly, the more realistic simulation of vegetation cover does not lead to improved simulation of surface fluxes. At all sites similar characteristics of the simulations are evident. During winter all three configurations simulate similar latent and sensible heat fluxes in line with observations (Fig. 10). Towards the start of the growing season the standard configuration of JULES with constant LAI = 2.0 overestimates latent heat flux due to an unrealistically large vegetation coverage. The simulations with phenology and crops have lower vegetation cover and simulate lower latent heat flux but are still noticeably greater than observations. At around the peak of crop cover, all simulations underestimate the latent heat flux and overestimate the sensible heat flux due to lower simulated LAI compared to observations.
Site-level simulations for the maize crop type are shown in Figs. 11 and 12. The crop parametrisation is reasonably successful in capturing the LAI and canopy height of maize at all evaluation sites. Similarly, GPP and yields are lower than observed although the seasonal pattern of GPP is close to observations. Overall, model simulations broadly capture the patterns of latent and sensible heat fluxes although, again, there are no major improvements in model performance with the explicit inclusion of crops. At Fermi, in 2006, the cropspecific simulation captures the observed evolution of LAI reasonably well with peak LAI slightly closer to observations than the standard JULES simulations. However, this, again, does not improve the simulation of heat fluxes.
All model configurations overestimate the partitioning of energy in to latent heat before the growing season begins and underestimate it during the crop growing season, despite widely varying LAI values. This could be due to the relatively weak LAI-surface conductance relationship found in JULES (Lawrence and Slingo, 2004). This is reflected in the low sensitivity to LAI between fixed phenology and seasonal phenology. In these simulations we would therefore not expect a large response to an alternative representation of crop LAI phenology. This comparison serves as a reminder that improving the realism of a model may not guarantee im-proved performance in the model in other aspects. The results also show that JULES (crop and standard configurations) is not able to capture the magnitude of observed GPP fluxes. This suggests that using the standard physiological parameters for C 3 and C 4 grasses is not appropriate when representing crops particularly as JULES does not include nitrogen fertilisation explicitly. Tuning of parameters that describe leaf nitrogen, for example, may improve fluxes of GPP and hence overall yields. It is also worth noting that the parameters used for the crop model in the site simulations are from the global set-up and hence probably not optimal for site simulations.

Discussion and conclusions
When designing JULES-crop we took a flexible approach in acknowledgement of the different requirements of the science community. This means the model can be used to address a range of science questions, for example, (a) to assess global climate impacts on crop functional types over long integrations with climate model output, (b) to represent a num- ber of crop cultivars of the same crop type at the site scale forced with weather observations and (c) to assess how crops may impact on biogeophysical feedbacks to climate including albedo, partitioning of turbulent fluxes and seasonality of LAI. In this paper we present results from a generic, crop functional type parametrisation implemented at both global and site scales to show how this model performs in an Earth system model context. Having the aim of generality necessarily means that the model loses out in terms of specificity. However, with further effort it should be possible to tailor the model set-up for more specific applications but with the requirement that attention is given to the choice of parameter values. Default values are provided here as a starting point for model development and initial evaluation.
These results demonstrate the importance of evaluating the performance of JULES-crop in a holistic sense, assessing both its ability to simulate land surface fluxes in addition to crop growth and development dynamics and to recognise that identified biases in performance are the result of the combined JULES-crop model, not just the added crop component. Adding a crop parametrisation has increased the complexity of JULES. However, this has not led to an immediate improvement in the model's simulation of surface fluxes, at least at the measurement sites examined. More effort needs to go into developing the parameter sets for crops within JULES, particularly for the existing set of plant functional type parameters which control productivity. Comparing the regional patterns of yield to observations gives useful insight into the existing limits of the model. It is clear that some important processes are missing, particularly irrigation (although this model development will shortly be submitted for release). Developing a nitrogen cycle for JULES (model development also in progress) should also improve the model simulations, as introducing nitrogen limitation has been shown to reduce overall productivity in Earth system models (Thornton et al., 2009). JULES-crop will still exclude many management factors which affect regional yields. Licker et al. (2010) estimated global yield gaps and showed they were greatest in tropical regions. Although not directly comparable with our simulations, this study shows us that JULES-crop simulations are likely to overestimate yields in tropical regions compared with observations. How- ever, we have deliberately not introduced a yield gap adjustment as it would not be physically based and as such would be difficult to apply to future simulations. It is, however, important to capture regional differences due to management as they will effect patterns in productivity and hence feedbacks to the climate. In an Earth system model context it is better to represent these management processes explicitly where possible, as they effect not only crop growth but also may well influence the local climate directly (e.g. irrigation; Sacks et al., 2009).
As a yield simulation model, there are encouraging signs that JULES-crop can simulate variability in yield associated with climate fluctuations. However, it is clear that JULEScrop overestimates the magnitude of this variability. Whilst the absence of irrigation is most likely a contributing factor to the overestimation of yield variability, the implication that the model is too sensitive to changes in environmental conditions should also be investigated further.
Including crops in JULES gives a more realistic seasonal cycle of leaf area index, which affects the seasonality of carbon fluxes (timing of peak flux and lower winter fluxes). This was seen at both the global and site levels. The impact of crops on the energy balance was to alter the partitioning of latent and sensible heat fluxes particularly in winter, which led to small impacts on temperature in some countries. These impacts were marginal at the country and site scales despite quite large differences in LAI. It is possible that the relationship between LAI and evaporation is too weak in JULES (Lawrence and Slingo, 2004), which may explain why a more realistic representation of LAI did not improve the energy fluxes. We may expect a higher sensitivity in fully coupled atmosphere models.
Crop production systems are by their very nature heavily influenced by humans. This represents a challenge to the JULES model which, to date, assumed vegetation to be static and, within each vegetation tile, homogeneous by the use of global constants for parameter values. The level to which this approach can be extended to crops is limited. Whilst some processes might be considered fundamental (i.e. photosynthesis) others can vary from place to place for the same crop (e.g. sensitivity of development rate to day length). Furthermore, human interference can alter the fundamental process, for example the application of fertiliser to increase leaf nitrogen contents impacting on photosynthesis. For applications of JULES-crop that rely on accurate yield simulations, the inclusion of either a yield gap variable or the factors that determine it such as fertiliser applications, pest control, and soil fertility should be a priority for future model development. Inclusion of winter wheat is also of high priority for JULES-crop. This is important for use of JULES-crop as a yield simulation model but also an Earth system model, as the additional presence of vegetation cover from autumn to spring would impact on surface characteristics (albedo, heat capacity, etc.).