Dynamic leaf energy balance: deriving stomatal conductance from thermal imaging in a dynamic environment

Spatiotemporal co-ordination of the stomatal response to light intensity over the lamina of different wheat leaves was assessed using a novel dynamic energy balance model.


Introduction
After a period of increased crop production over the past 50 years (Wik et al., 2008;Pingali, 2012;Tillman et al., 2015), increases in yield have fallen to ~1% per annum (Fischer and Edmeades, 2010;Ramankutty et al., 2018). To meet the increasing food demand, crop yield needs to increase at a rate of 2.4% per annum over the next few decades (Tilman et al., 2011;Ray et al., 2012Ray et al., , 2013Long et al., 2015). The expanding gap between crop production and demand will drive more acute food insecurity across the globe that will only intensify with global population growth (Long et al., 2015;FAO et al., 2018). Along with new agronomic practices, the development of innovative breeding techniques and novel trait discoveries are required to improve yield. Many of the current commercial wheat cultivars have been bred with yield as the main target trait, an approach that is now reaching its limit in term of potential improvement (Fischer and Rebetzke, 2018). Breeding for yield has resulted in a loss of diversity in wheat (Haudry et al., 2007) and in the indirect selection for traits that contribute to higher yields (Reynolds et al., 1999;Richards, 2000;Fischer and Rebetzke, 2018). The natural genetic diversity between cultivars or landraces provides the opportunity to discover desirable traits (e.g. resistance to pest or stress; high water use efficiency) that, when crossed with high-yielding varieties, could produce progenies with improved performance and yield under a given environment. However, identifying individuals with the desired eco-physiological and agronomic responses and traits requires the development of appropriate phenotyping tools (Furbank and Tester, 2011;Reynolds et al., 2018).
Plant phenotyping to characterize diversity in physiological traits often requires measuring large numbers of plants ideally under controlled (e.g. glasshouse) and sometimes under noncontrolled environmental (e.g. field) conditions to identify those plants with specific allelic combinations that result in the phenotypes with increased agronomic value (Fiorani and Schurr, 2013;Fahlgren et al., 2015). In the last decade, significant advances have been made in genomics resources and tools that have greatly increased our capacity to understand and detect genetic diversity. However, phenotyping approaches have not advanced at the same pace and have been reported in the literature as a major technical bottleneck (Furbank and Tester, 2011;Fiorani and Schurr, 2013;Fahlgren et al., 2015). High-throughput phenotyping platforms have been developed to address some of this shortfall, but can still be limited in the type and rapidity of the measurement that can be undertaken (Li et al., 2014;Rahaman et al., 2015). Further improvements in plant phenotyping methods to match the progress made in genetics are required to identify new traits increasing crop yield in the near future (Furbank and Tester, 2011;McCouch, 2013;Prashar and Jones, 2014).
High-throughput phenotyping of leaf traits often exploits imaging techniques that permit rapid and remote measurements of plants, to assess differences in development, morphological features, pigment concentration, and responses to the environment. Phenotyping platforms typically use a combination of monochromatic, RGB, thermal, and multispectral cameras providing signals that are often an indirect measure of the physiological plant response/traits (e.g. chlorophyll fluorescence as an indicator of the photosynthesis rate). In most cases, interpretation of the raw measurement signal requires extra steps to convert the raw data into biologically significant findings or responses (Fahlgren et al., 2015;Ghanem et al., 2015). For example, leaf temperature measured using thermography requires leaf energy balance equations to derive leaf transpiration (E) and stomatal conductance to water vapour (g sw ) (Jones, 2013). Leaf stomatal conductance is often considered an important trait for future yield improvements and is of great interested to physiologists, as stomatal behaviour influences traits such as photosynthetic CO 2 uptake, water loss, and leaf temperature (Reynolds et al., 1999;Richards, 2000;Rebetzke et al., 2001;Munns et al., 2010;Fischer and Rebetzke, 2018), all of which impact crop performance. Recent studies examining stomatal behaviour have highlighted that slow stomatal responses to changing environmental stimuli can limit photosynthesis rates (A) and result in unnecessary water losses, decreasing productivity and plant water use efficiency (Lawson et al., 2010;Lawson and Blatt, 2014;McAusland et al., 2016;Matthews et al., 2018). Reduced stomatal limitation of photosynthesis under a dynamic environment has been shown to impact plant growth and biomass significantly Way and Pearcy, 2012;McAusland et al., 2016;Vialet-Chabrand et al., 2017b), and therefore could provide a novel breeding target to improve yield (Lawson and Blatt, 2014). However, currently there are limited phenotyping tools capable of exploring variation in stomatal behaviour and kinetic responses (McAusland et al., 2013). Although leaf temperature and thermography have been used successfully in the past to identify mutants with altered stomatal apertures in response to drought, air relative humidity, [CO 2 ], or light intensity (Merlot et al., 2002;Wang et al., 2004;Costa et al., 2015;Takahashi et al., 2015), these studies have focused on steady measurements of leaf temperature. Deriving g sw from leaf temperature kinetics under a dynamic environment is a challenge that has currently not been addressed, but is crucial to understand time-integrated variations of plant carbon and water budget.
Despite the excellent reports on the practical use of infrared thermometry to estimate E and g sw (Jones, 1999;Leinonen et al., 2006;Guilioni et al., 2008;Maes et al., 2016;Jones et al., 2018), these are limited to steady-state assumptions that are difficult to support under a dynamic environment (Jones and Leinonen, 2003). Deriving g sw from thermography is generally performed by applying a reworked version of the Penman-Monteith equation describing the sum of radiative energy received or lost by the leaf, as well as mass transfer processes to the atmosphere (Monteith and Unsworth, 2012). Simplifications have been proposed in the literature that involve using reference material (dry/wet leaf replica with similar optical properties) to eliminate the need for longwave radiation and/or humidity measures (Leinonen et al., 2006;Guilioni et al., 2008). Although energy balance equations are a mechanistic description of the processes involved in changes of leaf temperature, some energy fluxes (e.g. longwave radiation from the surrounding environment) and properties of the material studied (e.g. absorbance and emissivity) are difficult to evaluate, particularly in the confined space of a phenotyping platform. A significant step toward the use of thermal imaging to derive stomatal behaviour would be to redefine the Penman-Monteith combination equation for use under dynamic environmental conditions and associate it with a statistical approach to estimate the parameters that cannot be assessed without specific laboratory equipment.
To simplify the usage of energy balance equations, the difference in temperature between a transpiring and non-transpiring leaf can be used to calculate the amount of energy lost by transpiration (latent heat). Several methods have been investigated to measure the temperature of a non-transpiring leaf using dry reference materials that replicate the leaf optical properties (green paper or fabrics) or application of grease to the leaf surface to prevent transpiration (Leinonen et al., 2006;Guilioni et al., 2008). Although these methods appear as an interesting option, there is no guarantee that the reference material has the same optical (e.g. absorbance and reflectance) and thermal properties (e.g. specific heat capacity and emissivity) as the intact leaf. This means that under changing environmental conditions, the leaf and reference material temperatures could behave differently, introducing bias and error in estimating g sw from thermographs (Jones and Leinonen, 2003;Jones, 2013). Instead of attempting to mimic leaf properties, a more robust approach would be to directly include differences in thermal and optical properties in the energy balance model and predict leaf thermal kinetics from a reference material. The energy balance model allows an estimation of energy loss by transpiration, which is dependent on the gradient of water vapour from the leaf to the atmosphere and the total conductance to water vapour (g tw , the ease with which water vapour diffuses through stomatal pores and the boundary layer). As g tw is the combination of two conductances in series, g sw and the boundary layer conductance (g bw ), it is possible to derive g sw from g tw if g bw is known. A common method used to estimate g bw consists of applying a cycle of heating/cooling to a leaf replica (similar shape) resulting in an exponential increase or decrease in temperature from which g bw can be derived using the value of the slope (see details in Jones, 2013). Although this method provides an estimate of the boundary layer, it can only provide a value after several minutes of response when the environmental conditions are relatively stable (Leuning et al., 1989;Leuning and Foster, 1990;Brenner and Jarvis, 1995;Stokes et al., 2006;Katsoulas et al., 2007). Here, we propose the use of a passive method relying only on the difference in energy balance observed under a dynamic environment of references with different thermal properties to estimate the boundary layer conductance.
Rearranging and solving the Penman-Monteith combination equation for g sw supposes that the energy budget is closed (i.e. that all the fluxes are taken into consideration in the thermal budget), and integrates the error on each thermal flux in the calculation of g sw . An alternative approach is to predict leaf temperature using an energy balance model with an inbuilt dynamic model of g sw , which, once fitted on the observed temperature, will provide the temporal kinetics of g sw , as well as key parameter values related to the leaf thermal response. Using Bayesian inference, the probability of each parameter value included in the energy balance model and the prediction errors can be characterized for different individuals, enabling a more precise quantification of the trait diversity associated with the parameters.
The aim of this study was to provide a method to process thermograms describing leaf temperature kinetics under fluctuating conditions and derive temporal responses of stomatal conductance (g sw ), using a new interpretation of energy balance equations. The results deliver important information on the diversity of the leaf responses to changes in light intensities, providing a breakthrough in data processing for phenotyping. Model performance and its application for plant phenotyping are presented, as well as a description of the biological importance of the response traits estimated.

Materials and methods
The Penman-Monteith equation (Penman, 1948;Monteith, 1965) combines the surface energy balance with mass transfer (the transport of water vapour from the leaf to the atmosphere) and was originally developed to compute transpiration from cropped surfaces using commonly measured weather data (solar radiation, air temperature, vapour content, and wind speed). Under the assumption of steady-state environmental conditions, the Penman-Monteith combination equation can be solved for leaf stomatal conductance (Monteith and Unsworth, 2012;Jones, 2013) and is therefore widely used in association with thermography to study plant response. Here, we further develop this equation by including temporal effects that affect thermal exchange between the leaf and its environment using differential equations. Although such equations are more complex to use and solve, they are essential to understand the dynamics of leaf temperature and water loss under a changing environment.

Energy balance model under steady-state environmental conditions
Under steady-state environmental conditions, leaf temperature is at equilibrium and therefore the sum of the radiative energy received and lost by the leaf is null ( Fig. 1; see Table 1 for symbol description and units). Any imbalance in the energy fluxes due to change in environmental conditions impacts the energy stored by the leaf, and results in variation in leaf temperature. Neglecting any metabolic storage (heat stored as chemical bond energy, e.g. photosynthesis), the energy balance equation reduces to: where R n is the net radiation (W m -2 ), C is the sensible heat transfer (W m -2 ), λ is the latent heat of evaporation of water (J kg -1 ), E is the evaporative flux (kg m -2 s -1 ), and S is the net physical storage (W m -2 , causing a change in leaf temperature and equal to 0 under steady-state conditions). After a change in environmental conditions, a delay is required to reach a new equilibrium depending on the material utilized, and therefore Equation 1 may over-or underestimate leaf temperature during the transition from one equilibrium to another. T air , air temperature; I (dir) , direct irradiance; I (dif) , diffuse irradiance; α l , leaf absorbance to shortwave irradiance; α r , reference absorbance to shortwave irradiance; g bh , boundary layer conductance to heat transfer; ρ, air density; C s , specific heat capacity of humid air; E, leaf transpiration. Factor 2 in the equations represent the input or output energy for the two sides of the leaf/reference.

Energy balance model under dynamic environmental conditions
Under dynamic environmental conditions, temperature kinetics of an object depend on the thermal and optical properties of the object, which describes how efficiently the energy is captured, exchanged, and stored. A differential equation describing the time dependency of temperature kinetics is available from the literature (Monteith and Unsworth, 2012;Jones, 2013): where ρ* and C* p are the density (kg m -3 ) and specific heat capacity ( J kg -1 K -1 ), respectively, of leaf tissue and l * is the leaf thickness (m). Steady-state energy balance equations are generally simplified by using the temperature of a reference material mimicking a leaf without transpiration, as described by Jones and colleagues (Leinonen et al., 2006;Guilioni et al., 2008;Jones et al., 2018). The difference in energy balance between the leaf and the reference material (placed in the same orientation as the leaf) allows a simplification of the energy fluxes (see below) to remove most of the effects due to the surrounding environment (e.g. longwave radiation). Including the thermal properties of the material (ρ*, C* p , l *) in the energy balance model can improve the accuracy of the model but also results in different denominators complicating the equation. Solving the differential equation that describes the difference in energy balance between the leaf and a reference material using Equation 2 can enable the use of any material with known properties as a reference. An advantage of this proposed approach is that instead of trying to mimic the properties of the leaf, the differences in thermal and optical properties between the two objects are directly included in the equations. It is possible to predict the temperature (T 1 , °K) of any material (e.g. leaf) using the temperature kinetics of a reference (T 2 , °K; e.g. black painted aluminium sheet) with knowledge of the optical and thermal properties of the material. Under a given environment, the temperature kinetics of an object are predicted using the difference in energy balance between the object and the reference, which is described by the following equations: where ρ* 1 and ρ* 2 are the density (kg m -3 ), C* p1 and C* p2 are the specific heat capacity (J kg -1 ), and l* p1 and l* p2 are the thickness (m) of the objects measured. Therefore, k 1 and k 2 represent the amount of energy per unit area required to change the temperature of the material by 1 °K (J m -2 K -1 ). It is important to note that a factor of 2 can be applied to the latent heat term (λE 1/2 ) if the object is transpiring from both surfaces. Rearranging the equation by moving the denominator to the left-hand side results in: This configuration allows further simplification by regrouping similar energy fluxes: If E 2 is for a non-transpiring reference material (e.g. an aluminium replica) and therefore has a value of 0, the equation would be: The net radiation (R n ) for a horizontal two-sided object is defined by its longwave and shortwave radiation exchanges with its environment ( Fig. 1): where subscript x is either 1 or 2, α x is the absorbance to shortwave radiations (incident and diffuse; I s , W m -2 ), θ the Stefan-Boltzmann constant (W m -2 K -4 ), ε x the emissivity, T x the temperature of the object studied (°K), and L d the longwave radiation (W m -2 ) received from the surrounding environment (e.g. soil, wall). The factor of 2 in 2θε x T x 4 correspond to the longwave radiation energy lost by each side of the object. Using a reference material placed under similar conditions to the object measured removes the need to quantify L d as it disappears from the equation when the difference in energy flux is calculated (see below). This assumption is valid when the reference material and the leaf are surrounded by a similar thermal environment, with an increasing risk of error when the distance between the two objects increases (e.g. large imaging area).
The sensible heat is the energy lost by conduction and convection and is defined as: where g bhx is the one-sided boundary layer conductance to heat transfer (m s -1 ), ρ the air density (kg m -3 ), C s the specific heat capacity of humid air (J kg -1 K -1 ), and T a is the air temperature (°K). The factor of 2 in the equation correspond to the energy lost by each side of the object. Air properties such as ρ and C s are described in Supplementary Table S1 at JXB online. Substituting Equations 8 and 9 in Equation 7 gives: Details required to calculate the different terms are provided in Supplementary Algorithm S1. Temperature kinetics can be predicted using the following generic equation: where: k 2 dT2 dt is the change in energy of the reference material per unit of time; I s (α 2 -α 1 ) the difference in absorbed shortwave radiations; 2θ(ε 2 T 2 4 ε 1 T 1 4 ) is the difference in energy lost by longwave radiations; 2ρC s [g bh2 (T 2 -T air )-g bh1 (T 1 -T air )] is the difference in energy exchange by conduction and convection with the atmosphere; and λE 1 is the energy lost by transpiration.
Such an equation can be used to calculate T 1 at any given time using an ordinary differential equation (ODE) solver. The solver required Equation 11 to be calculable at any time t supposing that the derivative of the reference temperature kinetics (dT 2 /dt) and the discrete estimates of the environmental variables [I s , relative humidity (RH), T air , P atm ] can be interpolated as a function of time. Using a smooth spline on recorded values of T 2 and each environmental variable removes high-frequency noise (making the ODE solver more efficient) and provides functions that for any time t return a value and its derivative ( Supplementary Fig. S1). Using an ODE solver requires an initial value for the predicted temperature that was included as an estimated parameter in the Bayesian inference described below. Equation 11 is a 'stiff ' equation requiring a solver with variable step size (e.g. 'cvode' at https://computation.llnl.gov/projects/sundials), as standard methods such as Runge-kutta were orders of magnitude slower to solve this equation.

Leaf transpiration model
In response to changes in leaf energy balance, leaf transpiration varies, regulating leaf temperature to, for example, maintain an optimal range for photosynthesis. Leaf transpiration (E) can be calculated using Equation 12: Patm a conversion factor from Pa to kg m -3 and RTleaf Patm a conversion factor from mol m -2 s -1 to m s -1 , and where P atm is the atmospheric pressure (Pa), T leaf the leaf temperature (°K), g tw the total conductance to water vapour (mol m -2 s -1 ), e s the leaf internal vapour pressure, and e a the air vapour pressure (Pa).
Total conductance to water vapour (g tw , mol m -2 s -1 ) is the sum in series of the boundary layer conductance to water vapour (g bw , mol m -2 s -1 ) and stomatal conductance to water vapour (g sw , mol m -2 s -1 ): where g bw depends on leaf morphology and wind speed, and g sw depends mainly on the number of stomata and change in aperture in response to the surrounding environment. Substituting Equations 12 and 13 in Equation 11 allows us theoretically to solve Equation 11 for g sw , assuming that the different derivatives included in Equations 11 can be extrapolated using splines adjusted on the temperature kinetics. Although estimating the derivative at different time points using a spline is an approach that could be considered, it should be done with caution as it can give unstable g sw values due to the sensitivity to measurement noise of the method. Indeed, errors in estimating the derivatives of temperature kinetics from noisy data will propagate in the estimation of g sw . In the case of the temperature kinetics of the reference material, this error is minimized, as the pattern of variation (exponential increase) is simple to approximate, which is not always the case for the leaf temperature kinetics. Therefore, the approach chosen here was to directly include a model to predict g sw in the energy balance equations to calculate the derivative and solve it using an ODE solver, leaving the errors propagated to the predicted temperatures. The prediction errors of the model can thus be taken into consideration when fitting to the observed data.

Stomatal conductance model
The temporal response of stomatal conductance to water vapour (g sw , mol m -2 s -1 ) can be predicted using a previously described dynamic model (Vialet-Chabrand et al., 2017b). In the case of a continuous variation in light intensity, differential equations are required to model g sw (Vialet-Chabrand et al., 2013. However, for a step change in light intensity, such an equation can be solved analytically, resulting in a sigmoidal equation, which is simpler and faster to calculate as a function of time (t, s): Where g 0 is the initial value of g sw at t 0 =0, G is the steady-state target of g sw , τ the time constant of the response (s), ϕ the initial time lag (s), and s l the slope of the slow decrease/increase in g sw (mol m -2 s -2 ). The term Å ã has been included to ensure a proper scaling of the sigmoidal equation between g 0 and G even when ϕ is positive. Without this correction at t=0, the term e −e φ τ+1 would be positive when ϕ was positive, introducing an offset in the estimation of g 0 and changes in the scale of the response. Plants can biologically control the energy lost by leaf transpiration by dynamically adjusting g sw or by adjusting the boundary layer conductance by changing the leaf shape (e.g. curling) over a longer period of time.

Boundary layer conductance
Using two non-transpiring references mimicking leaf shape and having identical thermal properties (same material), but different optical properties (e.g. black and white) enables Equation 9 to be simplified to determine the boundary layer conductance: where T w and T b are the temperature (°K) of the white and the black references, respectively. Knowing the thermal (k, ε w , ε b ) and optical (α w , α b ) properties of the references, it is possible to derive g bh (m s -1 ) from Equation 15 by adjusting its value to fit the predicted white reference temperature kinetics to the measured data. Depending on the experimental conditions, it is possible to assume that g bh is stable over the measurement period and therefore only one value needs to be estimated. For field applications, it is also possible to include a model predicting the boundary layer conductance from wind speed measurements and the dimensions of the references (Jones, 2013). This enables us to take into consideration the effect of wind speed on leaf energy balance under a dynamic field environment.

Predicting leaf temperature kinetics and stomatal conductance to water vapour
Model implementation Stomatal conductance kinetics were predicted by fitting a model combining Equations 9-13 on leaf temperature kinetics under a dynamic environment. Using environmental variables as inputs [photosynthetic photon flux density (PPFD), µmol m -2 s -1 ; RH, %; T air , °C; P atm , Pa], the model predicted stomatal conductance kinetics (Equation 12) in response to changes in light intensity, which, combined with estimation of leaf boundary layer conductance (Equation 13), was used to calculate the transpiration rate (Equations 10-11). The latent heat of evaporation was then calculated and included in the leaf energy balance (Equation 9) to predict leaf temperature variations. Equation 9 was applied twice using two different references (black and white) to calculate leaf temperature kinetics, constraining the estimation of g bh as in Equation 13. The two different predictions of leaf temperature were compared simultaneously with the observed lead temperature using Bayesian inference (see Supplementary Algorithm S1).

Bayesian inference
Bayesian inference was used to tune parameter values to fit model predictions to observed data. Estimating parameter values of such a model is challenging because each parameter has a range of values that can produce equivalent model outputs depending on the precision of observations and the interaction with other parameters. Using Stan, a statistical modelling platform (http://mc-stan.org/), intervals representing 95% of the probable parameter values (credible interval) were estimated for each parameter. The credible interval included the correlation between parameter and the error due to the noisy data. The model outputs were fitted on the observed data by exploring the parameter space using three Markov chain Monte Carlo (MCMC) with 500 iterations (250 iterations without adaption). All chains converged without divergent transition to the same parameter values (Rhat <1.1), with effective sample size values >200 (Carpenter et al., 2017). Further information on the use of Bayesian inference to find parameter values of ordinary differential equations can be found online at http://mc-stan.org/events/stancon2017-notebooks/ stancon2017-margossian-gillespie-ode.html. Originally, Stan did not include the possibility to use spline functions in the Bayesian model. The function was included as an external C++ function wrapping the spline function available from the Boost Library into the template system used by Stan.
A possible problem during the model fitting could come from comparing observed leaf temperature kinetics with predicted temperatures using two different references, which could result in bimodal parameter distributions. However, the convergence toward the observed temperature of both sets of predicted temperatures is mainly dependent on the value of the boundary layer conductance (g bh , Equation 13), which is adjusted during the Bayesian inference, constraining the value of g bh . Measurements of g bh at different times of the day did not differ significantly and a constant g bh was assumed for both references during the experiment to predict leaf gas exchange. Moreover, the conditions that generally influence the boundary layer conductance (e.g. air mixing provided by the fan) were relatively constant within our imaging area.

Model validation using a leaf replica with constant g sw
Under similar conditions to those experienced by plants under step changes in light intensity, accuracy and precision of the model outputs (leaf temperature and conductance) were tested using a reference built to mimic leaf transpiration with a constant conductance to water vapour (Fig. 2). The aim was to compare the conductance predicted by the model, measured using gas exchange and calculated using a physical diffusion equation (Lehmann and Or, 2015). The replica was located at the same position and with the same orientation (horizontal) as the leaves and therefore received the same light intensity (430 μmol m -2 s -1 ) as the measured leaves, whilst the aluminium references (black and white) were at a lower location and to one side of the image, and therefore received a slightly lower light intensity of 300 μmol m -2 s -1 . The replica was an aluminium plate covered on one side by black electrical tape with known absorbance and emissivity, and on the other side with a felt fabric enclosed in a plastic microporous sheet. The microporous sheet had pores of a known and standard diameter of 0.5 mm and 40 µm depth, arranged in a grid of 3 mm pitch (density: 160 pores inch -2 ). When the felt was saturated with water, transpiration was dictated by the size and density of the pores, and the variation in environmental conditions (e.g. RH and wind speed). Conductance was calculated using pore dimension and a one end correction (Brown and Escombe, 1900;Cowan, 1978;Weyers and Meidner, 1991). The values from two thermocouples were used to measure the internal temperature of the replica attached on each side and compared with direct measurements from thermal imaging ( Supplementary Fig. S2). Temperature measured using thermal imaging resulted in average temperature kinetics from both sides and was used to fit the energy balance model and predict the conductance. The value of conductance determined by fitting the energy balance model on a step increase and decrease in light intensity was validated on an independent data set. In addition, the value was compared with conductance measured using a Li-Cor 6800 (Li-COR Biosciences, Lincoln, NE, USA) and a 9 cm 2 chamber under steady-state conditions (PPFD, 0 µmol m -2 s -1 ; [CO 2 ], 400 µmol mol -1 ; T air , 20 °C; RH, 60%; flow, 500 µmol s -1 ), and calculated from the anatomy of the transpirating surface (density and size of the pores).

Thermal imaging
Thermal images were recorded using bespoke designed software that continuously records the environmental variables, as well as the temperatures, and performed an emissivity correction for the regions of interest (see below). Black and white references were cut from the same piece of aluminium (width, 2 cm; length, 10 cm; thickness, 0.95 mm) and used as references in each captured picture. The temperature of the black reference was also measured with a surface thermocouple, allowing a comparison of the precision and accuracy of the estimates provided by the thermal camera. A temperature correction was performed using the method described in the manual of the thermal camera in which a crumbled piece of aluminium foil was used to correct each image for the influence of the reflected thermal emission from ambient sources. Images were saved every 3 s and each was an average of 100 raw images on a pixel-by-pixel basis, removing random noise. Leaf isolation was performed by thresholding the leaf area on the captured picture with the highest temperature difference between the air and the black reference (i.e. contrast).

Set-up validation
The thermal camera used in these experiments was a FLIR A655sc (FLIR system AB, Täby, Sweden) including an uncooled microbolometer detector (resolution, 640×480 pixels; spectral range, 7.5-14.0 µm; noise equivalent temperature difference, <30 mK) and autocalibrates with surrounding temperature. The light was provided by two identical LED light sources (LX601C, Heliospectra AB, Göteborg, Sweden) located on each side of the thermal camera. After a step change in light intensity, air temperature increased and induced artefact in temperature measurements using the thermal camera. To compensate for temperature drifts, a process called non-uniformity calibration was performed every minute, which calculated a new table of correction factors using a miniature black body that moved in front of the detector (Olbrycht and Więcek, 2015). Additionally, temperature from a black aluminium reference was measured using the thermal camera and compared with surface thermocouple measurements to check if the thermal kinetics were accurately captured ( Supplementary Fig. S3).

Plant material
Winter wheat (Triticum aestivum L.) plants were grown in 20-well propagators under well-watered conditions in peat-based compost (Levingtons F2S; Everris) for 2 weeks, and were vernalized for 6 weeks in a cold room at 5 °C. Plants were potted in 200 ml pots and moved to a controlledenvironment chamber (Photon Systems Instruments, Brno, Czech Republic) with 300 µmol m -2 s -1 of light intensity (10 h/14 h) provided by white LEDs, temperature of 22/18 °C, and an RH of 50/65% (day/night). The third fully developed leaf of the main tiller was used for measurements.

Thermal kinetics in response to a step change in light intensity
Plants were dark acclimated for 1 h and then thermal pictures were recorded for 10 min. After this period, plants were subjected to a step change in light intensity from 0 µmol m -2 s -1 to 430 µmol m -2 s -1 for 1 h, followed by an opposite step from 430 µmol m -2 s -1 to 0 µmol m -2 s -1 also for 1 h. Photon flux was determined with a quantum sensor (SKP 215; Skye Instruments Ltd, Llandrindod Wells, UK) placed near the reference materials, and was converted to energy using a radiation conversion factor measured with a spectroradiometer (model SR9910-PC, Macam Photometrics Ltd, Livingstone, UK). During the experiments, air temperature and RH (Supplementary Figs S3, S5) were recorded simultaneously with the thermal measurements. To ensure even illumination, leaves were attached to a support (length/width, 13/34 cm) that held them flat and horizontal to the light source ensuring no contact with support material. A mixing fan was installed inside the room to provide good air mixing and rapid thermal exchanges, and to prevent creation of local temperature and humidity gradients. The air mixing also increased boundary layer conductance and improved the reactivity of thermal changes, enabling the more efficient capture of temperature variations with the thermal camera.

Validation of the energy balance equation and experimental approach
Various materials (e.g. white reference and leaf replica) were used to validate the experimental approach by evaluating the error between modelled and observed temperature kinetics. Figure 3 shows how close the temperature kinetics predicted using the black reference were to the observed responses from the different materials subjected to a step increase and decrease in irradiance. The thermal kinetics of the white and black references displayed stable temperature during the initial dark period, followed by an exponential increase after the increase Fig. 2. Schematic of the leaf replica used to validate the energy balance equations and experimental set-up. The black tape had an emissivity of 0.97 and an absorbance of 0.96, the aluminium plate provides support with a rapid thermal response, the felt was used as water storage to maintain transpiration through time, and the porous plastic resulted in a constant conductance. The conductance was only dependent on the size of the pores and their distribution across the surface. Two thermocouples were integrated to measure the temperature of the plate on each side. in light intensity and an exponential decrease when light was returned to zero (Fig. 3A). The energy balance model accurately described the thermal kinetics of the white reference [root mean square error (rmse) 0.05 °C] and the observed difference between the white and black temperature during the illumination period that was due to the difference in energy absorbed by the two references. The temperature kinetics of the dry leaf replica followed the same pattern of temperature variations and was accurately predicted by the energy balance model using the same black reference (Fig. 3B, rmse 0.05 °C). The differences in temperatures values observed here were due to differences in light intensity received by the two objects located at different positions under light (300 µmol m -2 s -1 for the black reference, 430 µmol m -2 s -1 for the dry leaf replica) rather than absorbance as for the references. Transpiration from the wet leaf replica resulted in a large decrease in measured (and predicted) temperature compared with the dry leaf replica; however, a conserved pattern of temperature variations was observed (Fig. 3C). The energy balance model accurately described the thermal kinetics of the wet leaf replica by including energy loss by transpiration (Fig. 3C, rmse 0.05 °C). Parameter values of the energy balance equations were estimated by fitting the predicted temperature kinetics of the wet leaf replica on observed data. Using an independent data set and the same parameter values, the model accurately predicted the wet leaf replica thermal kinetics, validating the model predictions (Fig. 3D, rmse 0.08 °C). The pore conductance of the leaf replica was calculated from the size and distribution of the pores at 0.207 mol m -2 s -1 , and was confirmed by infrared gas analysis at 0.210±0.0008 SD mol m -2 s -1 . Pore conductance derived from thermal kinetics using Bayesian inference had a value of 0.207±0.0008 SD mol m -2 s -1 (Fig.  3C, D), demonstrating the predictive power and accuracy of the energy balance model.

Spatiotemporal variability of thermal kinetics in response to a step change in light intensity
Description of leaf temperature kinetics Large spatiotemporal variation in leaf temperature kinetics was observed between and within leaves of wheat (Fig. 4). Modelled temperatures were derived from the temperature kinetics of a black reference (T black , dashed black line). Observations described the temperature kinetics of (A) a white reference, (B) a dry leaf replica, and (C) and (D) a wet leaf replica. Parameter values were adjusted using Bayesian inference to minimize the error between observed and modelled data in (A), (B), and (C). Parameter values estimated in (C) were used to predict the temperature kinetics in (D) to test model predictions using an independent data set. The root mean square error (rmse) was calculated as an estimator of the quality of the model outputs.
A false-colour pallet with a continuous gradient in colour (black to white; see Fig. 4) was used to highlight the temporal changes in temperature that were observed following a step increase in light intensity. Both leaves examined behave differently, with areas reaching the highest and lowest temperature at different times within the response. A second false-colour scale separating data into discrete (colour) bands illustrates the large spatial heterogeneity within the leaves and revealed that different parts of the leaf behaved differently over time. To assess the spatial variation in temperature kinetics over the leaf, the leaves were divided into four separate areas from the tip to the bottom of the leaf. The selected areas displayed differences in temperature of between ~0.5 °C and 1.5 °C during the experiment, as well as significant differences in the temporal responses (Fig. 5A). During the initial dark period (first 10 min), leaf temperature differences (as high as 0.7 °C) were mainly influenced by the amount of transpiration, whereas the initial temperature increased when the light was switched on (first 5 min), and depended on light intensity and leaf characteristics (i.e. thickness, density, and specific heat capacity). This increase in temperature was only counterbalanced once stomata started to respond to the light, opening and increasing leaf evaporative cooling by transpiration. Towards the end of the light period, leaf temperature was mainly influence by the surrounding environmental conditions, tracking the decrease in air RH and increase in air temperature that together influence the leaf to air vapour pressure deficit (VPD I ; Supplementary  Fig. S4). During the second dark period (last hour), the large differences in leaf temperature were due to the level of leaf transpiration achieved during the light period and the rapidity of stomatal closure.

Description of stomatal conductance kinetics
Using an energy balance model fitted on the previously described leaf temperature kinetics ( Fig. 5A; Supplementary  Fig. S5A), g sw was derived ( Fig. 5B; Supplementary Fig. S5B) precisely for all areas and leaves assessed (Supplementary Fig.   Fig. 4. Time-series of thermal images displaying leaf temperature spatiotemporal differences for two leaves (A and B) subjected to changes in light intensities (grey background, 0 μmol m -2 s -1 ; white background, 430 μmol m -2 s -1 ). Two different colour scales are used to highlight either temperature kinetics or heterogeneity over the leaf surface. Average temperature kinetics for (A) leaf 4 and (B) leaf 5 are visible in Fig. 5 Fig. 5. Spatial and temporal response of (A) leaf temperature (T leaf ) and (B) stomatal conductance (g sw ) to step changes in light intensity. Dark areas represent a period where light intensity was 0 µmol m -2 s -1 and the white area a period where light intensity was 430 µmol m -2 s -1 . Each colour represents a leaf, and each curve represents a different leaf position. S6, rmse 0.069 °C). During the experiment, differences in stomatal behaviour and regulation of transpiration were observed within and between the leaves (Fig. 5B). During the initial dark period, g sw values ranged between 0.001 mol m -2 s -1 and 0.092 mol m -2 s -1 , and at the end of the light period g sw ranged between 0.346 mol m -2 s -1 and 0.733 mol m -2 s -1 , illustrating the large diversity in the regulation of transpiration. During the second dark period, g sw continued to increase for some individuals at the beginning and did not return to its original values, displaying a general decrease, with values ranging from 0.004 mol m -2 s -1 to 0.051 mol m -2 s -1 .
Interpretation of stomatal conductance kinetics using an energy balance model By fitting the energy balance model on the observed leaf temperature kinetics, parameter values describing the temporal response of g sw and the thermal exchanges between the leaf and atmosphere were estimated using Bayesian inference (Fig.  6). Leaf-level boundary layer conductance (g bw ) was estimated at 1.284±0.007SD mol m -2 s -1 and no significant variation between leaves was observed. Boundary layer conductance values were sufficiently high so as not to be the main limiting process for gas diffusion during the experiment. An important physiological parameter to understand the thermal response of the leaf was the amount of energy per unit area required to change its temperature by 1 °K (k) which displayed significant differences but was within a contained range of values (775-1002 J m -2 K -1 ). The leaf shortwave absorbance (α l ) was significantly different between leaves, with values ranging between 0.59 and 0.75, respectively, for leaf 1 and 4. The credible intervals of the estimated parameter were such that it was possible to study the variation of stomatal behaviour not only between leaves but also within individual leaves. Steady-state targets of g sw were defined for each dark/light/ dark period (g 1 , g 2 , and g 3 ) and showed significant variation within and between the leaves (Fig. 6A-C). The initial steadystate g sw values during the first dark period (g 1 ) were significantly different between leaves, with up to 10-fold higher values observed for leaf 4 compared with leaf 1 (Fig. 6A). In the second dark period, values of g 3 showed less variation between individuals, with an average value of ~0.025 mol m -2 s -1 . Spatial differences in g 1 and g 3 of ~0.020 mol m -2 s -1 were observed within individual leaves. During the light period, g 2 values showed significant differences between leaves (Fig. 6B), including a positive gradient (~0.1 mol m -2 s -1 increase) along the leaf lamina (between the tip and base of the leaf) for leaves 2-4 and a negative gradient for leaves 1 and 5 (~0.02 mol m -2 s -1 decrease).
In general, temporal responses of g sw displayed an ~2-fold faster increase (Fig. 6D-E) than decrease (Fig. 6G-H), with significant differences observed between leaves. Initial lag time (ϕ i and ϕ d ) showed greater within-leaf variability than between leaves. Interestingly, leaves with low time constant values for an increase in g sw (τ i ) also showed low time constant values for a decrease in g sw (τ d ). The slow increase or decrease of g sw (S l ) observed at the end of the first exponential response after the light was switched on displayed large differences between leaves (Fig. 6F), with positive values (increase) for leaves 1, 2, and 5, and negative values (decrease) for leaf 4. This behaviour was relatively conserved across the leaf lamina.

Thermal kinetics under rapidly changing environmental conditions
To illustrate the robustness of the proposed model to variation in environmental conditions, five leaves were assessed under fluctuating air temperature and relative humidity ( Supplementary  Fig. S7) to produce complex leaf evaporative demand and temperature kinetics (Fig. 7A). Despite the more complex leaf temperature fluctuations and the fact that the g sw model was developed to consider mainly variation in light intensity, model predictions showed a similar precision ( Supplementary  Fig. S8, rmse 0.098 °C) to the predicted leaf temperature kinetics previously described in Fig. 5 (Supplementary Fig. S6, rmse 0.069 °C). Estimation of g sw showed a rapid increase after the light was switched on, followed by a plateau or a slow decrease, and a relatively slow decrease of g sw when the light was switched off (Fig. 7B). In general, parameter values derived from the energy balance model in Fig. 8 were within the same range as those observed previously in Fig. 6. The steady-state targets for g sw for each period of the kinetics (g 1 , g 2 , and g 3 ) showed no consistent pattern but significant differences between leaves. All leaves presented a slow decrease in g sw after the initial exponential increase, with one leaf exhibiting a 4-fold lower value of s l . Shortwave absorbances were slightly higher than those displayed in Fig. 6, with values ranging from 0.76 to 0.89. Time constants for an increase in g sw (τ i ) were about half the values of those estimated for a decrease in g sw (τ d ), with significant differences between leaves. Lag times (ϕ) for the increase or the decrease of g sw had a similar magnitude in both cases. The amount of energy per unit area required to change the temperature of the material by 1 °K (k) was only significantly different between leaf 2 and 3, with an average for all leaves of 1018 J m -2 K -1 (Fig. 8J).

Parameter importance for dynamic energy balance predictions
Under dynamic environmental conditions, variations in parameter values in the energy balance equations impact leaf temperature kinetics with a different magnitude at different periods throughout the kinetics. The impact of these values on temperature kinetics is displayed in Fig. 8, and illustrates how much parameter values need to be adjusted to achieve a ±0.5 °C temperature variation within the step changes in light intensity. The results illustrated the relative importance of the boundary layer conductance (g bw , Fig. 9A), which influenced the entire temporal response of leaf temperature when changed by about ±50%. Leaf shortwave absorbance (α l ) impacted leaf temperature during the light period but required variations in value that were biologically improbable, as alterations of 50% were required to drive a 0.5 °C change in temperature. Steadystate targets for g sw (g 1 , g 2 , and g 3 ) impacted leaf temperature kinetics during the period for which they are defined as well as the initial part of the following period, with variations that were within the observed range of values shown in Figs 6 and 8. During the light period, increasing the leaf temperature by 0.5 °C required a 3-fold increase in the value of the slow increase or decrease in g sw (s l ). Parameters influencing the temporal response of g sw (Fig. 9G-J) only had a transient impact on leaf temperature but did not require large variations in values. For example, increasing the initial time lag (ϕ i ) or the time constant (τ i ) by ~1 min after a step increase in light intensity resulted in a 0.5 °C increase in leaf temperature, which lasted for 10 min (Fig. 9G, H). During the second dark period, increasing the initial time lag (ϕ d ) or the time constant (τ d ) by ~3 min resulted in a 0.5 °C decrease in leaf temperature, which lasted for at least 30 min (Fig. 9I, J).

Discussion
Despite the success of thermometry in selecting plants with improved yield or altered response to drought (Raskin and Ladyman, 1988;Reynolds et al., 1999;Merlot et al., 2002), a major limitation of this technique is the need for stable environmental conditions to interpret the temperature differences (Reynolds et al., 1999;Jones et al., 2009;Rischbeck et al., 2017;Prado et al., 2018).The results presented here provide strong evidence that thermography can be used to derive g sw under a dynamic environment, opening up a new avenue for plant phenotyping and selection. This is particularly relevant given the growing evidence that temporal responses of g sw limit Fig. 7. Temporal response of (A) leaf temperature (T leaf ) and (B) stomatal conductance (g sw ) to step changes in light intensity under fluctuating environmental conditions (air temperature and relative humidity). Dark areas represent a period where light intensity was 0 µmol m -2 s -1 and the white area a period where light intensity was 430 µmol m -2 s -1 . Each colour represents a leaf, and each curve represents a leaf section.
Our approach to describe the combination of leaf energy balance and mass transfer does not depend on any specific reference material or environmental conditions, making this method a versatile tool for thermography. The equations described in this study have been applied using aluminium references to predict, with high accuracy, leaf temperature kinetics in wheat, allowing g sw to be derived under a fluctuating light environment. Model outputs were validated using a leaf replica with a known conductance to water vapor subjected to the same fluctuating conditions. Although leaf replicas have been used in the past to study the importance of mass transfer in leaf energy balance (Zwieniecki et al., 2016;, our results highlight the potential of the low-cost leaf replica and references to validate leaf temperature predictions obtained using energy balance equations. The ability to derive g sw from thermography under a fluctuating environment using a relatively simple set-up opens the way to field measurements in the future, but will require further improvements to take into consideration parameters such as leaf orientations and local variations in L d . In this context, our approach could be combined with hemispherical references placed at different positions in the field (Jones et al., 2018) to take into account the different leaf angles and surrounding thermal conditions present in a field canopy. Indeed, our model can be easily adapted to an existing experimental set-up or used to reprocess existing data if enough information is available. Previous work has attempted to use dynamic energy balance equations (taking into consideration the temporal dimension) but are often incomplete in their implementation or use approximations to describe the continuous interactions between the leaf and a fluctuating environment (Bajons et al., 2005;Schymanski et al., 2013;Page et al., 2018). The ability of the model presented here to operate under any environmental condition is mainly due to the unique combination of equations integrating the continuous variations of the environmental variables, the impact of stomatal behaviour, and ultimately leaf temperature.
Combining stomatal conductance and energy balance models opens up new possibilities to use leaf temperature kinetics to improve our understanding of the mechanisms involved in stomatal responses to the surrounding environment. However, combining both models requires an estimate of boundary layer conductance (g bw ) that dictates part of the gas and heat exchange. Therefore, our model simultaneously estimates g sw , g bw , and leaf temperature in fluctuating environmental conditions using a set-up that is easy to use and low cost, while providing precise estimates. The values of g bw used to calculate g sw from the energy balance equations produced g sw values that were close to those of the leaf replica, suggesting that our Fig. 9. Sensitivity analysis representing the variation of parameter values required to change leaf temperature by ±0.5 °C during step changes of light intensity. Parameter values from leaf 2 and area 1 were used as an illustration (Fig. 5). Dark areas represent a period where light intensity was 0 µmol m -2 s -1 and the white area a period where light intensity was 430 µmol m -2 s -1 . Differences in leaf temperature were only achieved over parts of the temperature kinetics depending on the parameter (red shaded area) and reached during these periods a maximum of ±0.5 °C. Original parameter values and corresponding curves (solid black line) are displayed in each plot, with the achieved temperature differences (T diff ) and the corresponding parameter values (θ). In some cases (C, E, F, G, I), the maximum temperature deviation of ±0.5 °C was not reached due to a parameter value reaching a boundary (e.g. 0); the value achieved was shown instead of ±0.5 °C.
estimates of g bw were accurate enough to derive a valid g sw value. Numerous research studies has been published on measuring leaf boundary layer conductance (Leuning et al., 1989;Leuning, 1990;Brenner and Jarvis, 1995;Stokes et al., 2006;Katsoulas et al., 2007) and although our method is derived from a similar theoretical basis, it does not require specific equipment (e.g. power source, heating element) and relies on environmental variations, making it simple to operate. In the case of portable gas exchange chambers, investigating stomatal responses requires artificial conditions that alter the thermal exchange between leaf and atmosphere (e.g. influence of wall temperature) and often use high boundary layer conductance to simplify the measurements. Using thermography enables the determination of g sw without artificially altering the conditions surrounding the leaf, meaning that measurements of leaf responses are closer to those observed in the field. The outcomes from our modelling approach included the temporal response of g sw , a biological trait that has attracted significant attention in recent years (Lawson et al., 2010;Lawson and Blatt, 2014;McAusland et al., 2016;Meinzer et al., 2017;Matthews et al., 2018;Deans et al., 2018) because of its impact on both photosynthesis and water use efficiency.
Nocturnal stomatal conductance (observed after dark acclimation) was significantly different between individuals, suggesting differences in the regulation of water loss during the night period. Nocturnal transpiration is involved in essential physiological processes such as nutrient transport (Zeppel et al., 2014) and could represent 30% of daytime water consumption in crops (Claverie et al., 2018). The initial g sw values observed here under dark conditions represented up to ~10% of the value reached under the light period, confirming its importance for the regulation of plant water budget. Nocturnal regulation of g sw has been shown to be developmentally and genetically controlled in wheat Schoppach et al., 2016) and has been proposed as a breeding target to produce plants with improved tolerance to drought (Schoppach et al., 2014Coupel-Ledru et al., 2016). These results highlight the potential of our method for future breeding programmes aiming to screen populations with differences in nocturnal transpiration.
Under the light period, g sw values were similar to those reported previously in the literature for wheat (McAusland et al., 2016;Taylor and Long, 2017). However, stomatal responses reported here were more rapid than those reported for wheat by McAusland et al. (2016) and displayed a slow decrease in g sw after the initial exponential response, which was also significantly different between leaves. These two independent estimations of the temporal response of g sw for wheat were performed under different environmental conditions (controlled versus uncontrolled) and we believe that the rapid increase of leaf temperature after the light was switched on (due to the increase in incident energy), simultaneously with a decrease in air RH and an increase in air temperature experienced by the whole plant in the imaging area, may have contributed to this change in stomatal behaviour. Interestingly, Franks and Farquhar (2007) observed similar enhancement in wheat of the initial rate of increase in g sw and a lower final value at high VPD I (2 kPa) compared with low VPD I (1 kPa). Mott et al. (1999) proposed that the difference in response depends on hydraulic interactions among stomata that are mediated by transpiration-induced changes in epidermal turgor, which could explain our observations. In addition, after a step increase in light intensity, we observed an initial lag (~2-3 min) before g sw changed significantly, confirming the observations in wheat made by Taylor and Long (2017). This initial lag has been shown to impact photosynthesis in different species (Vialet-Chabrand et al., 2013;McAusland et al., 2016;Taylor and Long, 2017) but also leads to a rapid increase in leaf temperature due to the limited cooling by transpiration. Schymanski et al. (2013) highlighted the potential heat damage and hydraulic failure a leaf may experience during the first minutes of a sun fleck, which could be avoided by maintaining g sw at a relatively high value when the leaf returns to shade conditions. During our experiment, the lag time and rapidity of response for a decrease in g sw were substantially larger than for an increase in g sw , ensuring that g sw remained at a high value for a longer time period. Additionally, the sensitivity analysis performed on the model revealed that a small increase in parameter values ϕ i/d and τ i/d , controlling the rapidity of g sw variation, resulted in maintaining a low leaf temperature for a greater duration after a decrease in light intensity. These observations are compatible with a conservative strategy to maximize carbon fixation by unit of water loss (Ooba and Takahashi, 2003) and limit heat damage (Schymanski et al., 2013) under a fluctuating light environment. Even if there is asymmetry in the response between the rapidity of increase and decrease in g sw , it is interesting to note that leaves with a slow increase in g sw also had a slow decrease, suggesting a co-ordination of both responses as proposed by McAusland et al. (2016). Overall, the diversity observed here for only a few individuals is a promising result that could be used as breeding targets to improve yield by reducing the stomatal limitation on photosynthesis as suggested by Lawson and Blatt (2014), as well as to enhance thermal regulation and water-saving strategies to limit heat damage and increase drought tolerance in hot and dry environments.
Under a dynamic environment, thermograms revealed variation in temperature kinetics over the leaf lamina that could be attributed to patchy stomatal behaviour (Weyers and Lawson, 1997;Lawson and Weyers, 1999). Compared with previous research studying the spatial heterogeneity of g sw using thermograms (Jones, 1999;Saudreau et al., 2017;Sweet et al., 2017;Page et al., 2018), we assessed the temporal response of g sw for different sections of wheat leaves following a step increase in light intensity to characterize the rapidity and the magnitude of the stomatal response. The spatial gradient observed for the magnitude of the g sw response (up to an ~50% increase from the base to the tip of the leaf) stressed the importance of measuring g sw at a similar leaf position for each individual when using a portable gas exchange chamber (Lawson and Weyers, 1999). Estimations of the rapidity of the stomatal response are less influenced by the position on the leaf blade, suggesting that at this scale in wheat, any patchy effect is mainly due to variation of stomatal density rather than different stomatal behaviour. Prytz et al. (2003) monitored spatial and temporal stomatal conductance of Avena subject to low humidity and suggested that the synchronicity of stomatal behaviour at different leaf positions may be a result of the anatomy of monocotyledonous leaves in which the main veins running parallel with the leaf blade transport a hydraulic signal synchronizing different areas of the leaf in the longitudinal direction. More specifically, in wheat, Buckley and Mott (2000) provided strong evidence for long-distance hydraulic interactions coordinating the stomatal response in different areas of the leaf that could explain our observations. Previous research using infrared gas exchange to measure stomatal behaviour has shown that the temporal response of g sw displayed significant variations within (Matthews et al., 2018) and between species (McAusland et al., 2016) with comparable responses to our results. Previous studies using thermal signatures to assess stomatal kinetics have been limited to the first few minutes of a response to a step change in light intensity (Bajons et al., 2005;Page et al., 2018) or have only interpreted the relative changes in temperature (Prytz et al., 2003). Overall, the model presented here allowed us to investigate the parameters controlling the temporal response of g sw (τ, ϕ, S l ) that have been related to guard cell metabolism (Hills et al., 2012;Wang et al., 2014;Vialet-Chabrand et al., 2017a), and provided strong evidence of a co-ordinated rapidity of response across the leaf lamina in wheat despite the local variation in the magnitude of the response.
The findings and experimental approaches presented here have the potential to remove a major bottleneck in highthroughput phenotyping of stomatal-related traits, by allowing the interpretation of thermograms under a fluctuating environment. Time-series of thermograms associated with our new energy balance approach provided spatial and temporal characterization of stomatal conductance (g sw ) responses in wheat, highlighting the importance of co-ordinated stomatal responses across the leaf blade. The diversity and asymmetry of the temporal response of g sw observed after a step increase or decrease in light intensity can be interpreted as a strategy to maximize photosynthesis per unit of water loss and avoid heat stress under a fluctuating environment. Improving these aspects of stomatal behaviour is thought to be an important breeding target for future yield improvement, to which our work will directly contribute by allowing screening of large numbers of plants for biologically relevant traits. The techniques and methods presented here provide significant evidence for the accuracy of predicted g sw from leaf temperature kinetics, that can be easily transferred to existing equipment and will pave the way for further development in our understanding of stomatal behaviour in future field studies.

Supplementary data
Supplementary data are available at JXB online. Table S1. Equations required to derive environmental variables from the raw environmental measurements collected from various sensors. Fig. S1. Example of signal processing to remove high-frequency noise from infrared thermal measurement. Fig. S2. Example of temperature kinetics measured using thermal imaging. Fig. S3. Comparison of temperature measurements using an infrared thermal camera and a thermocouple. Fig. S4. Environmental conditions during step changes of light intensity. Fig. S5. Example of the performance of the energy balance model to reproduce leaf temperature kinetics and stomatal conductance. Fig. S6. Performance of the energy balance model to reproduce leaf temperature kinetics. Fig. S7. Environmental conditions during step changes of light intensity. Fig. S8. Performance of the energy balance model to reproduce leaf temperature kinetics.
Algorithm S1. Equations describing the energy balance model.