Integrating Wheat Canopy Temperatures in Crop System Models

Crop system models are generally parametrized with daily air temperatures recorded at 1.5 or 2 m height. These data are not able to represent temperatures at the canopy level, which control crop growth, and the impact of heat stress on crop yield, which are modified by canopy characteristics and plant physiological processes Since such data are often not available and current simulation approaches are complex and/or based on unrealistic assumptions, new methods for integrating canopy temperatures in the framework of crop system models are needed. Based on a forward stepwise-based model selection procedure and quantile regression analyses, we developed empirical regression models to predict winter wheat canopy temperatures obtained from thermal infrared observations performed during four growing seasons for three irrigation levels. We used daily meteorological variables and the daily output data of a crop system model as covariates. The standard cross validation revealed a root mean square error (RMSE) of ~0.8  ̋C, 1.5–2  ̋C and 0.8–1.2  ̋C for estimating mean, maximum and minimum canopy temperature, respectively. Canopy temperature of both water-deficit and fully irrigated wheat plots significantly differed from air temperature. We suggest using locally calibrated empirical regression models of canopy temperature as a simple approach for including potentially amplifying or mitigating microclimatic effects on plant response to temperature stress in crop system models.


Introduction
The application of crop system models is crucial for predicting the impact of climate change on crop yields and for deriving management strategies to optimize and stabilize agricultural production.Models that are able to accurately simulate the impact of extreme temperatures on crop yields are required for realistic results from scenario analyses and for complementing experimental studies on the performance of different genotypes under stress.Agrometeorological variables, such as minimum and maximum air temperatures, are one of the model key components controlling the simulated production, especially for rainfed production systems [1].Since the availability of meteorological data at sub-daily resolutions is limited, models are typically parametrized with data of a temporal resolution of one day.Thus, in general, models do not account for the diurnal cycle of input variables.
Currently, widely used crop system models are parametrized with air temperatures recorded at 1.5 or 2 m height or simulated daily maximum and minimum temperatures.Air temperature data are used to simulate crop development and crop physiological processes, such as photosynthesis, transpiration and light-use efficiency.Corresponding model routines use crop-specific temperature thresholds and temperature response functions.It is, however, well known that air temperatures may strongly differ from near-surface temperatures, which are modified by the soil and canopy physiological and structural characteristics.In particular, evaporation and transpiration rates affect microclimatic conditions [2][3][4], and, thus, might either amplify or mitigate the impact of temperature stress.Cooling down of the air surrounding a plant as a consequence of transpiration has been recognized as an efficient heat tolerance strategy [5].For a short, sheltered or dense crop canopy, the microclimatic conditions can be largely decoupled from overlying atmospheric conditions [6].Knowledge of the canopy microclimate, controlled by canopy aerodynamic properties and stomatal regulations, is, therefore, crucial for an adequate crop management [7].Crop system models simulating rice production compute spike temperatures for predicting the impact of heat on spike fertility and the resulting yield in order to account for the cooling effect of transpiration [8].Contrastingly, widely applied wheat, barley or rye models use strongly simplified approaches, such as critical air temperature thresholds, to derive yield reduction or temperature stress factors (see [9] for a recent review).Ignoring canopy-air temperature differences introduces uncertainty in crop yield responses [10] and results in unrealistic parametrizations of the temperature thresholds to capture the effects of heat stress [11,12].Most models are highly sensitive to input temperatures [13,14].Consequently, uncertainty in input temperatures leads to significant variability of model results and probably unrealistic results from scenario-based analyses of the impact of climate change on plant productivity.Using canopy temperatures in addition to air temperatures decreases such uncertainties and, in consequence, might improve the simulation of the effect of temperature as a defining, limiting and reducing factor in crop system models [5,[10][11][12]15,16].This requires, however, a systematic evaluation and adaptation of model routines [12] which are often based on empirical data, optimized for using air temperature, and often not able to account for the dependence of surface temperatures on crop type, developmental stage and canopy architecture.Further, using currently available model routines the data derived from radiometric measurements might overestimate the impact of drought and heat stress on crop yield since they do not capture the complete canopy profile [11].Assuming that the physiological basis of the complex interrelations between canopy temperature with stomatal conductance and plant water status as well as their expression in the field [17] is well-understood, the use of canopy temperatures in dynamic models will probably greatly facilitate the identification and selection of higher yielding and/or more resilient crop genotypes (e.g., [18][19][20][21]).
The availability of canopy temperature data sets is limited.In Germany, only a few agro-climatological weather stations record temperatures at 0.2 m height.While satellite-derived surface temperatures provide large area coverage, their temporal resolution is strongly limited.The model-based simulation of plant and canopy temperature dynamics, which are influenced by several feedback mechanisms with the atmosphere and the soil, is complex.Variables needed for the parametrization of these models are the optical and aerodynamic properties of the individual plants, the soil and the canopy, plant physiological activity and meteorological variables, such as wind speed at the canopy level and cloud cover.The acquisition of such data requires manual measurements by trained personnel and/or high expenses for the installation of suitable instrumentation techniques, for instrument maintenance and calibration.Further, several of these observations are often not standardized.Few crop system models simulate these dynamics (see [12] for a recent review), e.g., by employing simplified mechanistic models to derive canopy temperatures, by solving the local energy balance or by coupling crop models with advanced land surface models [22,23].However, most approaches rely on unrealistic assumptions, such as an exact closure of the energy balance or the use of a constant fraction of energy for the soil heat flux, or require a large amount of input data at high temporal resolution and long computational time.An alternative approach is the use of simple empirical equations, based on daily environmental and meteorological data, to simulate site-and crop-specific daily mean, maximum and minimum canopy temperatures.Multiple linear regression models are a widely used empirical approach to simulate data which are not available from measurements.Although strongly simplifying reality, they provide a straightforward and easy to implement approach to approximate local microclimatic conditions and, thus, account for their effect on crop yields in crop system models.
In this study, we aim at developing and testing the use of empirical regression models to predict canopy temperature dynamics.Models are based on radiometric canopy temperature measurements from two research sites located in northern Germany performed for three different irrigation treatments (no irrigation, medium deficit irrigation and full irrigation) during four growing periods (2010, 2011, 2013 and 2014).More specifically, we aim at (I) developing empirical models to predict maximum, minimum and mean daily canopy temperatures based on meteorological data and environmental variables available from the output of a crop system model operated at a daily time step and (II) studying potential effects of using derived canopy temperatures as input temperature in crop system models.
We hypothesize that diurnal canopy temperatures significantly differ from air temperatures for both, irrigated and deficit-irrigated wheat canopies.Such temperature anomalies impact the simulation of crop development and physiological processes and, therefore, need to be considered in crop system models.The canopy temperature of wheat is directly related to its water status [24][25][26][27][28].We therefore hypothesize that environmental variables available from widely applied crop system models, such as daily transpiration and soil evaporation, allow for approximating wheat canopy temperatures.

Research Sites and Experimental Set-Up
The experiments were conducted in Northern Germany at (I) the Hohenschulen Experimental Farm ("HS", 10.0 ˝E, 54.3 ˝N, 30 m a.s.l.(above sea level)) of the Kiel University and (II) a research site operated by the Thünen Institute ("BS", 10.26 ˝E, 52.18 ˝N, 79 m a.s.l.) located close to Braunschweig.(I) At the Hohenschulen research site ("HS"), in the humid climate of NW Germany, mean annual air temperature accounts with 8.4 ˝C and total rainfall averages 750 mm annually [29].The soil can be described as a pseudogleyic sandy loam (luvisol).Winter wheat (Triticum aestivum) was planted under a mobile rainout shelter with irrigation system in the years 2010 (sowing date: 24 September 2009, sowing density 300 plants/m 2 , cv.Dekan), 2011 (sowing date: 24 October 2010, sowing density: 300 plants/m 2 , cv.Dekan) [30], 2013 (sowing date: 31 October 2012, sowing density: 360 plants/m 2 , cv.Batis) and 2014 (sowing date: 2 October 2013, sowing density 300 plants/m 2 , cv.Batis).Plants were watered with an overhead boom irrigation system in 2010 and 2011 and by drip irrigation in 2013 and 2014.Three different irrigation treatments (W0 = no irrigation, W1 = medium drought stress, W2 = fully irrigated) with four replications were applied (Table 1).Using the rain-out shelter, drought stress was induced from the beginning of April (after danger of ground frost ceased).Canopy height (CH) and leaf area index (LAI) were measured on a weekly basis (LAI2000 plant canopy analyzer, Li-Cor Inc., Lincoln, NE, USA).Regular observations of the developmental stage were performed using the scale proposed by [31].Volumetric water content was measured in five depths (weekly: 5 and 35 cm, fortnightly: 65, 85 and 105 cm) using the MiniTrase Time Domain Reflectometry (TDR) System (Soil Moisture Equipment Corp., Santa Barbara, CA, USA).Air temperature (T air ), relative humidity and wind speed were measured at a reference height of 2 m within the plots.Net radiation was estimated using NR (net radiometer) Lite net radiometers (Kipp and Zonen, Delft, The Netherlands) within one plot of each treatment.The canopy temperature was derived from infrared (IR) radiometer measurements (in 2010 and 2011: SI 111 and SI 211 sensors, Apogee Instruments, Logan, USA; in 2013 and 2014: IR120 sensors, Campbell Scientific, Logan, UT, USA).The sensors were placed 0.5 m above the canopy using a nadir-viewing angle [30].Note that the position of the sensors was shifted at regular intervals to account for changes in the canopy height.Plot size averages 3.5 m ˆ3.6 m (12.6 m 2 ).The sensors were operated at a frequency of 10 s, and the data were averaged and logged at an interval of 1, 10 and 15 min in the years 2010 and 2011, 2013 and 2014, respectively.(II) At the Braunschweig research site ("BS"), operated by the Thünen Institute, rainfall averages 599 mm annually and the mean air temperature accounts with 9.2 ˝C [32].The soil is a luvisol of a loamy sand texture in the plough horizon, followed by a sand and gravel layer of more than 3 m in size [33].Soil water content was derived from a Trime Pico 64 TDR system (Imko GmbH, Germany) at 0-20 cm and 20-40 cm depth.In 2013 (sowing date: 29 October 2012, sowing density: 380 plants/m 2 , cv.Batis), plants were watered by circle sprinklers and manual irrigation (hand shower).Wind speed, net radiation and air temperature were measured at a station operated by the German Weather Service (DWD) located within the research field.The canopy temperature was measured with an IR 120 sensor (Campbell Scientific) at a frequency of 10 s.Data were logged and averaged using an interval of 10 min.The sensors were placed 0.5 m above the mean canopy height using a nadir-viewing angle.Plot size averages 5 m ˆ4 m (20 m 2 ).For this study, meteorological and canopy temperature data were averaged to hourly and daily values.

Canopy Temperature Data
The data set comprises canopy temperatures recorded by the infrared thermometers at 0.5 m above the canopy (cf.Section 2.1) in the years 2010, 2011, 2013 and 2014 from three different irrigation treatments, where W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation (Table 1, Figure 1).
Data gaps of a maximum length of four values within the hourly averaged canopy temperature data were interpolated using a spline function (R statistical software, zoo package).Subsequently, for complete diurnal cycles (24 measurements available), hourly data were aggregated to minimum, mean and maximum daily values for each irrigation treatment.Approximately half of the available data were measured at W2 irrigation levels (~46.5%),~30% and ~23.5% were measured at W0 and W1 irrigation levels, respectively (Table 1).Note that the plot or site effects on canopy temperatures are beyond the scope of this study and not considered.We did not perform an emissivity correction of measured temperatures.

Crop System Modeling
In order to derive time series of environmental variables describing canopy architecture and evapotranspiration, which affect microclimatic conditions, and, thus, canopy temperature, a model of soil-plant-atmosphere interactions, implemented in the HUME modeling framework [34], was fitted to observed soil water contents.It was then used to calculate the following environmental variables at a daily time step: leaf area index (LAI), crop height (CH), potential evapotranspiration (ETP), potential transpiration (Tpot), potential soil evaporation (Epot), actual evapotranspiration (ETA), actual transpiration (Tact) and actual soil evaporation (Eact).
The model framework comprises of submodels for calculating plant development, plant water uptake, evapotranspiration and vertical soil water transport.Basic approaches of the submodules can be found in [30,35,36].We refer to [30] for a detailed description of the model parametrization and equations.Plant growth data are derived from the interpolation of weekly LAI and crop height measurements.In a simple modeling approach the root distribution (RD) is calculated from the root length and depth growth simulations based on temperature sums [35,37] assuming a maximum total root length of 150 cm at anthesis and a maximum root depth of 130 cm.The higher root density in the upper soil layers of irrigated plots was considered by plot-wise adjustment of the root length distribution parameter [30].The vertical soil water transport is calculated by a numerical solution (using a variable internal time step length) of the water-content based formulation of the Richards equation.The submodel uses the functions published by [38] and revised by [39] to characterize the relation between soil diffusivity of water and soil water content (SWC) and the relation between unsaturated hydraulic conductivity and SWC.
The simulation of ETP is based on the Penman-Monteith equation [40] with aerodynamic resistances calculated from the approach presented by [41].Epot is derived from the fraction of global radiation reaching the soil surface, which is controlled by the LAI.Potential soil evaporation is reduced to Eact if the soil water potential values in the top layer undergo a critical threshold value of −2 × 10 −3 MPa in the top soil layer [42].Potential transpiration (Tpot) is the difference between the ETP and the sum of interception evaporation and Epot [30].The soil water uptake by the plants from each layer, expressed as a layer sink term [43], sums up to actual transpiration (Tact).The sink term of each layer is modeled by distributing Tpot to the rooted layers.The partitioning is controlled by the root length in each layer, modified by a root water uptake competition factor [44] and reduced by a layer-specific reduction factor depending on the layer-specific soil water potential [37].For the model parametrization literature values, measurements and results from an optimization method based on the Levenberg-Marquard algorithm (implemented in the modeling framework) are used.

Crop System Modeling
In order to derive time series of environmental variables describing canopy architecture and evapotranspiration, which affect microclimatic conditions, and, thus, canopy temperature, a model of soil-plant-atmosphere interactions, implemented in the HUME modeling framework [34], was fitted to observed soil water contents.It was then used to calculate the following environmental variables at a daily time step: leaf area index (LAI), crop height (CH), potential evapotranspiration (ETP), potential transpiration (Tpot), potential soil evaporation (Epot), actual evapotranspiration (ETA), actual transpiration (Tact) and actual soil evaporation (Eact).
The model framework comprises of submodels for calculating plant development, plant water uptake, evapotranspiration and vertical soil water transport.Basic approaches of the submodules can be found in [30,35,36].We refer to [30] for a detailed description of the model parametrization and equations.Plant growth data are derived from the interpolation of weekly LAI and crop height measurements.In a simple modeling approach the root distribution (RD) is calculated from the root length and depth growth simulations based on temperature sums [35,37] assuming a maximum total root length of 150 cm at anthesis and a maximum root depth of 130 cm.The higher root density in the upper soil layers of irrigated plots was considered by plot-wise adjustment of the root length distribution parameter [30].The vertical soil water transport is calculated by a numerical solution (using a variable internal time step length) of the water-content based formulation of the Richards equation.The submodel uses the functions published by [38] and revised by [39] to characterize the relation between soil diffusivity of water and soil water content (SWC) and the relation between unsaturated hydraulic conductivity and SWC.
The simulation of ETP is based on the Penman-Monteith equation [40] with aerodynamic resistances calculated from the approach presented by [41].Epot is derived from the fraction of global radiation reaching the soil surface, which is controlled by the LAI.Potential soil evaporation is reduced to Eact if the soil water potential values in the top layer undergo a critical threshold value of ´2 ˆ10 ´3 MPa in the top soil layer [42].Potential transpiration (Tpot) is the difference between the ETP and the sum of interception evaporation and Epot [30].The soil water uptake by the plants from each layer, expressed as a layer sink term [43], sums up to actual transpiration (Tact).The sink term of each layer is modeled by distributing Tpot to the rooted layers.The partitioning is controlled by the root length in each layer, modified by a root water uptake competition factor [44] and reduced by a layer-specific reduction factor depending on the layer-specific soil water potential [37].For the model parametrization literature values, measurements and results from an optimization method based on the Levenberg-Marquard algorithm (implemented in the modeling framework) are used.During the model fitting and parametrization, the soil texture was chosen for each plot separately [30] to account for site heterogeneities.
Environmental and meteorological variables were selected to formulate simple regression models which are able to predict daily minimum, mean and maximum canopy temperatures (T c,min , T c,mean and T c,max ).Selected covariates are observed at standard meteorological stations (net radiation (Rn), total incoming radiation (Rint), air temperature (T air ), humidity (RH) and saturation vapor pressure deficit of the air (VPD), rainfall (R) and wind speed (W)) or commonly available from the output of crop system models (as specified in the following).Using the HUME model (cf.Section 2.3), we simulated time series data of crop height and leaf area index and calculated the main components of the plant-water-cycle, such as the total and plant available soil water content, ETA (= Eact + Tact), ETP, Tpot and Tact, Epot and Eact (cf.Section 2.3).We transformed absolute model results to relative values by formulating a range of ratios based on the available possible combinations of these variables, characterizing plant and soil water relations independent from the absolute amount of available water.Different formulations describing the evaporation (Eact/ETP and Eact/ETA, Eact/Epot) and the transpiration (Tact/ETP, Tact/ETA, Tact/Eact, Tact/Tpot) were applied as potential model covariates.Furthermore, we tested the suitability of the interpolated LAI, the natural logarithm of the LAI (LAI log ) and crop height (CH) values as predicting variable and tested for potential interaction effects between modeled environmental (as specified above) and meteorological data.
By using a random sampling strategy (CreateDataPartition function [45]) the whole data set (685 observation days, cf.Table 1) was partitioned in a training (50%) and a testing data subset (50%).Optimized combinations of covariates for a multiple linear regression model (maximum R 2 ) were automatically calculated using a forward stepwise-based selection procedure (FWDselect function [46]).Meteorological, as well as environmental data are generally strongly inter-correlated.We assessed the multicollinearity of the co-variables by using the variance inflation factors (VIF; corvif function [48]) and consider VIF ď3 as a benchmark [50].Model reformulation and optimization were based on a quantile regression (QR) analysis (qr function [49]).The quantile regression is based on the conditional quantiles of the response variable distributions, thus, offering a more complete view of possible causal relationships between variables and useful for ecological applications with a limited number of available and strongly interacting variables [51].Results from the quantile regression showed that the contributions of the selected independent variables to the conditional distribution of the canopy temperatures can vary significantly at different levels of canopy temperatures.Such patterns are hidden in ordinary multiple least-square regression model analyses.The predictive quality of the explanatory variables when applied to the training and testing data was identified at the 5% level of significance (p < 0.05).

Impact Study
To highlight the impact of using canopy instead of air temperatures in crop system models we performed an impact study for the observation years 2010 (20 April-27 July), 2011 (1 April-27 July), 2013 (30 April-19 July) and 2014 (17 April-19 July).T c,mean , T c,min and T c,max for a well-irrigated wheat plot (W2) and for a not irrigated wheat plot (W0) are computed using the multiple linear regression models presented in this study (cf.Section 2.4) using data from the HS research site (cf.Section 2.1).Most crop system models apply air temperature-based "heat units" and use temperature thresholds to simulate the impact of temperatures on crop development, crop physiological processes and crop yield.To highlight the potential impact of the differences between air and crop temperatures on corresponding model results, we (I) calculate cumulative sum curves of the difference between modeled daily canopy temperature and the air temperature (∆T), (II) compute the number of days where modeled and measured canopy temperatures and air temperatures exceed temperature thresholds of 20 ˝C, 25 ˝C and 30 ˝C, (III) calculate extreme thermal unit (ETU) sums above the optimal temperature of 20 ˝C (according to [52]).
ETUs are calculated as the cumulative sums of the difference between the actual crop or air temperature, respectively, and the optimum temperature of 20 ˝C, where ETU = 0 if the actual temperature is below the optimum temperature [52].

Empirical Models of Daily Canopy Temperatures
In the following, we provide a brief description of the multiple linear regression models developed for estimating T c,mean , T c,min and T c,max (cf.Section 2.4) as a major result from our study.The relative importance of the covariates for the total estimated R 2 is described and discussed in Section 3.3.
Based on the results from the quantile regression analyses (data not shown) empirical canopy temperature models were fitted to observations during the pre-heading (corresponding to stages <50 of the Zadoks scale) and during the post-heading (>50 of the Zadoks scale) phase.During the pre-heading phase (phase I), the conditional quantile contribution of environmental and meteorological factors to the canopy temperature significantly varied from their contribution during the post-heading phase (phase II), leading to significant changes in estimated model coefficients (models of T c,min and T c,max ) and/or selection of model covariates (models of T c,min and T c,mean ).
For simulating the daily mean canopy temperature (T c,mean ), the covariates explaining more than 90% of the total variability are the mean daily air temperature ( ˝C), the incoming radiation (Rint (W¨m ´2)), the natural logarithm of the leaf area index (LAI log ), the vapor pressure deficit (VPD (hPa), only for phase I), the ratio of actual evaporation to ETP (Eact/ETP, only for phase II) and the product of the transpiration ratio, defined as the ratio of actual transpiration to the potential transpiration, and the VPD ((VPD*(Tact/Tpot)), only for phase II).We account for the changing contribution of the VPD, the Eact/ETP ratio and the product of transpiration ratio and VPD by using a dummy variable (DPhen), which equals 1 in the pre-heading and 0 in the post-heading phase (Table 2).
Table 2. Statistics of the multiple linear regression model developed using a forward stepwise-based selection procedure and quantile regression analyses for predicting T c,mean (T c,mean = mean daily canopy temperature) of a winter wheat canopy; SE = standard error, T air,mean = mean daily air temperature ( ˝C), Rint = incoming radiation (Wm ´2), LAI log = natural logarithm of the leaf area index, Dphen = phenology dummy variable with Dphen = 1 during phase I (pre-heading) and Dphen = 0 during phase II (post-heading, cf.Section 3.1), VPD = vapor pressure deficit (hPa), Tact/Tpot = transpiration ratio (Tact = actual and Tpot = potential transpiration (both in mm¨d ´1)), Eact/ETP = ratio of actual evaporation (Eact (mm¨d ´1)) to potential evapotranspiration (ETP (mm¨d ´1)).T c,max are modeled using Rint, T air,max , the natural logarithm of the leaf area index and the product of transpiration ratio and VPD.For the pre-and the post-heading phase significant changes in the model coefficients were calculated, thus, two model equations were derived (Table 3).For the pre-heading phase, T c,min can be approximated using T air,min and the canopy height (CH (m)).For the post-heading phase, it is simulated using the minimum diurnal air temperature, the VPD and the ratio of actual evaporation to ETP (Table 4).

Estimate
Table 3. Statistics of the multiple linear regression model developed using a forward stepwise-based selection procedure and quantile regression analyses for predicting T c,max (= maximum daily canopy temperature) of a winter wheat canopy; SE = standard error, T air,max = maximum daily air temperature ( ˝C), Rint = incoming radiation (Wm ´2), LAI log = natural logarithm of the leaf area index, VPD = vapor pressure deficit (hPa), Tact/Tpot = transpiration ratio (Tact = actual and Tpot = potential transpiration (both in mm¨d ´1)).

Variability of Canopy Surface Temperatures
In general, canopy temperatures of the drought treatment plots (W0) exceeded air temperatures (Figure 2), while the W1 and the W2 treatment canopy temperatures are lower (W1) and scatter more strongly (W2).Increased canopy temperatures result from the combined effect of decreased transpiration and evaporation levels of the plants and the underlying soil, respectively (cf.Section 3.3 for a discussion).The scattering of the W2 treatment data set can be mainly attributed the larger number of observations compared to the W1 and W0 treatments (cf.Table 1).Daily averaged surface temperatures (T c,mean ) range from ~8 ˝C to 26.5 ˝C (Figure 1) and daily averaged differences between canopy and air temperatures from ´4 ˝C to +3.8 ˝C (∆T mean ) (Figure 3).Compared with those of the irrigated plots, mean daily temperature differences between canopy and air (∆T mean ) values of the W0 treatment are clearly shifted towards higher values (Figure 3), indicating decreased stomatal aperture and lower transpiration levels under drought stress.The diurnal range of canopy temperatures measured at W0 plots clearly exceeds the diurnal range of air temperatures (Figure 4).
Compared with those of the irrigated plots, mean daily temperature differences between canopy and air (∆Tmean) values of the W0 treatment are clearly shifted towards higher values (Figure 3), indicating decreased stomatal aperture and lower transpiration levels under drought stress.The diurnal range of canopy temperatures measured at W0 plots clearly exceeds the diurnal range of air temperatures (Figure 4).

Predictive Ability of the Empirical Canopy Temperature Models
The standard cross validation revealed a root mean square error (RMSE) of ~0.8 °C for estimating Tc,mean, 2 °C (Phase I) and 1.5 °C (Phase II) for estimating Tc,max and 1.2 °C (Phase I) and 0.8 °C (Phase II) for estimating Tc,min, and could generalize well to the training data set (Table 5).
While the multiple linear regression model suggests a positive effect of radiation and air temperature, Tc,mean decreases with increasing LAI (Table 2).We calculated a negative slope for the VPD during the pre-heading phase and for the Eact/ETP ratio and the VPD-scaled transpiration ratio during the post-heading phase.The negative effect of the LAI and the VPD can be related to an increased transpiration cooling with increasing evaporative demand of the air and an increasing amount of transpiring plant tissues during the pre-heading stages.After the canopy has reached maximum leaf area and maximum height, during the post-heading stage, environmental conditions are characterized by decreased availability of soil water and increased air temperatures.Thus, for subsequent developmental stages, daily transpiration and evaporation help explain the variability of Tc,mean.An increasing fraction of actual evaporation (Eact) from the soil, e.g., after rain events, decreases the amount of energy available for heating the canopy (sensible heat fluxes).The increased latent heat flux arising from the soil further lowers the near-surface temperature.The negative slope of the VPD-scaled transpiration ratio can be related to the cooling effect of a higher amount of plant available soil water, weighted over the root distribution parameters and scaled with the evaporative demand of the air, on canopy temperatures.
The day-to-day variability of the maximum daily canopy temperatures can be simulated using incoming radiation, maximum daily air temperatures, the LAI and the VPD-weighed transpiration ratio (Table 3).Model coefficients differ for pre-and post-heading stages.The RMSE is higher for canopy temperatures simulated for the pre-heading stage (Table 5).The most likely cause is the higher overall range and day-to day variability of meteorological data typical for spring season conditions.

Predictive Ability of the Empirical Canopy Temperature Models
The standard cross validation revealed a root mean square error (RMSE) of ~0.8 ˝C for estimating T c,mean , 2 ˝C (Phase I) and 1.5 ˝C (Phase II) for estimating T c,max and 1.2 ˝C (Phase I) and 0.8 ˝C (Phase II) for estimating T c,min , and could generalize well to the training data set (Table 5).
While the multiple linear regression model suggests a positive effect of radiation and air temperature, T c,mean decreases with increasing LAI (Table 2).We calculated a negative slope for the VPD during the pre-heading phase and for the Eact/ETP ratio and the VPD-scaled transpiration ratio during the post-heading phase.The negative effect of the LAI and the VPD can be related to an increased transpiration cooling with increasing evaporative demand of the air and an increasing amount of transpiring plant tissues during the pre-heading stages.After the canopy has reached maximum leaf area and maximum height, during the post-heading stage, environmental conditions are characterized by decreased availability of soil water and increased air temperatures.Thus, for subsequent developmental stages, daily transpiration and evaporation help explain the variability of T c,mean .An increasing fraction of actual evaporation (Eact) from the soil, e.g., after rain events, decreases the amount of energy available for heating the canopy (sensible heat fluxes).The increased latent heat flux arising from the soil further lowers the near-surface temperature.The negative slope of the VPD-scaled transpiration ratio can be related to the cooling effect of a higher amount of plant available soil water, weighted over the root distribution parameters and scaled with the evaporative demand of the air, on canopy temperatures.
The day-to-day variability of the maximum daily canopy temperatures can be simulated using incoming radiation, maximum daily air temperatures, the LAI and the VPD-weighed transpiration ratio (Table 3).Model coefficients differ for pre-and post-heading stages.The RMSE is higher for canopy temperatures simulated for the pre-heading stage (Table 5).The most likely cause is the higher overall range and day-to day variability of meteorological data typical for spring season conditions.Table 5. Results of the multiple linear regression fit for all data and for the different irrigation treatments (W0-W2, cf.Section 2.1) using phenological subsets (I = pre-heading phase, II = post-heading phase) for modeling T c,mean (mean daily canopy temperature ( ˝C)), T c,max (maximum daily canopy temperature) and T c,min (minimum daily canopy temperature).R 2 and the root mean squared error (RMSE ( ˝C)) are given for the training (50% of all observations) and the testing data (50% of all observations).During the early vegetative phase, additional influences of soil temperatures and soil optical properties on radiometric temperatures, in particular at low LAI levels as observed in the drought treatment plots, cannot be ruled out.Short-term influences, such as the effect of rapid fluctuations of wind and radiation on maximum daily canopy temperatures which cannot be incorporated in the model if operated at a daily resolution are probably responsible for the underestimation of temperatures >30 ˝C (Figure 5).However, the RMSE is <2 ˝C, thus, simulated temperatures are clearly closer to the measured canopy temperatures than air temperatures, especially for the W0 treatment (Figure 2).Note that the number of observation days where air temperatures exceed 30 ˝C is limited.

Target
Minimum daily canopy temperatures, mostly representing night or early morning temperatures (data not shown), are simulated well (RMSE ~1 ˝C) using minimum air temperatures and by using crop height (for phase I), or VPD and the Eact/ETP ratio (for phase II) as additional explanatory variables (Table 4).During plant growth (phase I), the crop height is a controlling factor of the night-time temperature profiles, and, consequently, affects the radiometric temperatures measured by the infrared thermometers.As found for T c,mean , for subsequent developmental stages, components of the canopy water balance help explain the variability of T c,min .At night, transpiration is low or close to zero, thus, the actual evaporation is a suitable model covariate.Night-time or early-morning evaporation is driven by the evaporative demand of the air (VPD), giving reason for using VPD as a model covariate.
Note, that daily mean, maximum and minimum air temperatures were derived from hourly diurnal temperature measurements.The accuracy of models using temperatures characterized by higher uncertainty needs to be investigated in future studies.2-4).
Calculations of the relative importance of the covariates for the total estimated R 2 of the presented empirical models (Figure 6) highlight the dependence of their contribution on the crop water status: while, for the fully irrigated plots (W2), the air temperature contributes from ~60%-70% up to nearly 100% to the estimated R 2 of the multiple linear regression (MLR) model, the contribution decreases to <50% for Tc,mean and Tc,max of the W0 plots.In consequence, the percentage contribution of, e.g., incoming radiation and the VPD to the R 2 of the W0 canopy temperature models increases.Consequently, for simulating canopy temperatures under the influence of water deficits, solemnly using air temperature does not provide a suitable proxy for temperatures at the canopy level, and interactions between environmental and meteorological variables need to be incorporated.
A widely applied strategy to model temperature conditions within a crop canopy is the use of surface energy balance models [22,23].Since such models are optimized to simulate the partitioning of available energy in latent and sensible heat fluxes, derived canopy temperatures are unrealistic if underlying assumptions, such as the closure of the energy balance, are not fulfilled.They are computationally expensive and require data on land surface physical parameters and aerodynamic properties of the canopy.Anomalies of these variables with respect to "average conditions", which are common for water-limited conditions, limit the use of standard literature values.Consequently, the use of derived canopy temperatures in crop system models is limited.Contrastingly, suggested empirical models are a strong simplification of the complex processes related to the crop canopy temperature dynamics [2,4,17] and do not predict canopy temperatures over different scales in time  2-4).
Calculations of the relative importance of the covariates for the total estimated R 2 of the presented empirical models (Figure 6) highlight the dependence of their contribution on the crop water status: while, for the fully irrigated plots (W2), the air temperature contributes from ~60%-70% up to nearly 100% to the estimated R 2 of the multiple linear regression (MLR) model, the contribution decreases to <50% for T c,mean and T c,max of the W0 plots.In consequence, the percentage contribution of, e.g., incoming radiation and the VPD to the R 2 of the W0 canopy temperature models increases.Consequently, for simulating canopy temperatures under the influence of water deficits, solemnly using air temperature does not provide a suitable proxy for temperatures at the canopy level, and interactions between environmental and meteorological variables need to be incorporated.
A widely applied strategy to model temperature conditions within a crop canopy is the use of surface energy balance models [22,23].Since such models are optimized to simulate the partitioning of available energy in latent and sensible heat fluxes, derived canopy temperatures are unrealistic if underlying assumptions, such as the closure of the energy balance, are not fulfilled.They are computationally expensive and require data on land surface physical parameters and aerodynamic properties of the canopy.Anomalies of these variables with respect to "average conditions", which are common for water-limited conditions, limit the use of standard literature values.Consequently, the use of derived canopy temperatures in crop system models is limited.Contrastingly, suggested empirical models are a strong simplification of the complex processes related to the crop canopy temperature dynamics [2,4,17] and do not predict canopy temperatures over different scales in time and space.However, the aim of our study was to provide a practical tool in plant growth models for reliable canopy temperature predictions using a small number of input variables.Presented models were designed for simulating the winter wheat canopy temperature during stem extension, heading and ripening growth stages.For transferring them to other crops and/or observation scales, equations can be refitted using local radiometric temperature observations, micrometeorological measurements and the output of a day step model.Coupling such empirical models to existing crop system models allows for including realistic near surface and canopy temperatures for site-and scale-specific estimates of crop development and the impact of heat stress on crop yield, provided that these models reliably simulate the components of the soil water balance.We are aware that current agroecosystem models were designed to use air temperatures as input for simulating plant development and the impact of heat stress.This may imply the need for the recalibration of some parameters of temperature-dependent processes.However, the simulation of canopy temperatures is an important step towards a more realistic representation of temperature stress effects and the interaction between heat and drought stress in models [15], thus, for simulating the effect of temperatures as a limiting, defining and reducing factor.and space.However, the aim of our study was to provide a practical tool in plant growth models for reliable canopy temperature predictions using a small number of input variables.Presented models were designed for simulating the winter wheat canopy temperature during stem extension, heading and ripening growth stages.For transferring them to other crops and/or observation scales, equations can be refitted using local radiometric temperature observations, micrometeorological measurements and the output of a day step model.Coupling such empirical models to existing crop system models allows for including realistic near surface and canopy temperatures for site-and scale-specific estimates of crop development and the impact of heat stress on crop yield, provided that these models reliably simulate the components of the soil water balance.We are aware that current agroecosystem models were designed to use air temperatures as input for simulating plant development and the impact of heat stress.This may imply the need for the recalibration of some parameters of temperature-dependent processes.However, the simulation of canopy temperatures is an important step towards a more realistic representation of temperature stress effects and the interaction between heat and drought stress in models [15], thus, for simulating the effect of temperatures as a limiting, defining and reducing factor.

Case Study: Wheat Canopy versus Air Temperature
The shape of the cumulative sum curves of the difference between daily averaged canopy temperature and the air temperature (∆T) is highly dependent on the observation year (Figure 7).However, cumulative ∆T sums of the drought treatment plot (W0) almost constantly increase throughout the growing season, i.e., in 2014 (Figure 7).Contrastingly, we calculated negative cumulative sums of ∆T for the well irrigated plot (W2) at the end of the growing season (except in 2010).In 2011 and 2013, the canopy of the irrigated plot is almost constantly cooler compared to the air temperature, whereas, in 2010 and 2014, it was warmer until mid/end of June and cooler during the subsequent growing stages.Our data indicate increasing negative temperature anomalies for fully irrigated and positive anomalies for not irrigated wheat plots during the period ranging from

Case Study: Wheat Canopy versus Air Temperature
The shape of the cumulative sum curves of the difference between daily averaged canopy temperature and the air temperature (∆T) is highly dependent on the observation year (Figure 7).However, cumulative ∆T sums of the drought treatment plot (W0) almost constantly increase throughout the growing season, i.e., in 2014 (Figure 7).Contrastingly, we calculated negative cumulative sums of ∆T for the well irrigated plot (W2) at the end of the growing season (except in 2010).In 2011 and 2013, the canopy of the irrigated plot is almost constantly cooler compared to the air temperature, whereas, in 2010 and 2014, it was warmer until mid/end of June and cooler during the subsequent growing stages.Our data indicate increasing negative temperature anomalies for fully irrigated and positive anomalies for not irrigated wheat plots during the period ranging from the flowering stage until maturity.Obviously, during these periods, the higher level of canopy transpiration of fully irrigated crops can significantly reduce the effect of short-term heat events on plant physiological processes.Although modeled and measured absolute canopy temperatures from our study might strongly differ from those data observed at other crops, sites, and observation years, we assume that such differences between air and canopy temperatures of irrigated and non-irrigated plots are critical for the results of crop system models.the flowering stage until maturity.Obviously, during these periods, the higher level of canopy transpiration of fully irrigated crops can significantly reduce the effect of short-term heat events on plant physiological processes.Although modeled and measured absolute canopy temperatures from our study might strongly differ from those data observed at other crops, sites, and observation years, we assume that such differences between air and canopy temperatures of irrigated and non-irrigated plots are critical for the results of crop system models.If using maximum air temperatures (Tair,max), the number of days exceeding temperature thresholds of 20 °C and 25 °C is close to that of using maximum canopy temperatures from the fully irrigated treatment plot (Table 6).However, the number of days is significantly higher if using Tc,max W0.Note that measured temperatures for the drought treatment plot from 2014 are only available until beginning of May and that there is a data gap of seven days in 2011 (Table 1).Despite large interannual differences in the number of days exceeding sample temperature thresholds, our data highlight the need to account for air-canopy temperature differences in threshold-based modeling approaches.
To highlight the impact of temperature differences on summed thermal units, which are frequently used in models to simulate crop development, we computed ETU sums above the optimal temperature of 20 °C [52].Patterns of the sum curves clearly indicate lowest values if using Tc,mean W2 and highest values using Tc,mean W0 (Figure 8).Note that starting dates differ for each observation year.Absolute summed differences increase during the wheat flowering and ripening stages.It is well known that the difference between canopy and air temperatures increases at higher levels of crop water stress.Our results, however, suggest significant discrepancies in Tair,mean-and Tc,mean-based ETU calculations not only for the drought treatment but also for the well irrigated plot.These If using maximum air temperatures (T air,max ), the number of days exceeding temperature thresholds of 20 ˝C and 25 ˝C is close to that of using maximum canopy temperatures from the fully irrigated treatment plot (Table 6).However, the number of days is significantly higher if using T c,max W0.Note that measured temperatures for the drought treatment plot from 2014 are only available until beginning of May and that there is a data gap of seven days in 2011 (Table 1).Despite large interannual differences in the number of days exceeding sample temperature thresholds, our data highlight the need to account for air-canopy temperature differences in threshold-based modeling approaches.
To highlight the impact of temperature differences on summed thermal units, which are frequently used in models to simulate crop development, we computed ETU sums above the optimal temperature of 20 ˝C [52].Patterns of the sum curves clearly indicate lowest values if using T c,mean W2 and highest Crops show threshold responses to meteorological and environmental conditions, especially to temperature [53,54].Current crop system models use the sum of daily temperatures exceeding a base temperature and crop-specific temperature thresholds (cardinal temperatures) for deriving factors and/or linear and non-linear response functions controlling the simulation of phenology, growth (biomass accumulation), biochemical processes (nutrient uptake) and the harvest index.Well defined threshold responses are generally derived from growth chamber experiments performed under controlled conditions.We suggest avoiding the use of air temperatures recorded at 1.5 or 2 m height above the canopy for threshold-based modeling approaches since large positive or negative temperature anomalies to the near-surface temperatures can be expected.Such differences can have important consequences for accuracy of plant physiological processes and crop yield simulations.
Currently, few models are able to account for increasing canopy temperatures during drought stress by offering optional corrective temperature mechanisms based on water stress indices [12,55].However, next to the effects of drought stress on canopy temperatures our study further highlights the importance of considering negative temperature anomalies (transpiration cooling) found for fully irrigated crops, when simulating crop phenological development and the impact of heat stress effects on crop yields.

Conclusions
Canopy temperatures of water-deficit and fully irrigated wheat plots significantly differ from air temperature measurements.Such temperature anomalies need to be considered in crop system models.Variables available from daily standard meteorological observations (vapor pressure deficit, minimum, maximum and mean temperature, incoming radiation) and from widely applied agro-ecosystem models (potential and actual transpiration, evaporation, leaf area index, and crop height) allow for approximating wheat canopy temperatures using simple empirical linear regression models.Significant differences in model coefficients for the pre-and the post-heading phase can be related to the stronger influence of the components of the crop and the soil water balance during the post-heading stage.The use of canopy temperatures as an input for crop system models will lead to the need for adapting temperature thresholds and the temperature sum calculations.However, canopy temperatures provide a more realistic measure of microclimatic conditions and their potentially amplifying or mitigating effect on plant response to temperature stress.

Figure 1 .
Figure 1.Density plot of available surface temperature observations (Tc, daily averages of hourly values) for the different observation years and irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation).

Figure 1 .
Figure 1.Density plot of available surface temperature observations (Tc, daily averages of hourly values) for the different observation years and irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation).

Figure 2 .
Figure 2. Scatter plot of mean, maximum and minimum canopy temperatures (Tc) and mean, maximum and minimum air temperatures (Tair) derived from the aggregation of hourly measurement data to daily values for three different irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation, r = Pearson's correlation coefficient).The dashed black line indicates the line of identity.

Figure 3 .
Figure 3. Density plot of the mean daily difference between crop and air temperatures (∆Tmean) for the different observation years and irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation).

Figure 2 .
Figure 2. Scatter plot of mean, maximum and minimum canopy temperatures (T c ) and mean, maximum and minimum air temperatures (T air ) derived from the aggregation of hourly measurement data to daily values for three different irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation, r = Pearson's correlation coefficient).The dashed black line indicates the line of identity.

Figure 2 .
Figure 2. Scatter plot of mean, maximum and minimum canopy temperatures (Tc) and mean, maximum and minimum air temperatures (Tair) derived from the aggregation of hourly measurement data to daily values for three different irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation, r = Pearson's correlation coefficient).The dashed black line indicates the line of identity.

Figure 3 .
Figure 3. Density plot of the mean daily difference between crop and air temperatures (∆Tmean) for the different observation years and irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation).

Figure 3 .
Figure 3. Density plot of the mean daily difference between crop and air temperatures (∆T mean ) for the different observation years and irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation).

Figure 4 .
Figure 4. Scatter plot of the daily range of canopy temperatures (Tc, hourly values) and air temperatures (Tair, hourly values) for the complete observation period and for all irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation).

Figure 4 .
Figure 4. Scatter plot of the daily range of canopy temperatures (T c , hourly values) and air temperatures (T air , hourly values) for the complete observation period and for all irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation).

Figure 5 .
Figure 5. Scatter plot of the measured and the simulated maximum daily canopy temperatures (Tc,max) for the training and the testing data set during the pre-heading (I) and the post-heading phase (II) for all irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation).The dashed gray line shows the line of identity.Equations of the multiple linear regression models and r 2 (= coefficient of determination) are given (abbreviations as for Tables2-4).

Figure 5 .
Figure 5. Scatter plot of the measured and the simulated maximum daily canopy temperatures (T c,max ) for the training and the testing data set during the pre-heading (I) and the post-heading phase (II) for all irrigation treatments (W0 = no irrigation, W1 = medium deficit irrigation and W2 = full irrigation).The dashed gray line shows the line of identity.Equations of the multiple linear regression models and r 2 (= coefficient of determination) are given (abbreviations as for Tables2-4).

Figure 6 .
Figure 6.Estimated relative importance of covariates for simulating Tc,mean, Tc,max and Tc,min for the training data set during the pre-heading (I) and the post-heading phase (II).Abbreviations as for Tables 2-4.

Figure 6 .
Figure 6.Estimated relative importance of covariates for simulating T c,mean , T c,max and T c,min for the training data set during the pre-heading (I) and the post-heading phase (II).Abbreviations as for Tables 2-4.

Figure 7 .
Figure 7. Cumulative sums of canopy to air temperature differences (∆T ( ˝C)) for the drought treatment plot (dotted line) and the irrigated plot (dashed line) using mean daily modeled surface temperatures (∆T = T c,mean -T air,mean ) in 2010 (top left), 2011 (top right), 2013 (bottom left), 2014 (bottom right) (Hohenschulen research site).

Table 1 .
Summary of the canopy temperature data set (HS = Hohenschulen research site, BS =

Table 4 .
Statistics of the multiple linear regression model developed using a forward stepwise-based selection procedure and quantile regression analyses for predicting T c,min (= minimum daily canopy temperature) of a winter wheat canopy; SE = standard error, T air,min = minimum daily air temperature ( ˝C), CH = crop height (m), VPD = vapor pressure deficit (hPa), Eact/ETP = ratio of actual evaporation (Eact (mm¨d ´1)) to potential evapotranspiration (ETP (mm¨d ´1)).