Improving the representation of sugarcane crop in the Joint UK Land Environment Simulator (JULES) model for climate impact assessment

Bioenergy from sugarcane production is considered a key mitigation strategy for global warming. Improving its representation in land surface models is important to further understand the interactions between climate and bioenergy production, supporting accurate climate projections and decision‐making. This study aimed to calibrate and evaluate the Joint UK Land Environment Simulator (JULES) for climate impact assessments in sugarcane. A dataset composed of soil moisture, water and carbon fluxes and biomass measurements from field experiments across Brazil was used to calibrate and evaluate JULES‐crop and JULES‐BE parametrizations. The ability to predict the spatiotemporal variability of sugarcane yields and the impact of climate change was also tested at five Brazilian microregions. Parameters related to sugarcane allometry, crop growth and development were derived and tested for JULES‐crop and JULES‐BE, together with the response to atmospheric carbon dioxide, temperature and low‐water availability (CTW‐response). Both parametrizations showed comparable performance to other sugarcane dynamic models, with a root mean squared error of 6.75 and 6.05 t ha−1 for stalk dry matter for JULES‐crop and JULES‐BE, respectively. The parametrizations were also able to replicate the average yield patterns observed in the five microregions over 30 years when the yield gap factors were taken into account, with the correlation (r) between simulated and observed interannual variability of yields ranging from 0.33 to 0.56. An overall small positive trend was found in future yield projections of sugarcane using the new calibrations, with exception of the Jataí mesoregion which showed a clear negative trend for the SSP5 scenario from the years 2070 to 2100. Our simulations showed that an abrupt negative impact on sugarcane yields is expected if daytime temperatures above 35°C become more frequent. The newly calibrated version of JULES can be applied regionally and globally to help understand the interactions between climate and bioenergy production.


| INTRODUCTION
The use of bioenergy with carbon capture and storage is considered an essential mitigation strategy to limit global warming to below 2°C (van Vuuren et al., 2018). However, improving the representation of energy crops in biophysical models is still needed to further understand the interplay between climate and bioenergy production (Hanssen et al., 2020). While food crops have multiple and well-detailed models delivering robust climate impact assessments , bioenergy crop models still lag in number and detail (Surendran Nair et al., 2012). Yet, it is difficult to capture the bidirectional feedback between croplands and the atmosphere (Betts, 2005), leading to high uncertainty about the effects of large land-use change scenarios on climate and ecosystems (Hanssen et al., 2020).
Land surface models (LSMs) can be coupled to general circulation models (GCMs), providing an opportunity to simulate the impacts of climate on crop productivity while considering the feedback of croplands on atmospheric processes (Pitman, 2003). Initially developed to represent fluxes of carbon, water and energy at the landatmosphere interface, they have a wide range of applications in the fields of biogeochemical, hydrological and energy cycles at multi-scales, including in agriculture (Osborne et al., 2015). Although using sophisticated methods to solve the energy balance and carbon fluxes at multiple scales, the representation of bioenergy crops in LSMs is often limited to a broader generalization of vegetation based on crop functional type/plant functional type (PFT; Surendran Nair et al., 2012).
Sugarcane (Saccharum officinarum) is a key bioenergy crop with significant social, economic and environmental importance in many developing countries (Marin et al., 2016). Brazil is the largest producer accounting for ca. 40% of global stalk fresh mass production in 8.6 million hectares (IBGE, 2019). It is a strategic crop for bioenergy production and mitigation of climate change, where bioethanol and electricity from sugarcane biomass comprise a significant share of the energy matrix in Brazil. In this context, many process-based models were developed to support decision-making and strategic planning by simulating the complex interactions between Genotype × En vironment × Management (Dias & Inman-Bamber, 2020; Marin et al., 2015). Although widely tested and calibrated, the application of these models in large-scale modelling frameworks is still challenging because they were mainly developed to operate at the field scale.
In contrast, very few studies were dedicated to simulate sugarcane growth using LSMs, where only three models, namely Agro-IBIS , LPJmL (Lapola et al., 2009) and JULES (Black et al., 2012) had their performance assessed against regional yield data. Except for the Agro-IBIS model , the lack of data about sugarcane growth and energy fluxes was acknowledged as the major limiting factor to assess the biophysical consistency of simulations. While large-scale data can be used to evaluate LSMs results, many processes simulated within those models can only by validated at the local field scale (e.g. layered-canopy photosynthesis, carbon allocation to organ pools and senescence). Therefore, measurements from high monitored field experiments provide the appropriate data to isolate, calibrate and test the subroutines that solve these processes, as recently demonstrated for the JULES model with crops such as maize (Williams et al., 2017), soybean (Leung et al., 2020) and Miscanthus (Littleton et al., 2020).
The Joint UK Land Environment Simulator (JULES) is a community LSM that is applied both as a standalone model and as the land surface component in the UK Met Office's Unified Model Clark et al., 2011). JULES simulations are also part of the Coupled Model Intercomparison Project (CMIP), a foundational element for climate and earth systems sciences over the last years (Eyring et al., 2016). As a community model, JULES is freely available and it is in constant development where two parametrizations were recently developed to explicitly simulate crop growth and development, namely JULEScrop (Osborne et al., 2015) and JULES-BE (Littleton et al., 2020).
JULES-crop was initially presented by Osborne et al. (2015) and evaluated for wheat, maize, soybean and rice. Williams et al. (2017) conducted a comprehensive calibration and evaluation of JULES-crop for the maize crop and found the need for more studies focused on soil moisture stress, carboxylation rates and crop thermal requirements when using JULES-crop for C 4 species. The JULES-BE was developed to represent bioenergy crops and continuous harvesting schemes and was tested for | 1099 VIANNA et al. Miscanthus (Littleton et al., 2020). Although showing good performance to simulate water and carbon fluxes at one site in the UK, global observed yields of Miscanthus were far more variable than captured by the model, requiring further evaluation across different environments and crops to identify limitations and opportunities for improvement.
In this study, we aimed to calibrate and evaluate the JULES model, and its crop parametrizations JULES-crop and JULES-BE for improved representation of the sugarcane crop in climate impact assessments. A dataset of biomass, carbon fluxes, soil moisture and evapotranspiration (ET) collected across high monitored field experiments in Brazil was used to assist in calibration, statistical evaluation and biophysical consistency of the new parameters. As JULES is mainly applied to large-scale climate studies, we also assessed the performance of the new calibration in predicting the interannual yield variability observed in five contrasting regions where sugarcane was cultivated for 30 years in Brazil. Finally, we extended the 30-year simulations with future climate projections to understand the crop response to climate change scenarios using the new calibrations.

JULES-crop and JULES-BE
We utilized the version 5.2 of JULES model which is fully documented and publicly available at http://jules -lsm. github.io/vn5.2/. Best et al. (2011) provide the scientific description and set of equations employed to simulate the energy and water fluxes, whereas details about the carbon fluxes and natural vegetation PFTs are given in Clark et al. (2011). Furthermore, the JULES model includes many parametrizations that are in constant development to improve the representation of land surface processes. Recently, two parametrizations were made available to the scientific community aiming to explicitly represent growth and development, namely JULES-crop and JULES-BE (Littleton et al., 2020;Osborne et al., 2015). However, both parametrizations are not yet ready for fully coupled climate-surface simulations, as further evaluation across different environments and crops are still needed to refine and test its algorithms (Williams et al., 2017).
These employ different sets of equations for the allocation of plant carbon pools, leaf area index (LAI) and height as compared to the natural vegetation PFTs . In JULES-crop, crop development is simulated on each grid tile by the crop development index (DVI), which is mainly used to control biomass partitioning among plant carbon pools, specific leaf area (SLA) and senescence over time. JULES-crop also simulates carbon remobilization from reserves to harvestable parts, and employs allometric equations to determine canopy height (CANH) from stem biomass (e.g. equations A1-A7 of Williams et al., 2017). In JULES-BE, instead of DVI, a set of allometric relationships is employed to simulate carbon partitioning over time. A balanced-growth LAI is calculated to derive the step-rate CANH, which, in turn, is used to determine carbon allocation on leaves, stems and roots (e.g. equations S1-S8 of Littleton et al., 2020).
When switched-on, these parametrizations take the daily rates of water and carbon fluxes (e.g. carbon assimilation, respiration, soil water movement and ET) calculated by the main JULES engine and return key vegetation components required for the next step, such as CANH and LAI. Conversely, LAI and CANH must be prescribed as model input when none of these parametrizations is selected. Other parametrizations can also be selected to simulate LAI and CANH, for example, the dynamic vegetation or the phenology parametrizations . However, crop-related outputs in these configurations are limited, whereas carbon allocation and phenology are mainly simulated by allometric relationships determined by the general PFT parametrizations.

| Study sites
A total of 11 field experiments across Brazil were compiled in this study with detailed information about sugarcane grown under different climate and soil conditions (Table 1; Figure 1). Some sites contained multiannual sequential data (plant-cane and ratoons [re-growth]) and different water treatments (e.g. rainfed vs. irrigated), totalling 25 growing seasons. Therefore, the dataset was organized in a code system to facilitate the identification for each combination of site location, crop season (plant-cane and subsequent ratoons) and water regime (irrigated and rainfed). For example, the first plant-cane season (0) at the Piracicaba-ESALQ (EQ) site under irrigated (I) water conditions is represented as EQ-0-I, whereas the 3rd ratoon (3) at the Luis Antonio (LA) site under rainfed (R) conditions is expressed as LA-3-R. At the Capivari and Valparaiso sites where two distinct planting dates were tested, the same code system also applies but adds a suffix to distinguish between the wet (W) and dry (D) periods. A summary description for all sites is presented in Table 1, whereas the full list of codes for the 25 field trials is given in Table S1.
Measurements of stalk dry matter were available in all locations, including other carbon pools (e.g. green and dry leaves, sucrose and roots). Most of the experiments T A B L E 1 Description of field experiments used for model calibration and evaluation (A, allometry; B, biomass; ET, evapotranspiration; GPP, gross primary productivity; R, root profile; SWC, soil water). Soil texture class, climate classification (Koeppen) and cultivar are given for each site. Unique codes are assigned for each growing season in Sites where one or more seasons were used in the evaluation process. a The soil of CO site has a physical constraint that limit root growth at about 40 cm depth (Typic Fragiudults). also had CANH and LAI observations, while root profile measurements were available only at the CT site (Laclau & Laclau, 2009). The EQ, PF and LA sites also contained soil moisture and ET measurements along with biomass and allometric data in nine growing seasons. In addition, the gross primary productivity (GPP) for sugarcane was also available for three growing cycles at the PF site. The data were collected for widely used sugarcane cultivars in Brazil, including the most planted cultivar RB867515. At all sites, the crop was grown under controlled conditions to avoid any effect of pests and diseases, nutrient deficiency or suboptimal management. Full details of the instruments and the methods used for collecting these measurements are shown in each citation shown in Table S1.
To evaluate the model at larger spatiotemporal scales, five microregions were selected where sugarcane is traditionally cultivated in Brazil ( Figure 1): two in the Southeast (Piracicaba and Presidente Prudente), two in the Northeast (Petrolina and Recife) and one in the Centre-West (Jataí) region. The selection of these locations was based on the availability of yield records in the SIDRA database from the Brazilian Institute of Geography and Statistics for the years 1980-2010 (IBGE, 2019), as well as aimed to cover contrasting edaphoclimatic conditions for testing the models. Any municipality with less than 80% of the timeseries  or within an area of high incidence of irrigation were not considered in our analysis (ANA, 2019; Dias & Sentelhas, 2018). Data with unrealistic repeated yield records for three or more consecutive years were also removed due to the inappropriateness in representing the interannual variability of yields.

| Model inputs and configuration
All simulations were forced with sub-daily meteorological data of downward shortwave radiation, downward longwave radiation, rainfall, air temperature, wind speed, air pressure, specific humidity, diffuse radiation and irrigation records. The weather data were obtained from meteorological stations located at each experiment and those from the National Institute of Meteorology for each corresponding microregion for the years 1980-2010 (Vianna et al., 2020). Any missing information was gap-filled using the WFDEI meteorological forcing dataset (Weedon et al., 2014).
For the simulations under future climate, we utilized the two climate scenarios SSP1 (RCP2.6, low-CO 2 ) and SSP5 (RCP8.5, high-CO 2 ) from the CMIP6 database for each microregion until the end of the century (Eyring et al., 2016). We selected five GCM projections based on output availability from the Inter-Sectoral Impact Model F I G U R E 1 Field experiments described in Table 1 (circles) and the five sugarcane producing microregions selected in this study (triangles). The light-green area represents the counties where sugarcane production is above 100,000 t year −1 (IBGE, 2019).
Intercomparison Project: UKESM1-0-LL, MRI-ESM2-0, MPI-ESM1-2-HR, IPSL-CM6A-LR and GFDL-ESM4. The meteorological variables were bias-corrected to the 1980-2010 weather station baseline using the quantile mapping bias-adjustment methodology of Lange (2019). As the CMIP6 scenarios were not on a sub-daily time scale, we assumed the diurnal time course of events is preserved in our projections.
Measured soil water retention curves were utilized to derive the Brooks and Corey (1964) model parameters (e.g. Figure S1) for the corresponding sites. When the soil water retention curves were missing, we employed pedotransfer functions adjusted for Brazilian soils to find the appropriate soil water retention parameters (Marthews et al., 2014). Similarly, soil density, organic carbon and textural components were used to determine the heat capacity and conductance. A maximum soil depth of 10.8 m with the vertical discretization of 14 layers was assumed in all simulations, following the proposed improvements by Harper et al. (2021). The soil information for each of the five microregions were sourced from previous studies using the predominant soil type classification (Dias & Sentelhas, 2018;Vianna & Sentelhas, 2016).
The planting and ratooning dates were set accordingly to the given information for each field experiment (Table 1). To represent the main growing conditions in Brazil, our 30-year simulations and future projections considered a 12-month growing cycle with the ratooning date set at the mid of the milling period for each microregion: February for Recife and June for the other regions (Zheng et al., 2022). A spin-up time of 180 days before the planting/ratooning was assumed in each location to minimize the effect of the initial soil moisture content on the simulation results as described in Figure S2.

| Model calibration
The model was calibrated with field and literature data following three main steps. First, the PFT parameters representing C 4 grass species in JULES were refined to reproduce the response of sugarcane to carbon dioxide, temperature, and low-water availability (CTW-response). At this step, JULES-crop and JULES-BE parametrizations were switched-off, and the CANH and LAI data obtained from the field experiments were interpolated to continuous daily values and prescribed as model input at this step ( Figure S3). Therefore, the only processes considered at this step were those related to carbon fixation, plant respiration and stomatal conductance which are simulated by the JULES engine . This allowed us to calibrate the water and carbon fluxes without bias or compensation associated with the JULES-crop and JULES-BE parametrizations. JULES can simulate the effect of limited nitrogen (N) in natural vegetation, but we have not considered the N-response in this study because (i) JULES-crop is not yet adapted to capture the effect of N throughout crop growth (Williams et al., 2017) and (ii) no observational data to support such calibration was available in our dataset.
After the CTW-response calibration, the second and third steps focused on fitting the parameters of JULEScrop and JULES-BE parametrizations, respectively. We followed the same methodology employed by Williams et al. (2017) for JULES-crop and Littleton et al. (2020) for JULES-BE, aiming to simulate the phenology and allometry patterns observed in the sugarcane fields. To ensure JULES-crop and JULES-BE parametrizations applied the same CTW-response adjusted in the first calibration step, the trait-based physiology feature of JULES was switchedoff in all simulations (l_trait_phys = F). The same multilayered canopy radiation scheme (can_rad_mod = 6) was selected for both parametrizations, which was extensively tested by Williams et al. (2017).
We used data from five growing seasons to calibrate the models: first to third ratoon seasons at EQ site, second ratoon of PF site and the third ratoon season of LA site (Table 1). These sites were selected because of the large amount of observational data required to derive the allometric relationships for sugarcane, prescribe LAI and CANH in step one and to make it possible to compare our performance results with previous sugarcane modelling studies (Marin et al., 2015;Vianna et al., 2020).
While the calibration of JULES-crop and JULES-BE was largely supported by field-level observations, the CTW-response of sugarcane was adjusted mainly with previously reported carbon assimilation and transpiration response to CO 2 and air temperature (Peixoto & Sage, 2017;Stokes et al., 2016), and maintenance respiration (Liu & Bull, 2001). We also used the response functions of recently tested sugarcane dynamic models, namely DSSAT-CANEGRO (Jones & Singels, 2018), APSIM-sugar (Dias et al., 2019) and SAMUCA (Vianna et al., 2020) to adjust and evaluate our simulations. When a parameter could not be derived directly from the literature or field measurement, the Nelder-Mead optimization method was applied to obtain the optimum values that minimize the absolute error in the calibration set (Nelder & Mead, 1965).

| At controlled field conditions
All the remaining experimental data not used in calibration were employed to evaluate the parametrizations at field scale by means of the statistical indexes performance: coefficient of determination (r 2 ), accuracy (d), root mean squared error (RMSE) and the Nash-Sutcliffe model efficiency (NSE; Wallach et al., 2018). These indexes were calculated after each calibration step to verify the isolated impact of JULES-crop or JULES-BE parametrizations on the water and carbon fluxes simulations. The output variables selected to evaluate model statistical performance were stalk dry mass (SDM), LAI, CANH, soil water content (SWC), ET and GPP. We also compared the fitted allometric relationships of JULES-crop and JULES-BE with observations at field trials and with other sugarcane models to assess the calibration consistency in representing sugarcane growth.

| At regional scale and under climate change conditions
This analysis aimed to assess the performance of the new calibrations in representing the spatial patterns of average yields and its interannual variability in five contrasting producing regions. It is widely known that process-based crop models are not able to consider all the driving factors affecting crop yields in real conditions. This concept is defined as the yield gap, and was also determined for sugarcane in Brazil using the IBGE database and a number of crop models (Dias & Sentelhas, 2018;Marin et al., 2016). Therefore, to compare our simulated SDM with IBGE yield records, we multiplied our SDM simulations by the water-limited yield gap factors (Ya/Yw) reported for each corresponding microregion in Dias and Sentelhas (2018). Both simulated and observed yields were expressed and compared in a dry basis assuming a constant moisture content of 75% and a harvest index of 80% (Marin et al., 2016).
To evaluate the model's ability in capturing the interannual variability of yields, we employed the same method described by Müller et al. (2017) to detrend our time series and obtain the correlation indexes between simulated and observed SDM time series for each microregion. To assess the uncertainty and consistency of our simulations, we also provide a local sensitivity analysis on the calibrated model parameters, and the effects of different ratooning dates and climate inputs provided by the five GCMs historical data in Supporting Information. Finally, the ensemble results of SDM simulations considering both models and climate projections were assessed using yield trends, and the main climatic drivers of yield changes were isolated and quantified.

| Modelling the CTW-response of sugarcane with JULES
The default C 4 species parameters values of JULES represented well the characteristic response curve observed for sugarcane photosynthesis and stomatal conductance to atmospheric CO 2 concentrations ( Figure 2). Compared to the current atmospheric CO 2 concentration (~410 ppm), the effect of atmospheric CO 2 enrichment on photosynthesis is small for C 4 species. This process, which is largely controlled by the PEP-carboxylases enzyme activity, is empirically represented in DSSAT-CANEGRO and SAMUCA models assuming that the CO 2 effect is F I G U R E 2 Sugarcane photosynthesis (a) and transpiration (b) response to atmospheric CO 2 concentration as simulated by the C 4 species parameters set of Joint UK Land Environment Simulator (JULES) model (solid black line in a, and black boxplots in b) in comparison with the response functions of SAMUCA and DSSAT-CANEGRO models (a) and measurements taken by Stokes et al. (2016) in a glasshouse experiment (red boxplot) (b). neglected above 270 ppm (Figure 2a). The same set of parameters also showed a consistent decline in crop transpiration at increasing CO 2 concentrations, with a rate of −0.07% ppm CO 2 −1 between 400 and 1080 ppm at the EQ site conditions. This mimicked the transpiration reduction of about 28% (±8%) observed in sugarcane plants cultivated under 720 ppm of CO 2 in a glasshouse experiment ( Figure 2b). Based on these results, the parameters related to the plant response to CO 2 concentration were not altered in this study. The response curve of the maximum rate of carboxylation of Rubisco (V cmax ) to temperature observed for sugarcane is characterized by a steep increase at suboptimal temperatures above 10°C, and a rapid drop in a relatively narrow thermal phase between 35°C and 47°C ( Figure 3a). We refined the five parameters (tupp_io, tlow_io, q10_leaf, neff_io and nl0_io, Table S2) that control the response of V cmax to temperature in JULES using field observations ( Figure 3a). The shape of the curve obtained for sugarcane was asymmetric with a peak at around 35°C, with lower V cmax below 10°C and above 42°C. This agreed with the overall pattern adjusted in trapezoidal equations used in other sugarcane models to scale photosynthesis; nevertheless, our calibration showed a more realistic response when compared to V cmax data (Figure 3a).
Maintenance respiration of sugarcane is also temperature dependent as demonstrated in the literature (Figure 3b). Using nitrogen content in sugarcane roots, stems and leaves from the literature (Table S2), we were able to fit the response curve of maintenance respiration to temperature in JULES. In the model, leaf dark respiration is calculated at the canopy layer level by multiplying V cmax by a fixed dark respiration parameter (fd_io , Table S2). Plant maintenance respiration is then calculated by relating the whole canopy dark respiration with the current carbon pools at roots, stem and leaves, and the mass ratio parameters of nitrogen-to-carbon in roots (nr_nl_io) and stems (ns_nl_io) as related to leaves (Table S2).
Because of the dependency between V cmax and respiration response curves employed in JULES, our adjustment was best fitted to the observations within the 10 to 30°C temperature range. This also simulates a theoretical drop in respiration rates above 40°C, similar to the DSSAT-CANEGRO model, while there is no observational evidence of maintenance respiration decline at this temperature threshold for sugarcane (Figure 3b). The growth respiration in JULES is set to a fixed proportion of assimilated carbon (r_grow_io), for which we set the value of 0.242 gC gC −1 that is employed in sugarcane modelling studies.
Drought stress is common in sugarcane fields; however, observational data capturing the direct effect of low water availability in sugarcane photosynthesis and biomass are scarce. Two of the field experiments included in this study, CT and UN (Table 1), aimed to capture the effect of drought stress on sugarcane growth but did not show any significant differences between irrigated and rainfed treatments. At UN, the final SDM were 29.1 and 31.6 t ha −1 in the rainfed and irrigated treatments, respectively, whereas minor differences were also found between irrigated and F I G U R E 3 Response of maximum rate of carboxylation (V cmax ) (a) and respiration (b) to temperature. Black lines represent the adjustment made with Joint UK Land Environment Simulator (JULES) using measured data (filled markers) for sugarcane (Liu & Bull, 2001;Peixoto & Sage, 2017;Sage et al., 2014). The grey dashed line in (a) is the maize crop calibration of JULES-crop provided by Williams et al. (2017). Coloured lines in (a and b) are, respectively, the trapezoidal function and the respiration rates used by the DSSAT-CANEGRO, APSIM-sugar and SAMUCA models. Solid and dashed lines in (b) represent, respectively, the crop respiration when roots and leaves are the dominant carbon pool (10% stems) and when the stem is the major pool (75% stems). APSIM-sugar does not simulate maintenance respiration; therefore, it is not included in (b). rainfed SDM, LAI and CANH in the CT trials ( Figure S14). Sugarcane models well-adjusted to rainfed conditions generally assume that the photosynthetic capacity of sugarcane starts decreasing when the relative available soil water (RASW) is between 60% and 40% (Figure 4a). In JULES, the carbon assimilation of vegetation is linearly reduced with soil moisture when the RASW drops below a threshold, determined by the parameter fsmc_p0_io. Therefore, we adjusted the fsmc_p0_io = 0.6 to reproduce a drought stress response of recently assessed sugarcane models in Brazil (Figure 4a).
The root-intersect densities measured at the CT site were employed to fit the relative root density profile function of JULES for sugarcane ( Figure 4b). As JULES does not simulate the effect of soil moisture in root architecture, we made no distinction of root parameters under rainfed and irrigated conditions. In addition, to simulate the shallower root system of 40 cm observed at the CO site, the parameter rootd_ft_io was set to 0.05 only for this site (dashed line in Figure 4b).
The performance of JULES with refined CTW-response for sugarcane showed high precision (r 2 = 0.94) and accuracy (d = 0.98) for soil moisture simulations (Figure 5a), with an RMSE = 0.027 cm 3 cm −3 and NSE = 0.93 (Table 2). ET simulations showed a wider difference compared to observations, with a precision (r 2 ) of 0.36, accuracy (d) of 0.74 and an RMSE = 1.24 mm day −1 . We found an average overprediction of 0.48 mm day −1 in ET simulations (Table 2), which is evidenced mainly at the LA site ( Figure S5d) and can be associated with uncertainties in root water uptake simulations in high sand content conditions (e.g. 74% of sand). Still, the seasonal pattern of ET simulations showed a good agreement across sites ( Figure S5). Simulations of GPP showed substantially better agreement than ET, with a precision (r 2 ) of 0.77 and accuracy (d) of 0.92. These results show a better performance when compared to simulations using the PFT parameters values for maize, in particular for GPP which showed an NSE reduction from 0.66 to −0.14 ( Figure S7; Table S3).

| Modelling sugarcane growth and development with JULES-crop and JULES-BE
After refining the CTW-response of sugarcane in JULES, we adjusted the crop-specific parameters of JULES-crop and JULES-BE to represent sugarcane development, carbon partitioning and allometry. Full details about this process are provided in Supporting Information along with the parameter values and rationale for adjustment. After calibration, both JULES-crop and JULES-BE parametrizations replicated the growth patterns of SDM, LAI and CANH for sugarcane at the field trials used for validation ( Figure 6). Precision (r 2 ) and accuracy (d) for SDM simulations were between 0.72 and 0.92, with an average RMSE of 6.4 t ha −1 for both models (Table 3). Simulations of CANH were similar for both JULES-crop and JULES-BE, which can be partly attributed to the allometric relationships adopted to simulate this same crop variable. JULES-BE also presented higher precision for LAI simulations than JULES-crop, but both parametrizations showed similar accuracy (d ~0.75). The performance of SWC, ET and GPP, was negatively affected when compared to the simulations where LAI and CANH were prescribed (Tables 2 and 3). Yet, SWC simulations presented high-performance indexes (r 2 ~0.88, d = 0.96) for both parametrizations, whereas the performance of ET simulations was more negatively affected for JULES-BE than JULES-crop. This is also evident for the GPP simulations, where JULES-crop showed superior performance indexes (r 2 = 0.78, d = 0.92) than JULES-BE (r 2 = 0.57, d = 0.81).
This can be mostly associated with the differences in the representation of LAI over time between parametrizations (Figure 6b,e). A peak followed by a steep decline in LAI is verified at around 200 days after planting for JULEScrop, while JULES-BE showed a substantially narrow LAI variability across sites and a systematic overprediction at earlier developmental stages. These can be explained by the approximations made to describe the biomass partitioning and SLA for JULES-crop, and because the allometric relationships employed by JULES-BE do not capture the whole variability observed in our dataset ( Figure S13). Simulations of SDM growth by JULES-crop show a delay when compared to observations at PF and CT sites, with a steep growth until the harvest date, which can also be attributed to the approximation of carbon partitioning functions of JULES-crop.
None of the field experiments with rainfed versus irrigated treatments showed differences in sugarcane growth (CT and UN), but the rainfed field of CO where the root system was limited to 40 cm depth presented lower SDM, LAI and CANH than other sites, as a result of moderate water stress. Both parametrizations were able to simulate lower SDM, LAI, CANH for the CO site, whereas no or a small variation was found for simulations where rainfed versus irrigated trials did not show a significant difference ( Figure S14). This may indicate both parametrizations' ability in capturing water stress in sugarcane simulations.
When compared to the performance of other sugarcane dynamic models evaluated in some of our sites ( Figure S16), the LAI simulations of JULES-crop and JULES-BE simulations showed slightly higher correlations (r = 0.69 and r = 0.77, respectively), but also higher RMSE. JULES-crop had a higher RMSE for LAI than other models (1.69 m 2 m −2 ), whereas the RMSE of JULES-BE was similar to the APSIM-Sugar model (1.35 m 2 m −2 ), but still higher than DSSAT-CANEGRO and SAMUCA.

| Spatiotemporal variability of sugarcane yields simulated by JULES-crop and JULES-BE
Simulations of interannual variability of SDM for sugarcane were similar for JULES-crop and JULES-BE for most of the years (Figure 7). JULES-crop showed a stronger signal than JULES-BE, as observed in some years (e.g. years 1997-2001 at Recife and Petrolina) and evidenced by the higher coefficients of variation (CV; Table 4). Although the correlation index (r) between simulations and detrended yields reported by IBGE were not high (0.36-0.56), we found significant correlation (p < 0.05) for at least one parametrization in all locations (Figure 7). JULES-crop showed the higher correlation as compared to JULES-BE in most locations, except for Presidente Prudente and Recife.

F I G U R E 5 Comparison between Joint UK Land Environment Simulator (JULES) simulations with non-explicitly crop growth and observations of soil water content (SWC) (a), evapotranspiration (ET) (b) and gross primary productivity (GPP) (c) at the evaluation trials.
The blue line represents the linear correlation between simulated and observed data, and the dashed red line is the 1:1 reference line. Site codes are given in Table 1 and Table S1.
T A B L E 2 Performance of Joint UK Land Environment Simulator simulations with non-explicitly crop growth for soil water content (SWC), evapotranspiration (ET) and gross primary production (GPP) at the evaluation trials (

F I G U R E 6
Stalk dry mass (SDM), leaf area index (LAI) and canopy height observed (markers) and simulated (lines) by Joint UK Land Environment Simulator (JULES)-crop (a-c) and JULES-BE (d-f) in the validation field trials across Brazil. Site codes are given in Table 1 and Table S1.  Petrolina is located in the Brazilian semiarid region, where crops experience extreme drought conditions leading to a high interannual yield variability for rainfed sugarcane (CV > 20%, Table 4). In Recife, the predictions of both parametrizations corresponded well to yield oscillations, notably between 1996 and 2008. Some events F I G U R E 7 Interannual variability of sugarcane stalk dry matter (t ha −1 ) simulated by Joint UK Land Environment Simulator (JULES)crop and JULES-BE (solid lines), compared with the interannual variability reported by IBGE from 1980 to 2010 (circles). The correlation between simulated and observed de-trended yields is given for JULES-crop and JULES-BE in r(C) and r(B), respectively, and asterisks indicate the p-values of significance (**p < 0.01 and *p < 0.05).

T A B L E 4 Average and coefficient of variation (CV) of stalk dry mass yields simulated for sugarcane by Joint UK Land Environment
Simulator (JULES)-crop (crop) and JULES-BE as compared to the reported yields by IBGE before (Yw) and after (Ya) the corrections to account for the yield gap (Ya/Yw) in each region were not captured well in our simulations, such as the year 1993 in Petrolina and Recife, 1987in Piracicaba and between 1998and 2001 in Jatai for JULES-BE. Yet, simulations of both parametrizations showed good agreement with IBGE data for Presidente Prudente, even with low interannual variability. The average SDM simulated by JULES-crop was nearly 10% higher than JULES-BE, which is explained by the differences in SDM growing curves at the end of the season (Figure 6). Sugarcane yields reported by IBGE were on average 48% lower than those simulated by the JULES-crop and JULES-BE models (Table 4). This difference was lower for Petrolina and higher for Recife, which are mainly explained by other limiting factors that the models are not yet ready to simulate, such as pest and diseases, nutrients deficiency or suboptimal managements. After correcting our simulated yields with yield gap factors (Ya/Yw), the RMSE of simulations dropped from 20.1 to 5.8 t ha −1 , improving the representation of yields reported by IBGE (Table 4 and Figure S18a). Our simulations also replicated the pattern of yield variability observed in different regions, with higher-to-lower CV% from the semi-arid region of Petrolina to favourable edaphoclimatic conditions in Southeast ( Figure S18b).

| Climate change impact on sugarcane yields simulated by JULES-crop and JULES-BE
Future yield projections for sugarcane using the calibrated JULES-crop and JULES-BE showed an overall low sensitivity to climate change (Figure 8). In general, all locations indicated a small positive trend for SDM yields for both scenarios (SSP1 and SSP5) until the end of the century (Table 5), except Jataí which showed a clear negative trend for the SSP5 scenario between 2070 and 2100 (−0.18 t ha −1 year −1 ). This was also verified in the simulations at Presidente Prudente for the same period and scenario, but at lower rate (−0.05 t ha −1 year −1 ). In contrast, simulations in Petrolina for the SSP5 scenario are pointing to the higher trends among locations and scenarios (+0.17 t ha −1 year −1 ). However, sugarcane yields in this region by 2100 would be still below the current yields of regions in Centre-South Brazil and with high interannual variability. Although yield projections of JULEScrop show larger variability as compared to JULES-BE, the same trend patterns were found for both models ( Figure S26).
The temperature was the main driver for negative trends observed in Jataí and Presidente Prudente for scenario SSP5 (Figure 9). Projected daytime temperatures above 35°C are expected to decrease the rate of carbon assimilation, as it falls on the rapid drop phase of the T-response curve for sugarcane (Figure 3). In Presidente Prudente, the effect of high temperature is only noticed in the end century (2080-2100), whereas in Jatai, the average of the temperature factor starts decreasing much earlier (2066). Simulations for Petrolina in scenario SSP5 also showed a negative impact of high temperature from 2072 (Figure 9), but the small positive trend projected for precipitation (+1.13 mm year −1 , p < 0.001) was enough to compensate the small positive trend observed in SDM simulations ( Figure 8; Table 5). This effect was likely amplified by the higher atmospheric CO 2 concentration in the SSP5 scenario (Figure 2; Figure S27), as it induces a better water use efficiency in sugarcane. Despite the rising temperatures in the SSP5 scenario, we did not find any trends in droughts and temperature that significantly affected the SDM simulations in Piracicaba and Recife regions.

| Adjusting the sugarcane response to environmental conditions with JULES
Refining the CTW-response has become important for improved agricultural representation in integrated assessment models (Ruane et al., 2017). While field experiments datasets were key to support calibration of crop phenology and biomass partitioning, pre-existing sugarcane model parametrizations and literature were important to adjust and evaluate JULES for conditions not fully covered in our datasets (e.g. CO 2 enrichment, Figure 2 and droughtstress, Figure 4). We verified that the default parameter values used in the photosynthesis-stomatal conductance method implemented in JULES  replicated well the response of sugarcane to CO 2 concentrations (Figure 2), thereby not requiring any change.
We used V cmax observations from two independent studies to refine the carbon assimilation response of sugarcane to temperature (Peixoto & Sage, 2017;Sage et al., 2014), where the shape of response function replicated the trapezoidal pattern simulated in other welltested dynamic models (Figure 3a). Sugarcane leaves and roots have a relatively higher respiration rate than stems (Liu & Bull, 2001); thus, the respiration rate per unit of crop biomass is expected to be higher at earlier developmental stages than at the maturation stage, when stem biomass is the dominant carbon pool. DSSAT-CANEGRO and SAMUCA models differentiate crop respiration rates by carbon pool, and JULES was also able to mimic this. However, the maintenance respiration response at high temperatures is a source of uncertainty, particularly because of the divergence verified among models' response curves and the lack of observations at high temperatures (Figure 3b).
Measurements of the direct effect of drought on sugarcane biomass and photosynthesis are scarce. In a field trial limiting the root of depth to 0.97 m, Smit and Singels (2006) found a significant increase in leaf senescence and lower developmental rates of sugarcane when the RASW dropped from 70% to 30%. This agrees with the other sugarcane models that we used for calibration of the drought response function in JULES (Figure 4a). Furthermore, both JULES-crop and JULES-BE were able to better capture the moderate drought effect observed at the CO site when the root system was limited to 40 cm depth ( Figure S15). This suggests that unless root penetration is not limiting, sugarcane can resist moderate dry spells.
Field observations from Laclau and Laclau (2009) also show that root plasticity may offset drought stress by better exploring the soil profile as compared to irrigated conditions. They also found a maximum root depth greater than 4 m in both irrigated and rainfed treatments (irrigated = 4.25 m, rainfed = 4.70 m), with unnoticeable difference on the root front velocity between both treatments. Although JULES-crop can simulate variable root depth throughout the growing cycle, this feature is switched off in previous studies by setting the parameter rt_dir_io to zero (Osborne et al., 2015;Williams et al., 2017). We also F I G U R E 8 Historical and future projections of stalk dry matter yield of sugarcane simulated in five traditional producing regions with Joint UK Land Environment Simulator (JULES)-crop and JULES-BE, considering five general circulation models (GCMs; UKESM1-0-LL, MRI-ESM2-0, MPI-ESM1-2-HR, IPSL-CM6A-LR and GFDL-ESM4) and two scenarios (SSP1 and SSP5). Solid lines are the median and the ribbons are the standard deviations of the five GCMs and both models' ensemble. The blue lines represent the average between JULES-crop and JULES-BE simulations driven by weather observations (e.g. Figure 7).

Crop BE Ensemble
Jatai ( (Figures 3 and 4) and the corresponding projections of daytime temperature and precipitation in five traditional producing regions, considering five general circulation models (GCMs; UKESM1-0-LL, MRI-ESM2-0, MPI-ESM1-2-HR, IPSL-CM6A-LR and GFDL-ESM4) and two scenarios (SSP1 and SSP5). Solid lines are the median and the ribbons are the standard deviations of the five GCMs and both models' ensemble. The dashed grey line in daytime T represents the 35°C threshold.
followed this because JULES-crop relates root depth as a function of carbon allocated to root biomass (rootc), which may lead to oversensitivity of drought stress particularly at early developmental stages. Other models simulate the deepening of the root front as a function of degree days only (Jones & Singels, 2018;Marin et al., 2015;Vianna et al., 2020), resulting in more realistic results. This should be considered an improvement for future versions of JULES for plant cane systems, while assuming a fixed root depth may best represent the root growth pattern in ratooning sugarcane systems (Smith et al., 2005).

| Modelling sugarcane growth and phenology with JULES-crop and JULES-BE
Both parametrizations were able to replicate the growing patterns observed at sugarcane field experiments ( Figure 6). Although requiring a higher number of parameters, JULES-crop provides more flexibility and details in sugarcane simulations as compared to JULES-BE. Some structural improvements required to refine sugarcane simulations with JULES-crop are related to the implementation of multiannual growing seasons and in carbon partitioning equations. We found high sensitivity to the carbon partitioning parameters ( Figure S19), which was also confirmed for maize (Prudente Junior et al., 2022). Cuadra et al. (2012) implemented the carbon partitioning equations of the DSSAT-CANEGRO model to simulate sugarcane growth in the Agro-IBIS model. Black et al. (2012) employed the same approach to simulate sugarcane with JULES in Ghana; however, their parametrization was not available to compare with our calibrations. Making the carbon partitioning equations of JULES-crop more flexible may be sufficient to refine simulations without the need of a dedicated parametrization for sugarcane. Sugarcane dynamic models also explicitly simulate other processes which are not considered in JULES-crop or JULES-BE. An example is the tillering process (tiller m −2 ), which controls plant population and then used to upscale stalk biomass and LAI (Jones & Singels, 2018). In both parametrizations, LAI is obtained by multiplying the active leaf carbon pools to the SLA, which is simulated with crop development in JULES-crop and considered a fixed value in JULES-BE. Four SLA observations from different sites at the early stages of sugarcane development (DVI < 0.2) are below 5 m 2 kg −1 , suggesting that SLA can be influenced by the tillering process that drives plant population at early canopy formation. After canopy closure, SLA in sugarcane is almost constant, which was replicated by JULES-crop ( Figure S12a) or represented by the sigl_io parameter of JULES-BE.
A per-site comparison reveals that the level of performance is similar between the calibration and evaluation sites. Although our simulations performed best for some specific variables and sites, we did not identify any particular site where the parametrizations were best fitted to all variables together (CANH, LAI and SDM) to characterize overfitting ( Figure S17; Figure 6). The adjusted allometric relationships also do not show evidence of overfitting to a single condition, as it clearly represents the overall pattern across the field trials ( Figures S12 and S13). This was also verified in the representation of carbon partitioning and leaf senescence for JULES-crop ( Figures S9 and S10). However, only two sites had sufficient data to adjust these processes (EQ and PF [4 trials]), requiring more information to investigate any possible overfitting for these processes in future studies.
Both parametrizations, but mainly JULES-BE, rely on allometric relationships to derive canopy information. We found that these allometric relationships hold for explaining the overall growth pattern of sugarcane growth; however, a large deviation is found across the field trials used in this study (Figures S12 and S13). This may suggest these relationships can be refined with genotype-specific data or improved to account for the response of plant expansion to drought stress and temperature (Jones & Singels, 2018;Vianna et al., 2020). This study aim was to provide a general parametrization for sugarcane in JULES and thereby we did not generate genotype-specific sets of parameters, which would require a substantially higher amount of data across genotype x environment interactions.
Crop-specific models are generally expected to have better performance than crop-generic models such as JULES-crop and JULES-BE, as they account for specific processes such as tillering, sucrose accumulation and source-sink dynamics (Marin et al., 2015). This can improve not only the overall performance of the models, but also the ability to represent other crop components that are relevant for growers (e.g. sucrose content). On the other hand, those models generally represent the crops physiology in a simplistic way, such as relaying on Radiation Use Efficiency and crop ET coefficients. Although JULEScrop and JULES-BE were not developed to simulate crop growth in such detail, still, our results showed that both parametrizations have comparable performance indexes to sugarcane models ( Figure S16) and were able to simulate stalk dry yields with good precision (r 2 > 0.70) and accuracy (d > 0.90; Table 3).
Despite showing better statistical performance in LAI simulations than JULES-crop (Table 3), a systematic overprediction of LAI at the early developmental stage was noticed in JULES-BE simulations (Figure 6e). This was also verified in the Littleton et al. (2020) results for Miscanthus and had a substantial impact on our ET and GPP simulations using JULES-BE in comparison with JULES-crop ( Figures S5 and S6; Table 3). However, the performance of simulating sugarcane stalk dry yields and canopy formation (LAI) had good agreement with observations while preserving the overall pattern of fluxes of soil moisture, ET and GPP.

| Assessing model ability in
simulating the spatiotemporal variability of sugarcane yields and the effect of climate change The highest yield simulations were found for Recife, where solar radiation and rainfall are high enough to sustain vigorous vegetative growth. In contrast, the Petrolina region located at similar latitude but in the Brazilian semi-arid region presented the lowest yields which are mainly explained by severe drought conditions. The other regions presented small differences among average yield values, being Jatai, Presidente Prudente and Piracicaba ranked from the high to low yields. These patterns agree with modelling studies of sugarcane under water-limited conditions at the same locations (Dias & Sentelhas, 2018;Marin et al., 2016;Vianna & Sentelhas, 2016). After considering the corresponding yield gaps factors (Ya/Yw), both parametrizations were able to replicate the average yield patterns reported by IBGE across the studied regions (Table 4; Figure S18). Pest and diseases, weed pressure, soil degradation due to long-term monoculture and inappropriate traffic of heavy machinery are generally associated with a decreased longevity of sugarcane ratooning system and yield gaps (Dias & Sentelhas, 2018;Marin et al., 2019).
The interannual yield simulations consistently captured the patterns of yield variability across regions and some strong oscillations observed in the time series. We found higher interannual yield variability for Petrolina, with lowest dips in the years 1993, 1997 and 1998. These periods coincide with two strong droughts that affected over 80% of the Brazilian Northeast region in the 1990s (Cunha et al., 2018). This also overlaps with the lower yields simulated in 1998 for Recife (Northeast), but was not captured in 1993 possibly due to a mismatch of regional planting/ratooning date adopted by farmers in that year ( Figure S24). The other regions showed lower interannual yield variability, except for the beginning of the 1980s which was also marked by a period of strong droughts in Southeast Brazil (Gozzo et al., 2019).
Yet, we did not find high correlation indexes (r) for our interannual yield variability simulations, ranging from 0.33 to 0.56 (Figure 7). However, our results agree to other studies that compared simulations of interannual yield variability for sugarcane with the IBGE time series. Cuadra et al. (2012) obtained a correlation index of 0.44 for the Agro-IBIS model in mesoregions located in Southeast Brazil, while Black et al. (2012) found an average correlation index of 0.51 in the same regions. In addition, Müller et al. (2017) showed that the interannual yield variability simulated at subnational level usually do not show high correlation for process-based models, which can be largely attributed to the fact that models do not capture all the driving factors for crop yields, but also to the assumptions made to approximate the real cropping system conditions.
In our simulations, we considered the soil, weather and ratooning dates to represent the prevalent conditions for each region following recent modelling studies for sugarcane (Dias & Sentelhas, 2018;Vianna et al., 2020). We also investigated how much the model parameters, ratooning dates and climate input can affect our results (Supporting Information); however, a global sensitivity analysis may still be required to fully understand the interaction between these elements and the realized model responses under plausible ranges of parameters and conditions (Pereira et al., 2021). While the introduction of yield gap factors helped explaining the spatial differences between simulated and observed yields, the effects of biotic factors or suboptimal managements on the interannual variability of yields still remain difficult to be determined. In this context, new techniques for the reconstruction of past events in cropping systems are increasingly important not only to identify yield anomalies, but also to reduce the uncertainties about land-use and management conditions for crop growth simulations (Zheng et al., 2022).
Climate change impact studies in sugarcane generally find a positive or neutral trend on stalks yields (Dias & Inman-Bamber, 2020). Nevertheless, we found that if temperatures above 35°C become more frequent in the future, a decline in net carbon uptake is expected to affect sugarcane yields if the crop does not adapt (Figures 3a and 9). This agrees with the recent findings of Flack-Prain et al. (2021), where temperature is the main driver for sugarcane yield decline in future projections for São Paulo, Brazil. In their study, the Soil-Plant-Atmosphere model was employed which simulates photosynthesis at sub-daily timestep. JULES also operates at sub-daily timesteps, allowing to account for the effect of maximum temperatures that generally occur when solar radiation is abundant and photosynthesis rates are high.
Uncertainties of crop growth projections can be substantially reduced by improving the representation of temperature response in crop models (Wang et al., 2017). Models that use the daily average temperature (T avg ) may not be able to capture this effect unless the photosynthesis response to temperature is modified to account for the daily temperature amplitude (Marin et al., 2015). Crop respiration was also a source of uncertainty in JULES simulations, which could reduce even further the net carbon uptake of sugarcane under high temperatures (Figure 3b). While yield gaps are generally attributed to suboptimal crop management, worth noting that some abiotic effects are still not included nor systematically tested in sugarcane models, such as ozone damage, waterlogging, heat stress and lodging (Gomathi et al., 2015;Moura et al., 2018). Quantifying the effects of extreme weather in sugarcane growth remains a challenge, as the appropriate field data to model the crop response under extreme conditions are rare and difficult to be obtained.
Our results were used to understand model capabilities and identify opportunities for improvement in future versions. While input parameters and observed data are lacking for the development of bioenergy models, we provide crop-related parameters that can support modelling studies for sugarcane. The newly calibrated version of JULES can be applied to help understanding of interactions between climate and bioenergy production in the present and under climate change scenarios.