tealeaves: an R package for modelling leaf temperature using energy budgets

Abstract Plants must regulate leaf temperature to optimize photosynthesis, control water loss and prevent damage caused by overheating or freezing. Physical models of leaf energy budgets calculate the energy fluxes and leaf temperatures for a given set leaf and environmental parameters. These models can provide deep insight into the variation in leaf form and function, but there are few computational tools available to use these models. Here I introduce a new R package called tealeaves to make complex leaf energy budget models accessible to a broader array of plant scientists. This package enables novice users to start modelling leaf energy budgets quickly while allowing experts to customize their parameter settings. The code is open source, freely available and readily integrates with other R tools for scientific computing. This paper describes the current functionality of tealeaves, but new features will be added in future releases. This software tool will advance new research on leaf thermal physiology to advance our understanding of basic and applied plant science.


Introduction
Plants grow, survive and reproduce under a wide variety of temperatures because natural selection endows them with adaptations to cope with different thermal regimes. Cushion plants in the alpine grow near the ground to stay warm, desert plants decrease absorptance to stay cool (Ehleringer et al. 1976) and plants keep stomata open, which can protect against extreme heat waves (Drake et al. 2018). Understanding these diverse mechanisms of thermal adaptation and acclimation may provide insight into how plants respond to increasing temperatures and how these responses influence ecosystem function with anthropogenic climate change (Rogers et al. 2017;Crous 2019). Because leaves are the primary photosynthetic organ in most plants, regulating leaf temperature is critical (Berry and Björkman 1980). Photosynthesis peaks at intermediate temperatures (Sage and Kubien 2007). When leaves are too warm, evaporation increases exponentially, photo-and nonphotorespiratory losses subtract from carbon gain (Jones 2014) and critical loss of function occurs about ~50 °C (O'Sullivan et al. 2017). When leaves are too cold, maximum photosynthetic rates decline and can lead to damage from excess solar radiation (Huner et al. 1993) as well as nighttime dew and frost formation (Jordan and Smith 1994). Natural selection should favor leaf morphologies and physiological responses that optimize leaf temperature in a given environment (Parkhurst and Loucks 1972;Okajima et al. 2012;Michaletz et al. 2016).
To understand leaf thermal physiology, plant scientists need mathematical and computational tools to model leaf temperature as a function of leaf traits and the environment. Balancing energy budgets is a powerful mathematical tool for understanding how leaf traits and environmental parameters influnce plant physiology that has been used for over a century (Raschke 1960). The equilibrium leaf temperature is that in which the energy gained from incoming solar and infrared radiation is balanced by thermal infrared radiation losses, sensible heat loss/gain and latent heat loss through evaporation (Gutschick 2016). Leaf angle, size and conductance to water vapour alter leaf temperature by changing how much solar radiation they intercept and how much heat they lose through sensible and latent heat flux. Likewise, enrvironmental factors such as sunlight, air temperature, humidity and wind speed influence heat transfer between leaves and the surrounding microclimate (Gutschick 2016). Under controlled conditions, leaf energy budget models are highly accurate (Schymanski and Or 2017). Hence, they can offer deep insight on plant thermal physiology by asking how temperature is affected by one factor in isolation or in combination with another.
Leaf energy budget models have many applications, but perhaps their most widespread use is in modelling optimal leaf size and shape. The boundary layer of still air just above and below the leaf surface determines sensible and latent heat transfer and is proportional to leaf size (Gates 1968). All else being equal, larger leaves have a thicker boundary layer, slowing heat transfer and decoupling leaf temperature from air temperature. This likely explains why, for example, many warm desert species have small leaves (Gibson 1998). Using leaf energy budgets, Parkhurst and Loucks (1972) further predicted that leaves should be small in cold air and large under warm, shaded conditions. More recently, Okajima et al. (2012) extended these models, showing that small leaves maximize photosynthetic rate under high insolation and warm temperatures, but large leaves increase water-use efficiency in shadier habitats. Wright et al. (2017) used energy budget models to show that dew and frost formation may select against large leaves at high latitudes. Energy budget models also help explain variation in leaf shape, such as lobing and dissection, because heat transfer is determined by effective leaf width (aka characteristic leaf dimension; Taylor 1975) rather than total area. Effective leaf width is 'the diameter of the largest circle that can be inscribed within the margin' (Leigh et al. 2017). Lower effective leaf width reduces leaf temperature under natural conditions in the sun (Leigh et al. 2017) and is under selection in sunny, drier habitats (Ferris et al. 2015). Besides leaf size and shape, energy balance models are useful in understanding many plant processes and traits (Gates 1965), such as transpiration (Gates 1968), optimal stomatal conductance (Buckley et al. 2014), stomatal arrangements (Foster and Smith 1986), leaf thickness (Leigh et al. 2012), response to sunflecks (Schymanski et al. 2013), carbon economics (Michaletz et al. 2016) and water-use efficiency (Schymanski and Or 2016).
Despite the utility of leaf energy budget models, there are a dearth of open source, customizable, computational tools to implement them. The plantecophys package implements a similar energy budget model (Duursma 2015). However, the model is simplified for faster computation needed in ecosystem and global land surface models (Leuning et al. 1995). Therefore, it does not incorporate features such as different boundary layer conductances on each leaf surface, nor can users easily change default parameters for specialized cases. The Landflux website also has an Excel spreadsheet for leaf energy budgets (Tu 2019), but it is prohibitively time-consuming and not reproducible to use spreadsheets for large-scale simulations. Because computational tools are limited, potential users must develop models anew and learn the numerical methods necessary to find solutions. Ideally, there should be a platform in which novices can model leaf temperature to solve an interesting problem without having to write their own model and learn complicated numerical algorithms. At the same time, we need a platform that can be easily modified for experts that want to extend existing leaf energy balance models.
The goal of this paper is therefore to develop software that models leaf temperature as a function of leaf traits and the environment with physical realism. This software should be open source so that the methods are transparent and code can be modified by other researchers. Secondly, it should be readily available to novice modelers yet customizable by those working on more specific problems. Finally, it should easily integrate with other advanced tools for scientific computing. To that end, I developed an R package called tealeaves to model leaf temepature in response to a wide variety of leaf and environmental parameters. The source code is open source and available to modify; it is easy to use with default parameters, but also customizable; and because it is written in R, the output from tealeaves can be analysed and visualized with the vast array of computational tools availble in the R environment.
Leaf energy budgets consist of incoming and outgoing energy fluxes. Incoming energy includes radiation from solar (aka shortwave) and thermal infrared (aka long-wave) sources. Outgoing energy includes losses of infrared radiation, sensible heat and latent heat flux due to evaporation (Fig. 1). Note that sensible and latent fluxes can be positive when leaves are warmer than the air or condensation (dew) occures, respectively. When leaves reach a thermal equilibrium with their environment-generally within a few minutes-these incoming and outgoing energy sources balance one another. Formally, one solves for the leaf temperature (T leaf ) that balances the energy budget: (1) R abs is the absorbed radiation, S r is thermal infrared radiation loss, H is sensible heat flux and L is latent heat flux. All of these values have units W m −2 . Environmental and leaf parameters like irradiance, air temperature and leaf absorbtivity determine how much energy is absorbed and radiated. In the daytime, radiation coming into the leaf generally exceeds that coming out, resulting in a leaf temperature above air temperature. Sensible heat loss (gain) occurs when the leaf is warmer (cooler) than the surrounding air, but the rate is influenced by other parameters like wind speed and leaf size. Latent heat loss occurs when the water vapor pressure of the surrounding air is lower than that in the leaf, driving evaporation. Relative humidity of the air and leaf conductance to water vapour determine how strong the vapor-pressure deficit is and how much resistance to water transport occurs. Leaves gain energy when water condenses as dew. The leaf energy balance model works by taking a set of environmental and leaf parameters, then finding the leaf temperature at which they balance one another. The equilibrium leaf temperature can be well above air temperature under high irradiance, low wind speed and/or low conductance to water vapor; leaf temperature can be near or even below air temperature under low irradiance, in small leaves and/or leaves with high rates of evaporation. The primary aim of this paper is not to extend leaf energy budget theory, but to describe the tealeaves package which implements existing models. Therefore, the mathematical details behind the model are provided in an Appendix at the end of this paper. Tables 1 and 2 list all mathematical symbols in parameter inputs and calculated output values. Supporting Information- Table S1 lists current default parameter values and realistic ranges with references to the literature. These equations describe the current tealeaves implementation. As 2) The leaf radiates some energy. 3) Heat is lost through convection when the leaf temperature is greater than air temperature. 4) Latent heat is lost through evaporation, which is driven by the vapor-pressure deficit between the air and leaf interior and the leaf conductance to water vapor. Calculations used the following leaf parameter values in this example: d = 0.1 m; Calculations used the following environmental parameter values in this example: Table 1 for symbol definitions. future releases will alter some assumptions and incorporate new features, I mention future modifications in the Discussion. In this section, I describe how leaf energy budget models are implemented in R and provide worked examples.

Solving in R
R is a fully open source programming language for statistical computing that allows users to develop their own packages with new functions. tealeaves takes three sets of parameter inputs: leaf parameters, environmental parameters and physical constants (see Table 1). The package provides reasonable defaults, but users can input new values to address their question, as I demonstrate in the next section. With one or more parameter sets, tealeaves uses the uniroot function in R base package stats to find the T leaf that balances the leaf energy budget (Equation 1). It outputs the equilibrium T leaf and energy fluxes in a table for analysis and visualization.
Unlike previous leaf energy models, tealeaves ensures that calculations are technically correct by assigning stadard SI units with the R package units Pebesma et al. (2016). Every parameter and calculated value must have correctly assigned units. If units are not properly defined, tealeaves will produce an error because it is unable to convert values. For speed, calculations can optionally be made without units and will yield correct results if provided values have correct units. To ensure accuracy, these unitless functions are tested against their counterparts with units using the testthat package (Wickham 2011). Other R packages that contributed to tealeaves are crayon (Csárdi 2017), dplyr , glue (Hester 2019), furrr (Vaughan and Dancho 2018), future (Bengtsson 2019), ggplot2 (Wickham 2016), magrittr (Bache and Wickham 2014), purrr (Henry and Wickham 2019a), rlang (Henry and Wickham 2019b), stringr (Wickham 2019), tidyr .

Worked examples
In this section, I provide two worked examples; more complex worked examples are found in the Supporting Information. The first illustrates that it is straightforward to use tealeaves with a few lines of code with default settings. The second shows that it is also possible to model T leaf across multiple leaf trait and environmental gradients for more advanced applications.

Example 1: a minimum worked example.
The box below provides R code implementing the minimum worked example with default settings.

Extended examples
To see the range of possible applications for tealeaves, I ran four additional sets of simulations. The first models the leaf-toair temperature differential for different leaves sizes across a gradient of air tempuratures; the second models the leaf-to-air temperature differential across a gradient of incident solar radiation for different stomatal conductances; the third models the leaf-to-air temperature differential for different-sized leaves under free, mixed and forced convection; and the fourth models the effect of stomatal ratio on evaporation under free and forced convection. These extended examples are documented more fully in the Supporting Information with accompanying R code.
To provide a sense of which leaf and environmental parameters affect T leaf the most under 'typical' conditions, I varied stomatal conductance (g sw ), leaf size (d), stomatal ratio (SR), relative humidity (RH), solar radiation (S sw ) and wind speed (u) over a wide range of realistic values while holding all other values constant at their default setting [see Supporting Information- Table S1].

tealeaves's source code is open to all
A development version of tealeaves is currently available on GitHub (https://github.com/cdmuir/tealeaves). A stable version of tealeaves will be released on the Comprehensive R Archive Network (CRAN, https://CRAN.R-project.org/package=tealeaves). I will continue developing the package and depositing revised source code on GitHub between stable release versions. Other plant scientists can contribute code to improve tealeaves or modify the source code on their own installations for a more fully customized implementation.  Fig. S1B]. Leaves can operate below air temperature at low light, but above it at higher light [see Supporting Information- Fig. S1E]. Stomatal ratio has only a modest effect on leaf temperature [see Supporting Information- Fig. S1C]. Most users will want to modify these default settings and simulate leaf temperature over a range of leaf and environmental parameters, so these results are not generalizable to all cases. By design, tealeaves easily allows users to define multiple simultaneous trait and environmental gradients (see Worked Example 2 and Fig. 2).

Discussion
Scientists have used energy budgets to model leaf temperature for over a century (see Raschke (1960) for historical references). Despite many advances in our understanding of the environmental and leaf parameters that affect heat exchange (Gutschick 2016), there exist few computational tools to implement complex energy budget models. The tealeaves package fills this gap by providing a platform for modelling energy budgets in a transparent and reproducible way with R (R Core Team 2019), a freely available and widely used programming language for scientific computing. Unlike previous software, tealeaves removes ambiguity by forcing users to specify proper SI units through the R package units (Pebesma et al. 2016). Neophytes with little experience modelling leaf temperature may get started quickly without having to develop their model de novo, while specialists can modify the code to customize tealeaves to their specificiations. tealeaves also readily integrates with the vast array of data analysis and visualization tools in R. These features will enable wider adoption of leaf energy budgets models to understand plant biology. However, as I discuss below, the current version of tealeaves has several important limitiations that can be addressed in future releases. Previously, researchers wanting to implement sophisticated leaf energy budget models that required numerical solutions had to write their model and learn a numerical algorithm to solve it. Most often, these solutions are not published and/or are not open source. This slows down research for non-specialists by introducing unnecessary barriers and can be error-prone. For example, the current tealeaves model relies on previous work by Foster and Smith (1986). Without a platform like tealeaves, extending their work required developing the mathematical and computational tools de novo every time. Also, the published version of Foster and Smith (1986) contains several small errors and typographical inconsistencies in the equations. While these are most likely mistakes made during typesetting and publication, without open source code, it is very challenging to determine if these mistakes also occurred in their computer simulations. Transparent, open source code does not prevent mistakes, but makes it easier for the community to discover mistakes and fix them faster.
The tealeaves model is more complex than most other leaf energy balance models in the literature, which makes it more flexible and realistic but more computationally intensive. For example, the R package plantecophys (Duursma 2015) uses the isothermal net radiation approximation of Leuning et al. (1995). This reduces the number of iterations and speeds up computation by assuming that long-wave radiative flux from the leaf is proportional to air rather than leaf temperature. Other recent models (e.g. Buckley et al. 2014;Schymanski and Or 2017) assume free convection or do not account for virtual temperature differences in moist air (Equation 14) (Okajima et al. 2012;Duursma 2015). Other models, except Okajima et al. 2012, also do not account for different surrounding temperatures above and below the leaf surface (T sky and T air , respectively) when calculating thermal infrared radiation (see Equation 3). These simplifications may be adequate for many applications and the benefit in speed may be necessary for scaling up. However, recent comparisons of simple and complex energy balance models in the field found that the complex models, like that implemented in tealeaves, fit the observed distribution of wheat leaf canopy temperatures better (Webber et al. 2017). Future development of tealeaves will make it easier to implement simpler models. A full range of models, from simple to complex, will serve the broadest range of applications and allow users to test when simplifying assumptions are justified. To facilitate this, future versions of tealeaves will allow users to provide functions rather than values for some parameters. For example, a less computationally intensive model could be contructed by using a function for calculating boundary layer conductances that assume free convection is negligible. Or users could substitute a more appropriate function for their application, such as cloudy rather than clear skies in the case of Equation 3.
Ultimately, the goal of tealeaves is to provide a platform for implementing realistic and fully customizable energy budget models. Such models may take too much computational time to be useful for large-scale ecosystem models, but they can help understand a wider range of fascinating and poorly understood leaf anatomical and morphological features. The Introduction lists several possible uses, but most of these problems cannot currently be solved with tealeaves alone. For example, many photosynthetic processes are temperature-sensitive, but it would require simultaneous modelling of leaf temperature, stomatal conductance and photosynthesis to predict optimal trait values. tealeaves intentionally does not specify these other models because the concept is that it should stand alone and be able to interact with many other models. For example, tealeaves could be combined with a stomatal response function such as the Ball- Berry model (Ball et al. 1987) or found through optimization (Buckley et al. 2014;Duursma 2015;Muir 2019). Similarly, the leaf temperature throughout a canopy could be modeled by running tealeaves over a gradient of light, wind, temperature and so forth. tealeaves should therefore be thought of as one component in an expanding ecosystem of interrelated tools for modelling plant physiology.
Currently, tealeaves has several limitations that I plan to address in future releases. It uses rather simple models of infrared radiation and direct versus diffuse radiation. Ideally, it would be better if users could supply their own functions to calculate these parameters from the total irradiance. The model also assumes leaves are horizontal, whereas leaf orientation varies widely. Following previous authors, I modeled heat transfer as a mixed convection (Equations 9 and 27), but this may not adequately describe real leaf heat exchange (Roth-Nebelsick 2001). tealeaves calculates equilibrium as opposed to transient behavior (Vialet-Chabrand and Lawson 2019), which may takes several minutes to reach. Finally, the model assumes a single homogenous leaf temperature rather than using finite element modelling to calculate leaf temperature gradients across leaves of different shapes. These are important limitations of the current software which can be addressed in future work.
In conclusion, tealeaves provides new software for leaf energy balance models in R. Leaf energy balance models are highly useful tools for understanding plant form and function and new computational tools will make these models more broadly accessible, advancing basic and applied plant science.

Supporting Information
The following additional information is available in the online version of this article- Table S1. Reasonable values for tealeaves parameter inputs with references to the primary literature. The current version of tealeaves uses a default value within the range of reasonable values. See Table 1 for a key to symbols. The temperature of smaller leaves is more closely coupled to air temperature. Each line represents a different leaf size (small, dashed line; medium, dotted line; large, solid line) and the leaf-to-air temperature differential (y-axis) over an air temperature gradient (x-axis). (B) Greater stomatal conductance cools leaves. Each line represents a different stomatal conductance (low, dashed line, 1 µmol m −2 s −1 Pa −1 ; medium, dotted line, 3 µmol m −2 s −1 Pa −1 ; high, solid line, 5 µmol m −2 s −1 Pa −1 ) and the leaf-to-air temperature differential (y-axis) over a gradient of incident solar radiation (x-axis). (C) Forced convection dominates in small leaves; free convection dominates in very large leaves. Leaf size is indicated by line type as in Panel (A). Vertical lines indicate approximate shifts from forced convection (Ar < 0.1), mixed convection (0.1 < Ar < 10) and free convection (Ar > 10).
Small leaves always experience forced convection, leading to lower leaf temperature compared to large leaves experiencing free convection. (D) Amphistomatous leaves (stomatal ratio ~ 0.5) evaporate more than hypo-or hyperstomatous leaves (stomatal ratio ~ 0 or 1, respectively), especially under free convection (low wind speed, u).
(SR = 0.5, x-axis) lose more water through transpiration than leaves with all stomata on one surface, leading to a lower leaf temperature. (D) Greater humidity (RH, x-axis) increases leaf temperature by limiting latent heat loss. (E) With low solar radiation (S sw , x-axis), leaf temperature is below air temperature; with high solar radiation, leaf temperature is greater than air temperature. (F) At greater wind speeds (u, x-axis), leaf temperature is more closely coupled to air temperature. The discontinuity represents the shift from laminar to turbulent flow. For reference, the dashed line is the air temperature in all simulations. All calculations used the following leaf parameter values unless they varied: d = 0.1 m; α s = 0.5; α l = 0.97; g sw = 5 µmol m −2 s −1 Pa −1 ; g uw = 0.1 µmol m −2 s −1 Pa −1 ; SR = 0.5 . All calculations used the following environmental parameter values unless they varied: P = 101.3246 kPa; r = 0.2; RH = 0.5; S sw = 1000 W m −2 ; T air = 298.15 K; u = 2 m s −1 . See Table 1 for symbol definitions.

Sources of Funding
This study received startup funds from the University of Hawai'i.

Conflict of Interest
None declared.
where Nu forced = aRe b (10) a, b, c, d are constants that depend on whether flow is laminar or turbulent and the direction of flow in the case of free convection (see below). In general, when the Archimedes number (also called the Richardson number) Ar = Gr/Re 2 0.1 , free convection dominates; when Ar = Gr/Re 2 10, forced convection dominates (Nobel 2009). The Nusselt number coefficients can be found in Monteith and Unsworth (2013). For forced convection, flow is laminar if Re < 4000, a = 0.6, b = 0.5 ; flow is turbulent if Re > 4000, a = 0.032, b = 0.8. These cutoffs for leaves are lower than for artificial surfaces because trichomes and other anatomical features of leaf surfaces induce turbulence more readily (Grace and Wilson 1976). For free convection, flow is laminar. For the upper surface when T leaf > T air or the lower surface when T leaf < T air , c = 0.5, d = 0.25 . Conversely, for the lower surface when T leaf > T air or the upper surface when T leaf < T air , c = 0.23, d = 0.25.
Grashof and Reynolds numbers approximate the ratio of bouyant or inertial (numerator in equations below), respectively, to viscous forces (denominator in equations below). They are calculated as: G is the gravitational constant and u is wind speed. The diffusion coefficient for momentum in air (D m ) is calculated for a given temperature following the same procedure above for heat diffusion (D h ; see Equation 8). The virtual temperature is calculated according to Monteith and Unsworth (2013) assuming that the leaf airspace is fully saturated while the air is has a vapour-pressure deficit of p air : T v,leaf = T leaf /(1 − (1 − ε)(p sat /P)).
ε is the ratio of water to air molar masses. The saturation water vapour pressure p sat as a function of temperature is calculated using the Goff-Gratch equation (Vömel 2016). The vapour pressure of air is calculated from the relative humidity (RH) as p air = RHp sat .

Latent heat flux and evaporation
Latent heat flux is the product of the latent heat of vaporization, the total leaf conductance to water vapour and the water vapour gradient: The latent heat of vapourization (h vap ) is a linear function of temperature. tealeaves calculates h vap using parameters estimated from linear regression on data from Nobel (2009) The water vapour pressure differential (d wv ) from the inside to the outside of the leaf is the water vapor pressure inside the leaf, which is assumed to be saturated (p leaf = p sat ), minus the water vapor pressure of the air (p air ), calculated as described above. This value is converted from kPa to mol m −3 using the ideal gas law: d wv = p leaf /(RT leaf ) − RHp air /(RT air ). (18) R is the ideal gas constant. The total conductance to water vapor (g tw ) is the sum of the parallel lower (usually abaxial) and upper (usually adaxial) conductances g tw = g w,lower + g w,upper .
The conductance to water vapor on each surface is a function of parallel stomatal (g sw ) and cuticular (g uw ) conductances in series with the boundary layer conductance (g bw ). The stomatal, cuticular and boundary layer conductance on the lower surface are: g uw,lower = (g uw /2)[R(T leaf + T air )/2].
Note that the user provides the total leaf stomatal and cuticular conductance to water vapur in units of µmol m −2 s −1 Pa −1 , which are then converted to units of m s −1 using the ideal gas law. Stomatal conductance is partitioned among leaf surfaces depending on stomatal ratio (SR); cuticular conductance is assumed equal on each leaf surface. The corresponding expressions for the upper surface are: g sw,upper = (g sw SR)[R(T leaf + T air )/2] g uw,upper = g uw,lower .
Partitioning total stomatal and cuticular conductance to water vapor across surfaces is not typically part of leaf energy budget models. However, stomatal (Muir 2015) and cuticular (Karbulková et al. 2008) conductance can differ dramatically on each surface, so it could be useful modify these parameters, especially if leaf temperature models are integrated with photosynthetic models (Muir 2019).
The boundary layer conductances for each surface differ because of free convection (Foster and Smith 1986) and are calculated similarly to that for heat (Equation 7): The diffusion coefficient for water vapor in air at a given temperature (D w ) is calculated using the Equation 8, except that is D w,0 is substituted for D h,0 . Each surface has its own Sherwood number (Sh): As with Nu, Sh is calculated assuming mixed convection: Evaporation rate (mol H 2 O m −2 s −1 ) is the product of the total conductance to water vapour (Equation 19) and the water vapour gradient (Equation 18):