The impact of climate change and climate extremes on sugarcane production

Sugarcane production supports the livelihoods of millions of small‐scale farmers in developing countries, and the bioenergy needs of millions of consumers. Yet, future sugarcane yields remain uncertain due to differences in climate projections, and because the sensitivity of sugarcane ecophysiology to individual climate drivers (i.e. temperature, precipitation, shortwave radiation, VPD and CO2) and their interactions is largely unresolved. Here we ask: how sensitive is sugarcane yield to future climate change, including climate extremes, and what are its key climate drivers? We combine the Soil‐Plant‐Atmosphere model with detailed time‐series measurements from experimental plots in Guangxi, China, and São Paulo State, Brazil. We first calibrated and validated modelled carbon and water cycling against field flux and biometric data. Second, we simulated sugarcane growth under the historical climate (1980–2018), and six future climate projections (2015–2100). We computed the ‘yield‐effect’ of each climate driver by generating synthetic climate forcings in which the driver time series was alternated to that of the historical median. In Guangxi, median yield and yield lows (i.e. lower decile) were relatively insensitive to forecast climate change. In São Paulo, median yield and yield lows decreased under all future climates projections ( x¯ = −4% and −12% respectively). At Guangxi, where moisture stress was low, radiation was the principal driver of yield variability (yield‐effect x¯ = −1.2%). Conversely, high moisture stress at São Paulo raised yield sensitivity to temperature (yield‐effect x¯ = −7.9%). In contrast, a number of other modelling studies report a positive effect of increased temperatures on sugarcane yield. We ascribe the disparity between model predictions to the representation of key phenological processes, including the link between leaf ageing and thermal time, and the role of ageing in driving leaf senescence. We highlight climate sensitivity of phenological processes as a key focus for future research efforts.


| INTRODUCTION
Sugarcane (Saccharum spp.) is grown across 26 million ha of land in more than 100 countries throughout the tropics and subtropics (FAO, 2019). In 2017, the three biggest global producers (Brazil, India and China) collectively harvested over one billion tonnes of cane. In Brazil, production is concentrated in the southeast and northeast of the country (Adami et al., 2012), while in China, production is largely focused in the southern region (Zhao & Li, 2015). The harvested crop is used both as a source of food and renewable energy (i.e. bio-fuel). A 'cash crop', sugarcane is economically and socially important (Ferreira, 2012;Goldemberg, 2007;IBGE, 2017). Farmed area continues to expand, driven by higher market prices amid a growing demand for ethanol (Marin et al., 2011;Song et al., 2016). From 2010 to 2017, annual sugarcane yields (kg ha −1 ) in the top producing countries varied by 5%-16% (FAO, 2019). The role of climate in driving inter-annual variation in sugarcane yield is confounded by differences in management, technological advancements, soils, crop variety and geography. Climate change projections across the major sugarcane-producing regions include rising surface temperatures (IPCC, 2014). Shifts in precipitation regime are also predicted, but are highly uncertain, and differ between the major producing regions (Magrin et al., 2014). Quantifying the climate sensitivity of sugarcane is therefore critical to predicting future crop yields and managing for resilient production.
Sensitivity analyses have offered important insights into the response of sugarcane yield to individual climate drivers (i.e. temperature, precipitation, CO 2 , shortwave radiation and vapour pressure deficit (VPD; Gouvêa et al., 2009;Marin et al., 2013). These analyses used crop models to quantify how yield varies with changes in temperature, precipitation and CO 2 relative to historical observed climates (Jones & Singels, 2018;Marin et al., 2013;Webster et al., 2009). For example, applying the DSSAT-Canegro model (V4.5_C2.2) to sugarcane sites in South Africa and Brazil, Jones and Singels (2018) report increases in stalk dry mass of between 3% and 7% in response to a 3°C increase in temperature from the baseline. By separating the yield effects of individual climate drivers, these studies are able to support targeted mitigation strategies (e.g. irrigation where future yield is constrained by precipitation). However, such sensitivity analyses do not account for observed covariance between climate drivers, nor the magnitude of the climatic shift predicted locally. In addition, the impact of projected changes in shortwave radiation and humidity (VPD) on sugarcane yield has received less attention (Gouvêa et al., 2009), compared to that of temperature, precipitation and CO 2 . This is despite (i) the negative impact of increased precipitation on incident shortwave radiation (via cloud cover), and (ii) the positive correlation between VPD and surface temperature.
An alternative approach to computing the climate sensitivity of crop yields is to force models with future climate projections (Eyring et al., 2016;Ruan et al., 2018;Taylor et al., 2012). The use of dynamically downscaled global climate models has allowed local yield predictions to be generated across a range of future climate change scenarios, including the Intergovernmental Panel on Climate Change's Special Report on Emissions Scenarios (SRES), and the Representative Concentration Pathways (RCPs; Van Vuuren et al., 2011). For example, using the downscaled HadCM3 global model under the SRES A1B climate scenario, de Carvalho et al. (2015) predict a 20% and 27% reduction in sugarcane yield by 2071-2100 for two different municipalities in North-eastern Brazil, as a result of projected declines in precipitation. However, to our knowledge, modelling studies have not yet included the more recent, Shared Socioeconomic Pathways (SSPs), which capture the wide range of uncertainty around future sustainable development, regional rivalry, inequality and fossil-fuelled development (Riahi et al., 2017).
Uncertainty in sensitivity analyses is derived from both the future climate data used, and from model representation of important eco-physiological processes. Two key crop models used to simulating sugarcane growth are DSSAT-Canegro and APSIM-Sugarcane. Both models couple detailed representations of sugarcane development with pre-existing crop model platforms (DSSAT, Decision Support System for Agrotechnology Transfer; APSIM, Agricultural Production Systems sIMulator). With respect to photosynthesis, DSSAT-Canegro (V4.5_C2.2) assumes an optimal daily mean temperature threshold of between 20°C and 40°C, beyond which photosynthesis declines (Ebrahim et al., 1998;Jones & Singels, 2018). In APSIM-Sugarcane, photosynthesis is optimal between mean daily temperatures of 15°C and 35°C, and declines to zero at 5°C and 50°C (Keating et al., 1999). However, Sage et al. (2013) report that at a fixed VPD, instantaneous net photosynthesis in sugarcane (µmol m −2 s −1 ) increases with leaf temperature until ~33°C and declines thereafter. By using the daily mean temperature, both models ignore diurnal variation in temperature, and as such, do not capture the non-linear photosynthetic response to midday temperature highs. Furthermore, both models simulate photosynthesis as a function of air temperature not leaf temperature. In doing so, neither model accounts for the linkages between canopy energy balance and photosynthetic rate. Canopy energy balance is dependent on air temperature as well as stomatal cooling, and thus plant water status.
With respect to CO 2 fertilization, both DSSAT-Canegro and APSIM-Sugarcane are limited to empirical functions. In DSSAT-Canegro, a zero-fertilization response is simulated above 270 ppm and in the absence of water stress (Jones & Singels, 2018). In APSIM-Sugarcane, CO 2 fertilization effects on transpiration efficiency and radiation use efficiency are calculated as a linear function (Ruan et al., 2018;Webster et al., 2009). Neither model explicitly simulates the effect of changes in CO 2 supply on instantaneous leaf carboxylation rate. Moreover, both these approaches empirically parameterize the coupling between water demand, supply and photosynthesis, and as such do not allow for observed plant water-use optimization under soil moisture stress (Bonan et al., 2014).
To resolve current uncertainty, process-based models which couple energy, water and carbon fluxes (i.e. photosynthesis) on hourly timescales are needed. For example, the Soil-Plant-Atmosphere (SPA) model (Williams et al., 1996) links hourly/half-hourly leaf-level photosynthesis (Farquhar & Von Caemmerer, 1982) and transpiration (Penman-Monteith equation) via a stomatal conductance algorithm that optimises leaf carbon assimilation, within the limits of plant water supply. In SPA, the optimum temperature for photosynthesis is 30°C, and simulated photosynthetic capacity is dependent on temperature response curves fitted via a polynomial relationship (McMurtrie et al., 1992;Rastetter et al., 1991).
To be useful to farmers, climate change impact predictions should capture shifts in median yield, and the risk that climate extremes cause marked reductions in sugarcane yield. Average changes in sugarcane yield are widely reported (Knox et al., 2010;Marin et al., 2013). For example, using the APSIM model in Southern China, Ruan et al. (2018) predict up to a 28% increase in sugar yield in the 2030s, up to a 44% increase in the 2060s, and up to a 41% increase in the 2090s under the RCP4.5 emissions scenario (relative to the years 1961-2010). Fewer studies have reported on inter-annual variation in yield (Biggs et al., 2013;Singels et al., 2014) and the impact of extreme climate events (Mall et al., 2016;Zhao & Li, 2015). Current estimates indicate drastic reductions in sugarcane yield of more than 50% relative to the baseline mean, as a result of projected drought events in North-eastern Brazil (de Carvalho et al., 2015). However, predicted shifts in the frequency and severity of potential yield reductions relative to baseline variability remain absent, as does a broader geographical context.
We aim to reduce current uncertainty on the impact of projected climate change and climate extremes on sugarcane production. We focus on experimental sites in two major sugarcane-producing regions where data for model calibration and validation are available; São Paulo State, Brazil (2005-2007rainfed), andGuangxi, China (2016-2018;irrigated). We apply the terrestrial ecosystem model SPA (Williams et al., 1996(Williams et al., , 1998, which includes a detailed process-based representation of photosynthesis and water transport simulated on an hourly timescale. The ability of SPA to model cropland systems has previously been demonstrated across a number of European and US sites, including C 4 crops (Sus et al., , 2013Wattenbach et al., 2010). Following model calibration and validation using detailed in situ data, we force SPA with historical climate data, and future climate data generated under six different SSPs as part of the CMIP6 project (MOHC-UKESM1.0; 2015-2100). We ask: Q1a. What is the sensitivity of median sugarcane yield in both study regions to future climate projections? b. What are the key climate drivers of median sugarcane yield variability?
Q2a. How sensitive are yield lows in both study regions to future climate projections?
b. What are the major climate drivers of low yields? Q3. What determines the sensitivity of sugarcane yield to projected climate change and differences in responses between study regions?

| Site characteristics
The Chinese field site (hereafter referred to as CH), located in Chongzou, Guangxi, China (22.52°N, 107.39°E, elevation 181 m), has hosted irrigation and fertilizer experiments from 2016 to 2018. Mean annual temperature is 23°C, with a mean annual rainfall of 1119 ± 181 mm. Rainfall is highest during June-August, and reduced during November-February. The planted cultivar is Liucheng-05136, which is common across Guangxi. The plant crop was grown between April and December 2016, with subsequent ratoon crops harvested in December 2017 and 2018. In October 2016, old leaves were stripped from the crop, as is common practice in the region. Leaf stripping did not occur in other years. Experimental plots were 64 m 2 (8 × 8 m), and the distance between adjacent plots was 2 m. We selected plots which were not nutrient limited (fertilizer application >250 kg ha −1 ). Plots selected differed between years, and received between 52-72 mm irrigation in 2016 (n = 6), 160-175 mm in 2017 (n = 24) and 360-385 mm in 2018 (n = 7). Further details are given in Hu et al. (2019).
The Brazilian field site (here after referred to as BR) was located in Luiz Antonio municipality, São Paulo State, Brazil (−21.63°N, −47.78°E). The planted varieties were SP81-3250, SP83-2847 and RB86-7515, typical of Southern Brazil. Mean annual temperature was 22°C, with a mean annual precipitation of 1517 ± 274 mm . Precipitation was highest during December-February and lowest during June-August. The site covered >400 ha, at an altitude of 552 m on Typic Haplustox soils. The distance between planted rows was 1.4 m. We utilize data from the second (April 2005-May 2006) and third (May 2006-May 2007 ratoon. The site was not subject to experimental treatments. Urea-based fertilizer equivalent to 56 kg N ha −1 , 150 kg K 2 O ha −1 of potassium and 400 kg ha −1 of limestone were applied 1 week post-harvest. Further site details are presented in Cabral et al. (2012Cabral et al. ( , 2013 and Cuadra et al. (2012).

| Field data
Field measurements of leaf area index (LAI), soil moisture, leaf-level gas exchange and crop yield were taken across three seasons at the CH field site (2016)(2017)(2018). LAI measurements were taken using a LI-COR LAI-2200C (LICOR Inc.) at each plot on a weekly-fortnightly basis. Volumetric soil water content was sampled at depths of 10, 20, 30, 40, 60 and 80 cm every 3-5 days using TRIME-PICO-IPH Time Domain Reflectometry (TDR). Leaf-level gas exchange measurements were conducted between September and November 2017 using a LI-6800 Portable Photosynthesis System (LI-COR Inc.). Measurements were taken across a range of PAR (325-1985 µmol m −2 s −1 ) and temperature (26-43°C) values. Post-harvest, canes were dried at 80°C to determine stalk dry mass. The brix content was sampled and used to determine total sugar yield.
At the BR field site, measurements of soil moisture, surface gas exchange, stalk biomass and LAI were taken across two seasons (2005)(2006)(2007). Soil water content was measured using 10 reflectometers (CS615, Campbell SI) installed in a vertical profile, with each sensor representing layers of 30 cm, down to 3 m depth. The BR field site hosted an eddy flux tower, the data from which was processed to provide half-hourly estimates of surface gas exchange (Cabral et al., 2013). Aboveground biomass was sampled approximately every 20-days and on harvest days. Samples were separated into stalk, live leaves and dead leaves, prior to being dried and weighed. Specific leaf area estimates were used to derive LAI estimates from leaf biomass.

| The Soil-Plant-Atmosphere model (SPA)
The Soil-Plant-Atmosphere model (SPA) is a hydrodynamic terrestrial ecosystem model (Williams et al., 1996). In SPA, a detailed eco-physiological scheme is used to calculate leaflevel photosynthesis (Farquhar & Von Caemmerer, 1982), and leaf-level transpiration is estimated using the Penman-Monteith equation. Canopy layers are partitioned between sunlit and shaded fractions. Transmittance, reflectance and absorption of long wave, near infra-red and direct and diffuse photosynthetically active radiation (PAR) is determined for each canopy layer (determined by Beer-Lambert's Law) and the soil surface using a radiative transfer scheme.
SPA simulates water fluxes along the soil-plant-atmosphere continuum, using a stomatal conductance model that explicitly links leaf gas exchange and plant hydraulic supply (Williams et al., 2001). The representation of stomatal conductance within SPA optimises leaf carbon gain per unit nitrogen and limits stomatal opening to keep leaf water potential above a critical value. Hydraulic supply is determined by the difference between soil water potential and minimum sustainable leaf water potential, together with the hydraulic resistance of the soil-to-leaf pathway (Meinzer & Grantz, 1990). Simulated latent heat is thus a function of soil moisture and plant hydraulics. Water inputs to the soil via precipitation take into account canopy interception, evaporation and drainage. Water retention through the soil profile is determined by soil texture (Saxton et al., 1986). Water evaporation from the soil surface is calculated as a function of soil water content and soil radiation balance (Amthor et al., 1994). Aboveground, hydraulic conductance is computed assuming resistance to xylem water supply increases with the height of the canopy layer (Williams et al., 1996).
The partitioning of carbon between plant organs (root, leaves, stem, storage organs) is dependent on development stage, where development stage relates to crop physiological age and morphological characteristics (Vries, 1989). Development rate in SPA is accelerated as accumulated air temperature increases. Leaf senescence occurs due to self-shading (when LAI exceeds a critical threshold) or leaf ageing (determined by developmental stage; Laar et al., 1997). Leaf senescence rate is also sensitive to weighted soil water potential (Ferreira et al., 2017). Following senescence, 50% of leaf carbon is remobilized, while the remainder enters the dead foliage pool. A full description of the crop module within SPA is presented in Sus et al. (2010). Further modifications were made to SPA such that (i) the partitioning of stem C into fibre and sucrose is modelled as a function of soil moisture stress, and (ii) leaf stripping events are captured within the model (described in full in Data S1).

| Climate drivers
Local, hourly meteorological data were used to run SPA during model calibration and validation for the specific field seasons at each site. Air temperature, solar radiation, relative humidity, wind speed and rainfall data were retrieved from weather stations located at the CH and BR field sites. Where site-level data were absent (i.e. record gaps), estimates from the nearest weather station were retrieved, bias corrected (via linear regression), and used to fill (Longzhou, China: 22.37°N, 106.75°E; Campo Fontenelle, Brazil: −21.98°N −47.34°E; São Carlos, Brazil: −21.02°N −47.88°E). Gaps in shortwave radiation estimates were filled using hourly interpolated data from the MERRA-2 re-analysis product, following bias correction (Gelaro et al., 2017). Remaining gaps <8 h were filled by spline interpolation. Historical meteorology data for each site  were retrieved from the MERRA-2 re-analysis product and interpolated to generate an hourly time series for this period.
Forecasts of future meteorological data were retrieved from the CMIP6 (MOHC-UKESM1.0) projections; SSP119 (SSP-based RCP scenario with very-low radiative forcing by the end of the century), SSP126 (SSP-based RCP scenario with low radiative forcing by the end of the century), SSP245 (SSP-based RCP scenario with medium radiative forcing by the end of the century), SSP370 (baseline scenario with a medium to high radiative forcing by the end of century), SSP434 (mitigation scenario with low radiative forcing by the end of the century) and SSP585 (SSP-based RCP scenario with high radiative forcing by the end of century; Good et al., 2019). As atmospheric CO 2 increases through projections SSP119, SSP245, SSP370 and SSP585, mean annual temperature typically increase ( Figure 1). VPD rises at BR, and is increasingly variable for both sites. PAR declines between projections SSP119 and SSP126, and increases between SSP370 and SSP585, mirroring precipitation trends. The daily estimates of minimum and maximum temperature, VPD, shortwave radiation, precipitation and wind speed were used to generate hourly time series (see Data S1). All projections were bias corrected against the MERRA-2 time series for the years 2015-2018. The effect of bias correction on future climate projections is presented in Data S1.

| Model calibration and validation
SPA was parameterized using field and literature-based estimates of vegetation parameters. For the BR field site, timeaveraged leaf N content (1.0 g N m −2 ) was calculated from field measurements. For the CH field site, leaf N content (1.1 g N m −2 ) was derived from time-averaged leaf N (%) measurements. Leaf mass per unit area estimates were site specific (95 ± 52 g m −2 BR) or based on literature estimates (70 g m −2 CH; De Vries et al., 1974). Minimum leaf water potential was set at −2.5 MPa reflecting field reports (Inman-Bamber & De Jager, 1986;Smit & Singels, 2006). The partitioning of photosynthates between stems, leaves and roots was calibrated at each site against LAI and yield measurements for the first year (CH 2016;BR 2005BR -2006 ; Table S1). Where root biomass data were absent (CH), we assumed an average below to aboveground biomass ratio of 0.12 ± 0.08 in line with published estimates (Laclau & Laclau, 2009;Otto et al., 2009). We evaluated SPA carbon and water flux estimates against field estimates of stalk biomass, sucrose yield, LAI and soil moisture from specific years not used in model calibration (CH 2017(CH -2018BR 2006BR -2007. We further validated model performance against independent field measurements not used in model calibration; (i) net ecosystem exchange (BR only) and (ii) leaf-level gas exchange (CH only).We compared model and field time series using mean, standard deviation, R 2 and RMSE where appropriate. Data used in model calibration and validation are summarized in Table S2.

| Model simulations
We conducted a series of model simulations to quantify the impact of climate, and individual climate drivers on sugarcane yield. First, we forced the calibrated model at the CH and BR sites with the historical  and projected (2015-2100) climate data. From these model simulations we retrieved estimated sugarcane yield (Y) for each year (n = 1103; CH 86 years × 6 future climate projections, plus 39 years historical climate; BR 85 years × 6 future climate projections, plus 38 years historical climate). Model simulations using historical climates serve as a baseline, for comparing shifts in average crop yield, and crop yield variability F I G U R E 1 Mean annual climate across historical and future climate scenarios retrieved for the Brazil and China field site. Historical estimates are derived from the MERRA-2 re-analysis product. Future scenario estimates, relating to different levels of radiative forcing are derived from CMIP6 (MOHC-UKESM1.0) simulations, bias corrected against MERRA-2 estimates for the years across which both products are available under future climates. Growing season length at BR is greater than 1 year, therefore the number of simulated yields is 1 year less than CH for each climate projection. Next, to quantify the impact of individual climate drivers on sugarcane yield, we generated synthetic climate forcings ( Figure S7). Within the historical climate data set , we computed the mean value for each climate driver, for each year (i.e. temperature, CO 2 , precipitation, PAR and VPD). For each climate driver, we calculated the median of the annual means. We identified the annual time series (on an hourly time-step) corresponding to the median for each climate driver. Seasonality in the median year time series did not differ widely from other years (see Data S1). For each climate projection, we substituted, one by one, the temperature, CO 2 , precipitation, PAR and VPD, with that of the median time series. We then ran the model with the new synthetic climate forcings (n = 5515; CH 5 climate drivers × 86 years × 6 future climate projections, plus 5 climate drivers × 39 years historical climate; BR 5 climate drivers × 85 years × 6 future climate projections, plus 5 climate drivers × 38 years historical climate) and retrieved the yield estimate for each alternation (Y tair , Y CO 2 , Y ppt , Y PAR and Y VPD ).

| Model output analysis
We compared variation in simulated sugarcane yield by plotting yield probability density functions for each site and for each climate projection. We then computed the descriptive statistics for each distribution (i.e. mean, median, standard deviation, 10th percentile, 90th percentile, minima and maxima; Q1a). We compared the difference in median yield between climate projections using a Mann-Whitney U test.
We quantified the relative importance of individual climate drivers in determining sugarcane yield (Q1b). We computed the effect size of individual climate drivers (temperature, CO 2 , precipitation, PAR and VPD) on sugarcane yield. For each year, under each climate projection, at each site, we calculated the difference (%) between Y, and Y tair , Y CO 2 , Y ppt , Y PAR and Y VPD (i.e. the climate driver yield effect). We compared the difference in yield effect between climate drivers using box plots (Q1b).
To investigate the effect of climate on yield lows, we compared the yield lowest decile as a percentage of the climate projection median (Q2a). We then identified years in which simulated yield was within the lowest decile and compared the effect of individual climate drivers on sugarcane yield using box plots (Q2b). Where appropriate, we also present further details on simulated carbon and water fluxes, to show which modelled processes drive yield variability. Modelled gross primary productivity (GPP), weighted soil water potential, evapotranspiration and LAI were retrieved from model simulations. The effect of individual climate drivers on modelled processes was computed as the difference between simulated values under projected climate time series and simulations where the respective individual climate driver was substituted with the historical median.

| Model validation for carbon and water cycles
Observed stalk biomass (BR) and partitioning to sucrose (CH) were well captured by SPA. Model estimates of sugar yield at the CH field site were within 1 SD of field observations for validation years (field observation 2017 = 656 ± 89 g C m −2 , SPA 2017 = 611 g C m −2 ; field observation 2018 = 625 ± 68 g C m −2 , model 2018 = 571 g C m −2 ). Stalk dry biomass at the BR field site was significantly correlated with SPA estimates during the validation year (R 2 = 0.97, p < 0.001, RMSE = 51 g C m −2 ). In addition, pre-harvest stalk dry biomass estimates at BR were within 1 standard deviation of field observations for the valida- SPA successfully simulated observed LAI across the growing season. Observed LAI increased during the early and mid-season, and declined before harvest (Figure 3). Temporal variation in LAI was well captured by SPA during validation years (CH, R 2 = 0.86, p < 0.001, RMSE = 0.44; BR, R 2 = 0.84, p < 0.001, RMSE = 0.55). SPA LAI estimates during the validation years were 8% higher than field observations at the CH site, and 2% higher at the BR (for LAI estimates >1 m 2 m −2 ).
Simulated carbon and water fluxes were consistent with independent field measurements not used in model calibration. SPA captured well the observed temporal variation in NEE at BR ( Figure S4a). SPA estimates of NEE were significantly correlated with field observations (measured at the BR field site only) during the calibration year (R 2 = 0.57, p < 0.001, RMSE = 3.1 g C m −2 day −1 ). During the validation year, SPA and field estimates remained significantly correlated; however, the model explained less variation in observations (R 2 = 0.29, p < 0.001, RMSE = 2.5 g C m −2 day −1 ). Mean residual NEE (measured NEE minus modelled NEE) was 1.5 g C m −2 day −1 during the calibration year, and 2.8 g C m −2 day −1 during the validation year. SPA also captured observed temporal variation in latent energy ( Figure S4c). SPA-simulated latent heat flux at BR was significantly correlated with field measurements during the calibration (R 2 = 0.69, p < 0.001, RMSE = 1.9 MJ m −2 day −1 ) and validation year (R 2 = 0.61, p < 0.001, RMSE = 2.2 MJ m −2 day −1 ). Mean residual latent heat flux (measured latent heat flux minus modelled latent heat flux) was 0.30 MJ m −2 day −1 during the calibration year, and −0.44 MJ m −2 day −1 during the validation year ( Figure S4b).

| Sensitivity of sugarcane yield to climate
In China there was a mix of negative and positive responses in simulated yield under future climate projections, relative to historical yields (Figure 4a). Median yield increased most under the SSP119 projection (+4%), and decreased most under the SSP370 projection (−3%) compared to the median yield under the historical climate of 762 g C m −2 (Table 1). Changes in mean yield were consistent with changes in median yield, increasing most under the SSP119 projection (+4%), and decreasing the most under the SSP370 projection (−3%), in comparison to mean yield under the historical climate (764 g C m −2 ). Coefficient of variation across years was similar for the historical and forecast projections (historical = 2.9%; forecast projections = 2.5%-3.3%). Simulated yield under future climates differed significantly from that under historic climate for all projections excluding SSP245 and SSP434 (Table S3).
At BR, yield largely decreased under future climate projections (relative to historical climates; Figure 4b). Median and mean yield declined the most under SSP585 (−8% and −10% respectively), and the least under SSP126 (0% and −1% respectively). Yield coefficient of variation was much larger in the forecasts than the historic value of 3.8% (forecast projections = 5.1%-10.7%), for instance, yield coefficient of variation more than twice as large in SSP119, SSP245, F I G U R E 2 Observed (data point) and simulated (line) sucrose accumulation at the CH field site (a) and stalk dry biomass at the BR field site (b). Error bars represent observed standard deviation. Sucrose and stalk dry biomass estimates prior to the vertical dashed line were used in model calibration. Estimates occurring after the dashed line were used in model validation F I G U R E 3 Observed (data point) and simulated (line) leaf area index at the CH field site (a) and the BR field site (b). Error bars represent observed standard deviation. LAI estimates prior to the vertical dashed line were used in model calibration. LAI estimates occurring after the dashed line were used in model validation SSP434 and SSP 585. This increase in variance was linked to a substantive drop in the minimum yields under many projections. The average reduction in the minimum yield was 26% across the six projections (Table 1). Simulated yield under future climates differed significantly from that under the historic climate for all projections excluding SSP126 (Table S3).

| Role of individual climate drivers in determining sugarcane yield variation
Across future climate projections in China, the yield effect of individual climate forcings was relatively modest. PAR had the largest effect on sugarcane yield at CH ( Figure 5), but increases in PAR under the SSP119 projection ( Figure 1) had a positive effect on median sugarcane yield of just 2.6%. Under the SSP245, SSP370 and SSP585 projection, lower and more variable PAR (Figure 1) had a negative effect on simulated sugarcane yield (x = −2.4%, x = −5.1%, x = −0.39%, respectively, where x is the median yield effect; Figure 5). These shifts in yield are of the same magnitude of the overall climate effect for each projection (Table 1), indicating the PAR drives the climate response in CH. The exception is SSP 585, where the rise in CO 2 concentration also make a contribution to increased yield, although the effect is small. Overall the effect of rising CO 2 at CH was to increase median yield by 1.0%-3.1%.

F I G U R E 4
The probability density function of simulated sugarcane yield across historical and future climate scenarios at the (a) China and (b) Brazil field site. Climate forcing used to drive the model was retrieved from the MERRA-2 re-analysis product (historical climate) and CMIP6 (MOHC-UKESM1.0) simulations (future climate scenarios) T A B L E 1 Descriptive statistics on simulated sugarcane yield (g C m −2 ) under historic  and future climate projections (2015-2100). To generate yield estimates, SPA was forced with climate data from the MERRA-2 re-analysis product (historic climates), and biascorrected CMIP6 (MOHC-UKESM1.0) scenarios The effects of individual forcings on yield forecasts in Brazil are much stronger than in China ( Figure 5). Temperature had the largest, and mostly negative, effect on sugarcane yield at the BR site, VPD and precipitation also contributed to differences in yield, mostly negative. Under the SSP119 projection, the median temperature yield effect was −5.7% and was similarly high across SSP245, SSP370 and SSP585 (x = −11.7%; x = −4.4%; x = −18.1% respectively), reflecting the shift in mean annual temperature (Figure 1). The median precipitation yield effect was most pronounced in SSP585 (x = −5.3%), while the median VPD yield effect was highest under SSP245 and SSP585 (x = −5.9%, and x = −9.0% respectively), reflecting the climate trend across projections (Figure 1). The effect of changes in atmospheric CO 2 on median-simulated yield at the BR was positive, ranging from 3.7% to 11.4%.
The overall impact across the full SSP forcings at BR was to reduce median yield by 1%-9% (Table 1). Compared to the individual effects revealed in the experiments noted above, we can conclude that the combined effects of temperature, VPD and precipitation are not simply additive. Even with the compensating effect of rising CO 2 , the interactive effect of changing climate on yield is less than the sum of the individual components.

| Climate sensitivity of yields lows
At the CH site, the climate sensitivity of yield lows was similar in sign and magnitude ( Figure S8) to the response of median yield ( Figure 5). There was no sensitivity of yield variance to the climate projections, nor a quantifiable shift in the value of the lower decile of yield relative to the median (Table 1). For historic and future climates, the lower decile of yields was 3%-4% lower than the median yield, reflecting consistent productivity under varied climate.
At the BR site, yield lows were more extreme under all of the future projections tested compared to historical yields ( Figure S8). For the historic case, the lower decile of yields was 5% less than the median. In the forecast cases, this reduction for the lower decile was 7%-18% compared to predicted median yields. Thus, the lower decile yield decreased under all future climate projections relative to simulated historical yields, and was more sensitive to climate shifts than median yield ( Figure S8). For instance, the temperature sensitivity of median yield reduction varied from 4% to 18% across four selected SSPs ( Figure 5) while lower decile yield reduction varied from 25% to 42%.

| Role of individual climate drivers in determining low yield years
At the CH site, PAR was the major driver of yield lows (i.e. lower decile), but precipitation caused the most extreme yield reductions ( Figure S8). The PAR yield effect was highest under projection SSP370 (5th Percentile, i.e. the median of the lower decile = −8.3%; minimum yield-effect = −9.8%), when PAR was lowest (Figure 1). However, the most extreme yield effect was associated with precipitation under the SSP119 projection (minimum yieldeffect = −11.0%). Projected precipitation reached extreme lows under both the SSP119 and SSP585 projections, but F I G U R E 5 The relative effect of variation in temperature, CO 2 concentration, precipitation, PAR and VPD on sugarcane yield under four future climate scenarios at the Chinese (orange) and Brazilian (grey) sites. To quantify the yield effect, we retrieved yield estimates (Y) from model simulations ran for each year, in each climate scenario. We then repeated the model simulations, substituting the temperature, CO 2 , precipitation, PAR or VPD time series for that which represents the historical median, and retrieved the yield estimate (Y tair , Y CO 2 , Y ppt , Y PAR and Y VPD respectively). The yield effect was then calculated as the difference between yields (i.e. Y and Y tair etc. relative to Y). Boxplots show the distribution (10th, 25th, 50th, 75th and 90th percentiles, plus outliers) of yield effects across the multiple simulations | 417 FLACK-PRAIN et AL. water constraints were alleviated via higher atmospheric CO 2 in the latter (Figure 1).

| Determinants of the disparity in climate sensitivity between Brazil and China
Yield sensitivity to projected climate change was higher at Brazil, compared to China, because the Brazil site was more prone to higher temperatures and greater soil moisture stress (Figure 1). In contrast, PAR proved relatively more important in determining yield at the China site because moisture stress was a lesser constraint and forecast temperature was lower relative to BR. The relative sensitivity of yield to changes in CO 2 was significantly higher at the Brazil site in comparison to China (p < 0.001; Figure S9). Yield sensitivity to CO 2 was affected by moisture stress (Figure 6a). The CO 2 fertilization response differed significantly between precipitation levels (p < 0.001). The CO 2 fertilization effect was strongest when moisture stress was highest (mean CO 2 yield effect by precipitation level <1000 mm = 59.4 ± 37.2 g C m −2 ; 1000-1800 mm = 38.0 ± 35.3 g C m −2 ; >1800 mm = 24.3 ± 26.5 g C m −2 ). The relative sensitivity of yield to changes in temperature was also heightened at the Brazil site in comparison to China (slope = −30.3 ± 1.5 g C m −2°C −1 and −2.5 ± 0.2 g C m −2°C−1 respectively). As with CO 2 , the yield response to air temperature was dependent on moisture stress, and differed significantly across precipitation ranges (p < 0.001; Figure 6b). The temperature yield effect was greatest at low precipitation (mean air temperature yield effect by precipitation level <1000 mm = −153.1 ± 100.6 g C m −2 ; 1000-1800 mm = −40.8 ± 72.9 g C m −2 ; >1800 mm = −2.4 ± 34.4 g C m −2 ). The relative sensitivity of yield to changing mean PAR differed significantly between China and Brazil (p < 0.001; slope = 0.85 ± 0.05 g C µmol −1 s −1 and 0.52 ± 0.06 g C µmol −1 s −1 respectively; Figure S9), and was dependent on precipitation level (p < 0.001; Figure 6c; mean PAR yield effect by precipitation level <1000 mm = −1.4 ± 15.2 g C m −2 ; 1000-1800 mm = −10.6 ± 25.4 g C m −2 ; >1800 mm = −23.1 ± 25.6 g C m −2 ). Yield was relatively insensitive to changes in PAR when the system was water limited (i.e. precipitation <1000 mm). The relative sensitivity of yield to changes in VPD differed significantly between sites (p < 0.001) and was heightened in Brazil. The VPD yield effect was strengthened when precipitation was low (mean VPD yield effect by precipitation level <1000 mm = −98.7 ± 82.6 g C m −2 ; 1000-1800 mm = −21.1 ± 42.5 g C m −2 ; >1800 mm = −16.3 ± 16.3 g C m −2 ).

F I G U R E 6
The effect of precipitation on the response of sugarcane yield to variation in (a) atmospheric CO 2 concentration, (b) temperature, (c) PAR and (d) VPD. Low (yellow <1000 mm), mid (grey 1000-1800 mm) and high (blue >1800 mm) levels of precipitation drive differences in climate yield effects. To quantify the yield effect, we retrieved yield estimates (Y) from model simulations ran for each year, in each climate scenario at the Brazil site. We then repeated the model simulations, substituting the CO 2 , temperature and PAR for that which represents the historical median, and retrieved the yield estimate (Y CO 2 , Y tair and Y PAR respectively). The yield effect was then calculated as the difference between yields (i.e. Y and Y tair etc.)

| Eco-physiological processes underpinning high yield climate sensitivity at Brazil
Low precipitation and high VPD limited crop eco-physiological processes in SPA through adjustments to plant water supply and atmospheric water demand respectively. When soil water input was low (i.e. low total precipitation), mean-weighted soil water potential across the growing season was reduced (Figure 7g, R 2 = 0.32, p < 0.001). As a result, total GPP and mean LAI declined (Figure 7a, R 2 = 0.18, p < 0.001; Figure 7d, R 2 = 0.17, p < 0.001). Total evapotranspiration across the growing season increased sharply with increases in VPD (Figure 7k, R 2 = 0.48, p < 0.001). The increase in evaporative demand was met by drawing more water from belowground, reducing meanweighted soil water potential (Figure 7h, R 2 = 0.33, p < 0.001).

F I G U R E 7
The response of modelled total GPP, mean LAI, mean-weighted soil water potential, total evapotranspiration and days to reach maturity (temperature only), to differences in precipitation, VPD and air temperature at the Brazil site. Climate effects on eco-physiological processes are derived from model experiments in which model estimates for a given process (P) were retrieved from model simulations ran for each year, in each climate scenario. We then repeated the model simulations, substituting the precipitation, VPD and air temperature for that which represents the historical median, and retrieved the process estimate (P ppt , P VPD and P tair respectively). The effect of each climate driver on each eco-physiological process was then calculated as the difference between model estimates (i.e. P and P tair etc.). Ecophysiological processes include climate effects on GPP (a-c), mean LAI (d-f), weighted soil water potential (g-i), evapotranspiration (j-l) and development stage (m) | 419

FLACK-PRAIN et AL.
Water stress caused by high mean VPD led to a decrease in both total GPP and mean LAI (Figure 7b, R 2 = 0.62, p < 0.001; Figure 7e, R 2 = 0.44, p < 0.001). High variability in the effect of individual climate drivers on crop eco-physiological processes is a result of simultaneous variation in other climate drivers. Increases in temperature enhanced the rate of crop development, but reduced total C assimilated. Total GPP declined as mean temperature increased (Figure 7c, R 2 = 0.70, p < 0.001), as did total evapotranspiration (Figure 7l, R 2 = 0.55, p < 0.001). Mean LAI across the growing season was increasingly reduced as mean temperature increased up to ~25°C, after which changes in mean LAI plateaued (Figure 7f, R 2 = 0.64, p < 0.001). The number of days taken to reach maturity declined as mean temperature increased (Figure 7m, R 2 = 0.94, p < 0.001).
The impact of air temperature on crop eco-physiology was effected through two key processes, namely (i) higher temperature enhanced the rate of crop and thus canopy development, and (ii) leaf turnover rate increased with ageing as a function of crop development rate. As a consequence, while early season LAI was typically higher under projected increases in air temperature (Figure 8b), late season LAI was reduced due to accelerated leaf turnover (Figure 8d). As a consequence, GPP was reduced (x = 17% or 678 g C m −2 ; Figure 8a). Despite no significant change in PAR incident on the canopy (x = 18.1 ± 9.4 vs. that under historical median air temperature x = 17.8 ± 9.6; p = 0.6), radiation use efficiency was also reduced (x = 19% or 0.11 g C MJ −1 ; Figure  8b), as a result of lower GPP.

| DISCUSSION
Our results demonstrate important differences in the sensitivity of median sugarcane yield to future climate change projections between two major producing regions in Brazil and China (CH −2.9% to +3.8%; BR −8.2% to +0.2%).
The key drivers of median sugarcane yield sensitivity differ between sites, reflecting local climatic constraints to eco-physiological processes (CH PAR mean yield-effect = −1.2%; BR temperature mean yield-effect = −7.9%). The sensitivity of yield lows (i.e. lowest decile) to future climate projections varies between sites and is greater for Brazil (CH −3.1 to +2.5%; BR −18.2 to −2.4%). The most extreme declines in simulated yield are a result of moisture stress and temperature (CH precipitation yield-effect −11.0%; BR VPD yield-effect = −78.3%; BR temperature yield effect = −71.1%). Yield sensitivity to elevated CO 2 , temperature, PAR and VPD was dependent on moisture stress (Figure 6). At Brazil, higher temperatures reduced crop yield despite accelerated crop development (Figure F I G U R E 8 Time series of simulated (a) GPP, (b) LAI, (c) radiation use efficiency (computed as total daily GPP divided by total daily PAR incident on the canopy) and (d) cumulative leaf litterfall across the growing season at the Brazil site. Ribbons represent the 5th and 95th percentile limits and lines represent the median, for projected climate time series (grey) and model simulations where temperature was substituted with that of the historical median (yellow) 7), due to the simultaneous acceleration of leaf turnover driving reduced C assimilation (Figure 8). Potential mitigation strategies and limitations to the approach used are discussed in Data S1.

| Divergence in the sensitivity of median sugarcane yield to future climate
Our results predict that the climate sensitivity of median sugarcane yield is higher for São Paulo, Brazil compared to Guangxi, China. However, low crop sensitivity to future climate change in Guangxi (<4%) is in contrast to previous estimates for the region. Using the APSIM-Sugarcane model, Ruan et al. (2018) predicted increases in sucrose dry matter of up to 46% across Guangxi, under the RCP4.5 and RCP8.5 projections. Likewise, our SPA-predicted yield declines of up to 8% in São Paulo are contrary to existing projections. Using the Agro-Ecological Zone model (AZM; Doorenbos & Kassam, 1979), dos Santos and Sentelhas (2014) predicted regional increases in sugarcane yield of 59%-82% by 2090.

| The impact of future climate change on yield lows varies between sites
Inter-annual variation in yield increased in São Paulo, but not Guangxi, under future climate change (Table 1; 10th percentile as a % of the median). Current meta-analyses focusing on C 4 crops, such as maize and sorghum, suggest yield variability is likely to increase in the future (Challinor et al., 2014). But other studies report a decrease in C 4 crop yield variability under future climate change as a consequence of elevated CO 2 (Marin et al., 2013;Walker & Schulze, 2008). Our results predict regional differences in yield variability are likely. We also predict more extreme yields in the bottom decile for São Paulo, but not for Guangxi under future climate. Everingham et al. (2015) report a similar divergence in the sensitivity of the 25th percentile to future climate projections for sites across Australia. Given the spatial heterogeneity in inter-annual variability trends, local-level predictions are needed to assess suitable climate adaptation and mitigation strategies, which can significantly reduce yield variation (Cardozo et al., 2018).

| Sugarcane yield at Guangxi is light limited under future climate projections
PAR is the principal driver of sugarcane yield variation in Guangxi ( Figure 5), indicating that radiation is the principal limit on yield in this relatively cloudy region (Figure 1). Our findings are supported by a yield gap analysis conducted in Southern China, which reported PAR as a key constraint to yield potential (correlation coefficient = 0.43; Zu et al., 2018). The yield effect of PAR at São Paulo is of a similar magnitude to that in Guangxi ( Figure 5); however, its relative importance is lower in comparison to other climate drivers. In line with our findings, 0.3%-6% shifts in sugarcane yield have been reported in response to projected changes in PAR for model experiments in Brazil and Swaziland, with other climate drivers such as temperature proving relatively more important (Gouvêa et al., 2009;Knox et al., 2010). With respect to crop eco-physiological processes, the enhanced sensitivity of simulated GPP to PAR when moisture stress is low (Figure 6c) is supported by inter-annual differences in eddy flux C estimates at the Brazil site (Cabral et al., 2013). The authors report higher sensitivity of gross ecosystem productivity to PAR when the humidity deficit was lower.

| Increases in temperature drive declines in sugarcane yield at São Paulo
Projected increases in temperature will have a substantive negative effect on sugarcane yield in São Paulo ( Figure 5). Similarly, Sonkar et al. (2019) report a 3%-15% decrease in stalk fresh mass and sucrose mass in response to a 1°C-4°C warming in Northern India using the Canegro-sugarcane model (though the effect was negated by higher atmospheric CO 2 ). Furthermore, Linnenluecke et al. (2020) use historical data to show the negative effect of mean daily maximum temperature increases on sugarcane output across Australia's main sugar growing regions.
However, our results are counter to a number of other crop model studies, in which accelerated canopy development under rising temperatures (across a range of baseline temperature) leads to greater light interception, and therefore higher yield (Gouvêa et al., 2009;Jones et al., 2015;Knox et al., 2010;Marin et al., 2013;Ruan et al., 2018;Singels et al., 2014). In SPA, as in both DSSAT-Canegro and APSIM-Sugarcane, temperature (or thermal time) drives crop development. Across models, leaf senescence is induced by shading and moisture stress (Singels et al., 2008). However, in SPA, leaf turnover is also dependent on ageing, which is a function of crop development rate ( Van Laar et al., 1997). Therefore accelerated crop development (as a consequence of increases in thermal time) led to increased late season leaf turnover rates (Figure 8d). As a result, simulated mean LAI across the growing season decreased with increases in temperature (Figure 7f). Simulated GPP declined in parallel (Figure 7c) driving lower yields.
We highlight that the ability of SPA to capture field LAI estimates during model validation years was strong (CH, R 2 = 0.86, RMSE = 0.44; BR, R 2 = 0.84, RMSE = 0.55). However the temperature ranges across which the model was tested were understandably limited. Previously modelling studies have struggled to capture LAI dynamics for sugarcane sites across a range of mean annual temperatures (Marin et al., 2015; APSIM-Sugarcane r = 0.45, RMSE = 1.35 m 2 m −2 ; DSSAT/ CANEGRO r = 0.67, RMSE = 1.10 m 2 m −2 ). Although current field evidence does support a broad increase in the number of dead leaves with thermal time (Smit & Singels, 2006), further work is required to gain a more detailed understanding of phenological responses to temperature beyond accelerated canopy development (Robertson & Donaldson, 1998).

| Yield sensitivity to CO 2 is dependent on moisture stress
Sugarcane yield is more sensitive to elevated CO 2 in São Paulo in comparison to Guangxi. The disparity between sites is a result of differences in water limitation, whereby São Paulo is subject to greater water constraints (Figure 1). Other modelling studies similarly report a positive impact of elevated CO 2 on sugarcane yield, but likewise differ in the extent of simulated CO 2 effects (Gouvêa et al., 2009;Singels et al., 2014;Sonkar et al., 2019). It has been suggested that beyond water-stress alleviation, elevated CO 2 has little impact on sugarcane (Stokes et al., 2016), as observed in other C 4 crops (Leakey et al., 2006;Ottman et al., 2001). Furthermore, modelling efforts focusing on maize, another C 4 crop, similarly project low yield sensitivity to rising CO 2 concentrations (Bassu et al., 2014). However, experimental evidence on the contrary has also been presented. De Souza et al. (2008) report high sugarcane yield sensitivity to elevated CO 2 even when water is not limited in an open top chamber experiment. In addition, Vu and Allen (2009) report leaf carbon exchange rate increases of up to 5%-9% under non-drought stressed conditions (720 ppm vs. 360 ppm CO 2 ). Similarly, we found that in the absence of moisture stress, simulated GPP increased by 4%-7% under elevated CO 2 (710-730 ppm vs. 368 ppm) at the Guangxi site.

| Moisture stress and temperature drive the most extreme yield lows
The magnification of extreme low yields forecast for São Paulo state is linked to forecast changes in precipitation, temperature and VPD. Precipitation and VPD determine the development of moisture restrictions in photosynthesis in the SPA model, through adjusting soil moisture and atmospheric demand. The sensitivity to moisture stress is heightened as this part of Brazil is warmer, less humid and drier than Southern China (Figure 1). Sugarcane is especially vulnerable to moisture stress due to its high water requirements (Singh et al., 2007). Our results align with pan-tropical modelling and field studies which report lower sugarcane yield in response to reduced precipitation in rainfed systems de Carvalho et al., 2015;Jones et al., 2015;Santillán-Fernández et al., 2016;Sonkar et al., 2019). As in our study ( Figure 5; Figure S8), Marin et al. (2015) found that predicted crop yields across Brazil were more sensitive to precipitation reductions than increases (crop yield +12% under +30% precipitation; crop yield −30% under −30% precipitation, predicted using the APSIM-Sugarcane model).

| CONCLUSION
The impact of future climate change and climate extremes on sugarcane production varies substantively across two major sugarcane-producing regions. In Southern China, we present evidence that median yield and the variance of yield with climate extremes will not change over future decades. But in São Paulo state, Brazil, we present evidence that median yield will likely decline, and extreme climate effects will reduce minimum yields even more strongly. Thus in Brazil, sugarcane production has a much larger climate risk. We suggest that mitigation strategies which target local climate drivers are needed to guard global sugarcane production from future climate risks. Further work should focus on the climate response of sugarcane phenology, in particular the linkages between climate, development rate and leaf senescence.